Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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























  • Prev by Date: Re: Joining two Tables or something like that.
  • Next by Date: Re: Have I found a bug?
  • Previous by thread: Re: NIntegrate bug in Mathematica 6?
  • Next by thread: diff eq