Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Evaluating a double integral with bivariate normal distribution yileds widely different results depending on the method used.

I define a bivariate normal distribution with {10, 3} and {8, 1.5} as mean and standard deviation for each and a correlation coefficient of 0.5.

f[x_ , y_] = 
  Simplify[PDF[MultinormalDistribution[{10,8}, {{9, 0.5 1.5 3}, {0.5 1.5 3, 2.25}}], 
    {x, y}]];
F[x_ , y_] = 
  Simplify[CDF[MultinormalDistribution[{10, 8}, {{9, 0.5 1.5 3}, {0.5 1.5 3, 2.25}}], 
    {x, y}]]
Fx[x_] = CDF[NormalDistribution[10, 3], x];

I wish to evaluate and plot a function ϕ which has the following definition :

ϕm[q_] = 
  (1/(1 - Fx[q]))*Integrate[Integrate[y*(x - q)*f[x, y], {x, q, 25}], {y, 2, 15}]

It can be done it in two ways. Either by using NIntegrate or by using the above double integral.

phim = Interpolation[ParallelTable[{q, ϕm[q]}, {q, 1, 25}]];
Plot[phim[q], {q, 2, 25}]

Using NIntegrate yields a completely different result.

ϕTs = 
  Table[
    {n, 1/(1 - Fx[n]) NIntegrate[(x - n) y f[x, y], {x, n, 25}, {y, 2, 15}]}, 
    {n, 1, 25, .5}];
ϕs = Interpolation[ϕTs, InterpolationOrder -> 1];
g2 = Plot[ϕs[x], {x, 1, 25}, AxesOrigin -> {0, 0}]

Why is that so?

share|improve this question

1 Answer 1

First, you need to use SetDelayed(:=),not Set(=) when we define functions.

Clear["Global`*"]
f[x_, y_] := Evaluate@Simplify[
    PDF[MultinormalDistribution[{10, 8}, {{9, 0.5 1.5 3}, {0.5 1.5 3, 2.25}}], {x, y}]];
F[x_, y_] := Evaluate@Simplify[
    CDF[MultinormalDistribution[{10, 8}, {{9, 0.5 1.5 3}, {0.5 1.5 3, 2.25}}], {x, y}]]
Fx[x_] := Evaluate@CDF[NormalDistribution[10, 3], x];

Now check it:

?f

f[x_,y_]:=2.02192*10^-8 E^(-0.0740741 x^2+x (0.296296 +0.148148 y)+(3.25926 -0.296296 y) y)

Second,it is hard for Integrate to get the answer,so just use NIntegrate.

ans1 = Table[{q, 
  1/(1 - Fx[q]) NIntegrate[(x - q) y f[x, y], {x, q, 25}, {y, 2, 15}]}, {q, 1, 25}]

{{1, 74.3559}, {3, 58.8851}, {5, 44.7569}, {7, 33.152}, {9, 24.6134}, {11, 18.7654}, {13, 14.8523}, {15, 12.2048}, {17, 10.3608}, {19, 9.02093}, {21, 7.92034}, {23, 6.09373}, {25, 0.}}

Last,check the ans,we use the definition of integral:Riemann integral http://en.wikipedia.org/wiki/Riemann_integral

myintegrate[q_] := Module[{dx = (25 - q)/100, a, b, c},
  a = Total@Flatten@Table[(x - q) y f[x, y], {x, q, 25, dx}, {y, 2, 15, (15 - 2)/100}];
  b = (25 - q) (15 - 2)/100^2;
  c = 1/(1 - Fx[q]);
  {q, a b c}
 ]
myintegrate /@ Range[1, 25, 2]

{{1, 74.356}, {3, 58.8849}, {5, 44.7562}, {7, 33.1502}, {9, 24.6101}, {11, 18.7609}, {13, 14.8471}, {15, 12.1997}, {17, 10.3568}, {19, 9.01942}, {21, 7.92463}, {23, 6.1129}, {25, 0.}}

Now we can believe that the results of NIntegrate are correct.

share|improve this answer
    
Thanks! This is exactly what I hoped! Plus, you give me interesting hints on usage of module. –  Xavier_B Jul 3 at 6:03
1  
@Xavier_B If you return to the site, you should upvote and accept this answer, since it worked for you. (Accepting the best, correct answer to your problem the normal way of marking such answers on the site. Plus it's the expected way of thanking those who gave you help.) –  Michael E2 Aug 1 at 17:37

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.