# Comparing NIntegrate with definite integration

I want to calculate without error a sum of double integrals with the following definitions for the functions within them:

fx[ x_ ] = PDF[NormalDistribution[10, 2], x];
Fx[ x_ ] = CDF[NormalDistribution[10, 2], x];
f[ x_ , y_ ] = Simplify[PDF[MultinormalDistribution[{10, 8},
{{2^2, 0.5*2*1.5}, {0.5 2 1.5, 2.25}}], {x, y}]];
\[Phi]m[q_ ] = Integrate[Integrate[(-q + x)*y*f[x, y],
{x, q, 25}], {y, 2, 15}]/(1 - Fx[q]);


When I compute the differential in q of [Phi][q] using two different methods, I get very different results. Compare:

D1[q_] = D[\[Phi]m[q], q];


With:

Td\[Phi]s = Table[{q,
1/(1 - Fx[q]) NIntegrate[-y f[x , y ], {x, q, 25}, {y, 2, 15}] +
fx[q]/(1 - Fx[q])^2 NIntegrate[(x - q) y f[x, y], {x, q, 25}, {y, 2, 15}]},
{q, 1, 23, 0.1}];
d\[Phi]s = Interpolation[Td\[Phi]s, InterpolationOrder -> 1]


Look at the difference in results between these 2:

Table[{D1[q], d\[Phi]s[q]}, {q, 4, 16}]

    {{0.0599831, -7.89327}, {0.213901, -7.6447}, {0.584625, -7.09157},
{1.25085, -6.1804}, {2.17537, -5.03748}, {3.22484, -3.88942},
{4.26853, -2.90709}, {5.23671, -2.14794}, {6.11377, -1.59297},
{6.91058, -1.19675}, {7.64478, -0.915004}, {8.33155, -0.71338},
{8.97959, -0.567722}}


If somebody can explain which is the best method (apart from the difference in computing time), I would be very grateful! Thanks. Xavier

-
I think your second expression just isn't right. It cant be that simple, there should be a bunch of D[Fx[q],q] terms and such. – george2079 Jun 24 '14 at 18:31
The expression of d[Phi]s is in fact the differentiation of [Phi]m in q when Mthematica has not yet replaced either f[x,y] nor Fx[q] by their respective definitions. – Xavier_B Jun 25 '14 at 8:24