Can't integrate the subtraction of two Gaussian CDFs

Questions
,
1

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.

latex-image-1.jpg
latex-image-1.jpg1177×220 35.2 KB

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.

2

This is extraordinarily dull, the result is not the integral, but the length of the integration interval, if one integrates {d, 0.01, 4.2} - guess what - the result is 4.2 - 0.0.1 = 4.19 and so on

In[116]:= NIntegrate[
 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, 4.2}, 
 PrecisionGoal -> 30, WorkingPrecision -> 60]

Out[116]= 4.19000000000000017742751712290782961645163595676422119140625

In[117]:= NIntegrate[
 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, 2.7}, 
 PrecisionGoal -> 30, WorkingPrecision -> 60]

Out[117]= 2.69000000000000017742751712290782961645163595676422119140625

because the max value of the CDF is 1 that means that somehow under integration the functions do not evaluate or have always their maximum value and so the result of integration is just 1 * (dmax - dmin).

Your usage of the Do is at least strange. It could be rewarding to make variables explicit - the dependency from e is working only through the priors table …

3

Now it works, specify d in the integration like d_?NumericQ forcing it into numeralization

comm377583.PNG
comm377583.PNG747×526 9.37 KB

see the notebook for a more conventional treatment of R and L and a more explicit treatment of variables, too. Gosh!

4

Thank you!

I was getting crazy with this error! Now it works perfectly!