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>