无法积分两个高斯 CDF 的减法 [Wolfram Mathematica]答案

【问题标题】：Can't integrate the subtraction of two Gaussian CDFs [Wolfram Mathematica]无法积分两个高斯 CDF 的减法 [Wolfram Mathematica]
【发布时间】：2014-10-27 15:51:38
【问题描述】：

我有一个函数，它是由三个 (k) 因素组成的乘积。每个因子是两个带有随机变量 R 和 L 的高斯 CDF 的减法。这些随机变量是根据 4 个参数定义的。

下面的代码显示了我如何绘制主函数（根据两个自变量 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]];

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

以下是我使用的代码。

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]];

但是，结果不是我所期望的。它是恒定的，在 3D 图中可以看出这不可能发生。有谁知道发生了什么以及如何获得功能的真正集成？提前致谢。

【问题讨论】：

标签： plot wolfram-mathematica gaussian integrate

【解决方案1】：

当您在函数 L 和 R 内计算 val 时，结果是符号的（因为 e 未定义）。因此，逻辑val < minVal 是不确定的，因此永远不会设置minVal（因此L 和R 每次都返回无穷大）

（还清理了一些其他的东西......）

 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]

【讨论】：