2014-10-27 30 views
1

我有一个函数,它是由三个(k)因子组成的乘积。每个因子是用随机变量R和L减去两个高斯CDF。这些随机变量是根据4个参数定义的。无法集成两个高斯CDF的减法[Wolfram Mathematica]

enter image description here

下面的代码显示了如何绘制的主要功能(根据两个独立变量d和e),以及如何随机变量被计算

sigma = 1; 
k = 3; 
priors = {}; 
AppendTo[priors, 1/k + e]; 
Do[AppendTo[priors, 1/k - e/(k - 1)], {c, 2, k}]; 

L[priors_, sigma_, d_, i_] := Do[ 
maxVal = -Infinity; 
Do[ 
    val = (2*sigma^2*Log[priors[[i]]/priors[[j]]] + d^2 (j^2 - i^2 + 2 (i - j)))/(2 (j - i) d); 
    If[val > maxVal, maxVal = val, Null]; 
, {j, 1, i - 1}]; 
Return[maxVal]; 
, {1}]; 

R[priors_, sigma_, d_, i_] := Do[ 
minVal = Infinity; 
Do[ 
val = (2*sigma^2*Log[priors[[j]]/priors[[i]]] + d^2 (i^2 - j^2 + 2 (j - i)))/(2 (i - j) d); 
If[val < minVal, minVal = val, Null]; 
, {j, i + 1, k}]; 
Return[minVal]; 
, {1}]; 

Print[ 
Plot3D[ 
    Product[ 
    If[R[priors, sigma, d, c] < L[priors, sigma, d, c], 0, 
    (CDF[NormalDistribution[(c - 1) d, sigma], R[priors, sigma, d, c]] - 
    CDF[NormalDistribution[(c - 1) d, sigma], L[priors, sigma, d, c]])] 
    , {c, 1, k}] 
, {d, 0.01, 5} 
, {e, -1/k, 1 - 1/k}, PlotRange -> {All, All, All}, AxesLabel -> Automatic]]; 

Resulting plot of the first code

现在,我想在d上集成函数(在与Plot3D相同的区域中,d = 0.01-5),并根据自变量e绘制结果。

enter image description here

下面是我使用的代码。

Print[ 
Plot[ 
    Integrate[ 
    Product[ 
    If[R[priors, sigma, d, c] < L[priors, sigma, d, c], 0, 
    (CDF[NormalDistribution[(c - 1) d, sigma], R[priors, sigma, d, c]] - 
     CDF[NormalDistribution[(c - 1) d, sigma], L[priors, sigma, d, c]])] 
    , {c, 1, k}] 
    , {d, 0.01, 5}] 
, {e, -1/k, 1 - 1/k}, PlotRange -> {All, All}, AxesLabel -> Automatic]]; 

Integrating over d

但是,得到的情节是不是我所期望的。这是不变的,在3D图中可以看出,这是不可能发生的。有没有人知道发生了什么,以及如何获得真正的功能集成?提前致谢。

回答

2

当您在函数LR内计算val时,结果为符号(因为e未定义)。逻辑val < minVal因此是不确定的,并因此minVal从未设置(以便LR每次都返回无穷大)

(清理了几个其他的事情以及..)

sigma = 1; 
k = 3; 
priors = Join[ {1/k + e} , Table[1/k - e/(k - 1) , {c, 2, k} ] ]; 
L[priors0_, sigma_, d_, i_, e0_] := Module[{priors, maxVal, val, e}, 
    Do[maxVal = -Infinity; 
    priors = priors0 /. e -> e0 ; 
    Do[val = (2*sigma^2*Log[priors[[i]]/priors[[j]]] + 
     d^2 (j^2 - i^2 + 2 (i - j)))/(2 (j - i) d); 
     If[val > maxVal, maxVal = val];, {j, 1, i - 1}];, {1}]; maxVal]; 
R[priors0_, sigma_, d_, i_, e0_] := Module[{priors, maxVal, val, e}, 
    priors = priors0 /. e -> e0; 
    Min[Table[(2*sigma^2*Log[priors[[j]]/priors[[i]]] + 
     d^2 (i^2 - j^2 + 2 (j - i)))/(2 (i - j) d), {j, i + 1, k}]]]; 
g[d_?NumericQ, c_, e_] := 
    Product[If[R[priors, sigma, d, c, e] < L[priors, sigma, d, c, e], 
0, 
    (CDF[NormalDistribution[(c - 1) d, sigma], R[priors, sigma, d, c, e]] - 
    CDF[NormalDistribution[(c - 1) d, sigma], L[priors, sigma, d, c, e]])], 
    {c, 1, k}]; 
Plot[NIntegrate[g[d, c, e], {d, 0.01, 5}], {e, -1/k, 1 - 1/k}, 
     PlotRange -> {All, All}, AxesLabel -> Automatic] 

enter image description here