MathGroup Archive 2008

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

Search the Archive

Re: NIntegrate[UnitStep[...]PDF[...],{x,...}] hard to integrate

  • To: mathgroup at smc.vnet.net
  • Subject: [mg93451] Re: NIntegrate[UnitStep[...]PDF[...],{x,...}] hard to integrate
  • From: Bill Rowe <readnews at sbcglobal.net>
  • Date: Sat, 8 Nov 2008 04:01:00 -0500 (EST)

On 11/7/08 at 6:02 AM, plindsay at mcs.st-and.ac.uk (peter lindsay)
wrote:

>2008/11/6 Bill Rowe <readnews at sbcglobal.net>

>>On 11/5/08 at 4:57 AM, erwann.rogard at gmail.com (er) wrote:

>>>f[x_]=(0.014999775454701741*E^(5.264*x) +
>>>E^(2.632*x)*(-0.012692488700608462 - 0.14964297032258167*x))/
>>>(12.579644604014753 + 7.102251209677398*E^(2.632*x) + E^(5.264*x))

>>>g[t_]:=NIntegrate[ UnitStep[t-f[x]] PDF[
>>>NormalDistribution[0,1/2],x], {x,-Infinity,Infinity} ]

>>>g[0] runs for a very long time before forcing the kernel to quit,
>>>let alone FindRoot[g[t]==0.25,{t,-1,1}]

>>A plot of f[x] shows it everywhere positive except the interval
>>from 0 to ~.8. This can be refined by doing

>>In[8]:= FindRoot[f[x], {x, .8}]

>>Out[8]= {x->0.84719}

>>So, UnitStep[t-f[x]] will be 0 outside this interval when t = 0.
>>This means g[0] amounts to integrating the PDF of a normal
>>distribution from 0 to .84719. For any distribution the integral of
>>the PDF is the CDF. The median of NormalDistribution[0,a] is 0. So,
>>g[0] must be

>>In[10]:= CDF[NormalDistribution[0, 1/2], .84719] - .5

>>Out[10]= 0.454903

>Mathematica 6 gives

>g[0] = 0.438019 for the first problem

I used Mathematica 6 to get the result above. But, I did not
actually solve for the root of the function near 0. Instead, I
took f[0] to be 0 based on the plot. Had I actually used
=46indRoot to get the root near zero then I get

In[20]:= low = FindRoot[f[x], {x, 0}]

Out[20]= {x->0.0211593}

So, using the approach I outlined above the answer to the first
problem would be

In[21]:= high = FindRoot[f[x], {x, 1}]

Out[21]= {x->0.84719}

In[19]:= CDF[NormalDistribution[0, 1/2], x /. high] -
  CDF[NormalDistribution[0, 1/2], x /. low]

Out[19]= 0.438026

which is reasonably close to the result you got.

The point is integrating the product of UnitStep[t-f[x]] and any
probability distribution function must be exactly equivalent to
the difference between the cumulative distribution function at
two points. These two points have to be the roots of t-f[x]. It
should be faster to find the roots then compute the difference
in the cumulative distribution function at the two roots than to
do the integration the original poster set up since the
cumulative distribution function does not need to solve a
general numerical integration problem.

I can make use of this to re-define g as

g[t_?NumericQ] :=
  Block[{x},
   CDF[NormalDistribution[0, 1/2], x /. FindRoot[f[x] - t, {x,
1}]] -
    CDF[NormalDistribution[0, 1/2], x /. FindRoot[f[x] - t, {x, 0}]]]

Once done I can now solve the second problem as

In[3]:= FindRoot[g[t] - .25, {t, 0}] // Timing

Out[3]= {0.094948,{t->-0.00124488}}

which is essentially the same answer as you obtained but clearly
much faster.

>and  t -> -0.00124255 for the second. Takes several minutes though.


  • Prev by Date: Storing and Loading Definitions / Emulating Associative Arrays
  • Next by Date: Re: Constructing a Label
  • Previous by thread: Re: NIntegrate[UnitStep[...]PDF[...],{x,...}] hard to integrate
  • Next by thread: Re: NIntegrate[UnitStep[...]PDF[...],{x,...}] hard to integrate