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.