Re: NIntegrate[UnitStep[...]PDF[...],{x,...}] hard to integrate
- To: mathgroup at smc.vnet.net
- Subject: [mg93464] Re: NIntegrate[UnitStep[...]PDF[...],{x,...}] hard to integrate
- From: er <erwann.rogard at gmail.com>
- Date: Sun, 9 Nov 2008 05:26:34 -0500 (EST)
- References: <gf3mj0$en2$1@smc.vnet.net>
On Nov 8, 4:36 am, Bill Rowe <readn... at sbcglobal.net> wrote:
> On 11/7/08 at 6:02 AM, plind... at mcs.st-and.ac.uk (peter lindsay)
> wrote:
>
>
>
> >2008/11/6 Bill Rowe <readn... at sbcglobal.net>
> >>On 11/5/08 at 4:57 AM, erwann.rog... 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.
Much helpful, thanks!
OK, but you have to assume that there are exactly 2 roots to f[x]-
t,right?
That may be the case for particular values of t, but for others I may
have, for example, a left interval (only 1 root) etc. So, what if I
need say g[t]==0.95, instead of 0.25? and what if I need to carry out
these computations for various definitions of f?
- Follow-Ups:
- Re: Re: Re: How do little quickest ?
- From: DrMajorBob <btreat1@austin.rr.com>
- Re: How do little quickest ?
- From: DrMajorBob <btreat1@austin.rr.com>
- How do little quickest ?
- From: Artur <grafix@csl.pl>
- Re: Re: Re: How do little quickest ?