Re: NIntegrate bug in Mathematica 6?
- To: mathgroup at smc.vnet.net
- Subject: [mg84222] Re: NIntegrate bug in Mathematica 6?
- From: Norbert Marxer <marxer at mec.li>
- Date: Fri, 14 Dec 2007 04:49:20 -0500 (EST)
- References: <fjoc9o$rit$1@smc.vnet.net>
On 12 Dez., 11:14, vlad <volodymyr.bab... at gmail.com> wrote:
> The following code in Mathematica 6:
>
> Clear[logdist1, pdfLog1, \[Mu]1, \[Sigma]1]
> logdist1 =
> LogNormalDistribution[(\[Mu]1 - \[Sigma]1^2/2), \[Sigma]1];
> pdfLog1[x_] := PDF[logdist1, x]
>
> \[Mu]1 = 0.05;
> \[Sigma]1 = 0.2;
>
> NIntegrate[Min[100, 100 b1]*pdfLog1[b1], {b1, 0, 5}]
> NIntegrate[Min[100, 100 b1]*pdfLog1[b1], {b1, 0, 500}]
> NIntegrate[Min[100, 100 b1]*pdfLog1[b1], {b1, 0, 10000}]
> NIntegrate[Min[100, 100 b1]*pdfLog1[b1], {b1, 0, +\[Infinity]}]
>
> Produces the following output:
>
> 94.1407
>
> 38.1789
>
> 38.1789
>
> 94.1407
>
> Note that the first and the last integrals have upper bounds of 5 and
> \
> [Infinity]
>
> The middle ones have bounds 500 and 10000
>
> All of the answers should be the same (we are way in the tail of the
> random variable density). I get no warnings or errors.
>
> Shouldn't Mathematica send me some warning that it has difficulty with
> convergence? Can I get Mathematica to send me a warning? If not,
> can I trust the numerical integration routines?
>
> Incidentally, Mathematica 5.2 give the correct answer of 94.1407 in
> all four cases.
Hello
With your input ...
Clear[logdist1, pdfLog1, \[Mu]1, \[Sigma]1]
logdist1 = LogNormalDistribution[\[Mu]1 - \[Sigma]1^2/2, \[Sigma]1];
pdfLog1[x_] := PDF[logdist1, x]
\[Mu]1 = 0.05;
\[Sigma]1 = 0.2;
... we can see that Mathematica knows how to integrate the
distribution both numerically and symbolically :
(Integrate[pdfLog1[b], {b, 0, #1}] & ) /@ {5, 500, 10000,
Infinity} // N
(NIntegrate[pdfLog1[b], {b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}
{1., 1., 1., 1.}
{1., 1., 1., 1.}
If we divide the range of integration integration and use the two
ranges [0, 1] and [1, b] then everything is OK.
(Integrate[100*b*pdfLog1[b], {b, 0, 1}] +
Integrate[100*pdfLog1[b], {b, 1, #1}] & ) /@ {5, 500, 10000,
Infinity} // N
(NIntegrate[100*b*pdfLog1[b], {b, 0, 1}] +
NIntegrate[100*pdfLog1[b], {b, 1, #1}] & ) /@ {5, 500, 10000,
Infinity}
{94.1407, 94.1407, 94.1407, 94.1407}
{94.1407, 94.1407, 94.1407, 94.1407}
Therefore the problem seems to be caused by the factor Min[100, 100 b]
in your integral.
Note that we can write this factor in different ways, e.g. (with Min,
Which, Switch, Piecewise, If) as the following plots show :
Plot[(Min[100, 100*#1] & )[b], {b, 0, 5}, PlotRange -> All,
ImageSize -> 200]
Plot[(Switch[100 < 100*#1, True, 100, False, 100*#1] & )[b], {b, 0,
5},
PlotRange -> All, ImageSize -> 200]
Plot[(Which[#1 < 1, 100*#1, True, 100] & )[b], {b, 0, 5},
PlotRange -> All, ImageSize -> 200]
Plot[(Piecewise[{{100*#1, #1 < 1}}, 100] & )[b], {b, 0, 5},
PlotRange -> All, ImageSize -> 200]
Plot[(If[#1 < 1, 100*#1, 100] & )[b], {b, 0, 5}, PlotRange -> All,
ImageSize -> 200]
We note that Integrate works with Min, Which and Piecewise, but NOT
with Switch and If.
N[(Integrate[(Min[100, 100*#1] & )[b]*pdfLog1[b], {b, 0, #1}] & ) /@
{5, 500, 10000, Infinity}]
N[(Integrate[(Switch[100 < 100*#1, True, 100, False, 100*#1] & )[b]*
pdfLog1[b], {b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}]
N[(Integrate[(Which[#1 < 1, 100*#1, True, 100] & )[b]*pdfLog1[b],
{b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}]
N[(Integrate[(Piecewise[{{100*#1, #1 < 1}}, 100] & )[b]*pdfLog1[b],
{b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}]
N[(Integrate[(If[#1 < 1, 100*#1, 100] & )[b]*pdfLog1[b],
{b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}]
{94.1407, 94.1407, 94.1407, 94.1407}
{0., 0., 0., 94.1407}
{94.1407, 94.1407, 94.1407, 94.1407}
{94.1407, 94.1407, 94.1407, 94.1407}
{55.9618, 55.9618, 55.9618, 55.9618}
In contrast NIntegrate does NOT work correctly with all of them and
gives three different answers :
(NIntegrate[(Min[100, 100*#1] & )[b]*pdfLog1[b],
{b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}
(NIntegrate[(Switch[100 < 100*#1, True, 100,
False, 100*#1] & )[b]*pdfLog1[b],
{b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}
(NIntegrate[(Which[#1 < 1, 100*#1, True, 100] & )[
b]*pdfLog1[b], {b, 0, #1}] & ) /@
{5, 500, 10000, Infinity}
(NIntegrate[(Piecewise[{{100*#1, #1 < 1}},
100] & )[b]*pdfLog1[b],
{b, 0, #1}] & ) /@ {5, 500, 10000, Infinity}
(Integrate[(If[#1 < 1, 100*#1, 100] & )[b]*
pdfLog1[b], {b, 0, #1}] & ) /@
{5, 500, 10000, Infinity}
{94.1407, 38.1789, 38.1789, 94.1407}
{0., 0., 0., 94.1407}
{94.1407, 38.1789, 38.1789, 94.1407}
{94.1407, 38.1789, 38.1789, 94.1407}
{55.9618, 55.9618, 55.9618, 55.9618}
I did not expect this.
I' m using Version 6.0.1.0 on Windows (32 - bit).
Best Regards
Norbert Marxer