Search code examples
plotwolfram-mathematicagaussianintegrate

Can't integrate the subtraction of two Gaussian CDFs [Wolfram Mathematica]

I have a function which is a product composed by three (k) factors . Each factor is a subtraction of two Gaussian CDF with random variables R and L. These random variables are defined according to 4 parameters.

enter image description here

The code below shows how I plot the main function (according to two independent variables d and e) and how the random variables are calculated

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

Now, I want to integrate the function over d (in the same region as the Plot3D, d=0.01-5) and to plot the results according to just the independent variable e.

enter image description here

Below is the code I've used. 

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

However, the resulting plot is not what I expect. It's constant and in the 3D plot it can be seen that this cannot happen. Does anyone know what is happening and what to do to obtain the real integration of the function? Thanks in advance.

Solution

  • When you compute val inside the functions L and R the result is symbolic (because e is not defined ). The logical val < minVal is thus indeterminate, and as a result minVal is never set (so that L and R return infinity every time )

    (cleaned up a couple of other things as well.. )

     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