Message Boards Message Boards

0
|
4915 Views
|
3 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Can't integrate the subtraction of two Gaussian CDFs

Posted 11 years ago

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.

3 Replies

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 ...

POSTED BY: Udo Krause

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

meaningful integration

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

Attachments:
POSTED BY: Udo Krause

Thank you!

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

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract