# Double integral over normal pdf gives inconsistent answer

I am using Integrate function to do the following double integral

Integrate[3.1 x^2.1 Exp[-y^2/2]/(2 π)^0.5, {y, -1, 3}, {x, 0, y + 35}]


The result is 1204.63.

When I check with the following integral

N[Integrate[Integrate[3.1 x^2.1, {x, 0, y + 35}, Assumptions :> y ∈ Reals]*
Exp[-y^2/2]/(2 π)^0.5, {y, -1, 3}]]


The result is 52774.4, which also agrees with what I got in Matlab. I thought the first one is a nice and clean way to do double integral than the second one, but it seems that the result is not what I want if the function w.r.t. x (x^2.1 in this case) does not have integer power.

I would be very thankful if anyone can answer my question.

• NIntegrate[3.1 x^2.1 Exp[-y^2/2]/(2 \[Pi])^0.5, {y, -1, 3}, {x, 0, y + 35}] is perhaps a nice compromise. Integrate is primarily for exact computation. It tries with inexact input, but in this case fails. Note that, with exact input, N@Integrate[31/10 x^(21/10) Exp[-y^2/2]/(2 \[Pi])^(1/2), {y, -1, 3}, {x, 0, y + 35}] also gives the right answer. Feb 14, 2015 at 18:22
• BTW, you don't seem to be asking a question. No "?" for example. Just observations about the behavior of Integrate. It's not clear whether you want an explanation or other workarounds or what. Feb 14, 2015 at 18:24
• Switch Integrate to NIntegrate. It seems that the machine precision values in your expression make it impossible to get the correct answer symbolically when posed in this way. Feb 14, 2015 at 18:29
• Thanks Michael And Oleksandr very much for your comment, it resolves my concern definitely! Sorry for being unable to make the question clear. I was using the first expression a lot to get numerical experiments without noticing the existence problem of inconsistency like this. I thought the second one is correct because it got the same result with Matlab, I was just stuck in the confusion and unable to figure out the problem with the first expression. BTW, when the power of x is integer instead of 2.1 in this case, the results of both expressions always agree. Thanks again :) Feb 14, 2015 at 18:42
• Related (perhaps a duplicate?): mathematica.stackexchange.com/questions/51809/… Feb 14, 2015 at 19:27

This does not answer the question why but I post for interest. The region of integration:

ir = ImplicitRegion[0 < x < y + 35 && -1 < y < 3, {x, y}];
RegionPlot[ir] The integral must be < Area[ir]f[38,0]:369981. (not a helpful bound). $0<x<38 \land -1<y<3$ would be a closer.

SetAttributes[dis, HoldFirst];
dis[u_] := {Style[#, Bold], ReleaseHold[#]} & @ HoldForm @ u;
f[x_, y_] := 3.1 x^2.1 Exp[-y^2/2]/(2 π)^0.5;
r1 = dis[NIntegrate[f[x, y], {y, -1, 3}, {x, 0, 35 + y}]];
r2 = dis[NIntegrate[f[x, y], {x, 0, 38}, {y, -1, 3}]];
r3 = dis[NIntegrate[f[x, y], {x, y} ∈ ir]];
r4 = dis[Integrate[f[x, y], {x, y} ∈ ir]]
Grid[{r1, r2, r3,r4}, Alignment -> Left, Frame -> All]


As has been observed by MichaelE2 and OleksandrR numerical integration yields correct result. Interestingly, Integrate works with region (presumably switching to numerical integration). The results are summarized below: • Yes, the last two entries are effectively equivalent. Integrate over a region will call NIntegrate with WorkingPrecision set to the minimum precision of the integrand and region, when the precision is less than infinity. Feb 15, 2015 at 13:38