# Integrate and NIntegrate yield different results for double integral

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?

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.

• Thanks! This is exactly what I hoped! Plus, you give me interesting hints on usage of module. – Xavier_B Jul 3 '14 at 6:03
• @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 '14 at 17:37