Matlab - Numerical Integration Question

Author

Numerical Integration Question

junoexpress

2006-12-27, 7:09 pm

Hi,

I am having a strange problem integrating a function.
I am intgerating the following function:

I(a,b) = Integral from a to b of F(x,Sigma)*sqrt(1+x)
where
F(x,Sigma) = 1/(Sigma*sqrt(2*pi)) * exp(-x^2/(2*Sigma^2)) PDF
for RV ~ N(0,Sigma)
where a is some fixed number such that -1 < a < 0
b is a number which is essentially infinity.

I find that when I make b very large (b>= 4), I(a,b) is approximately
zero for some of the a values (which it should not be), yet when b has
a smaller value (1<b<3), I(a,b) is about 0.99 for all of the a vvalues
(which it should be).

I am using the quadl command in MatLab, although I notice the same
phenomena with Mathematica: when I use 1<b<3 OR b = infinity, I get the
right answer, but with b > 1000, I also get an answer of zero along
with the warning message:

NIntegrate::ploss: Numerical integration stopping due to loss of
precision. \
Achieved neither the requested PrecisionGoal nor AccuracyGoal; suspect
one of \
the following: highly oscillatory integrand or the true value of the
integral \
is 0. If your integrand is oscillatory try using the option \
Method->Oscillatory in NIntegrate.

I varied the AccuracyGoal and while I lose the warning sign, I still
get a value of 0.

Thanks for any help you can provide as to explaining why this result
occurs.

Matt Brenneman

Roger Stafford

2006-12-27, 7:09 pm

In article <1167255020.411893.113530@f1g2000cwa.googlegroups.com>,
"junoexpress" <MTBrenneman@gmail.com> wrote:

>SNIP - Having trouble with numerical integration
> ......
> Matt Brenneman
--------------------------------
You are probably using 'quadl' with a very small value of sigma. If so,
you are dealing with an integrand that is almost identically zero for all
but values of x very near zero. Numerical integration routines have
trouble with that kind of situation. In their haste to provide rapid
answers they start out initially by sampling the integrand at various, but
not too numerous, points throughout its range. If they encounter what
looks like all zeros, they often give up quickly without doing further
testing, under the assumption that the integrand actually is identically
zero. You need to bring the limits of integration in sufficiently close
to zero on either side so that 'quadl' won't be fooled in this manner and
yet far enough away to provide an accurate answer.

A good way to test for this situation is to make a preliminary plot of
the integrand over the full planned integration range and see if the plot
shows a sharp spike at one point and zeros everywhere else. If so, the
range should be brought in close enough to avoid this.

Roger Stafford

junoexpress

2006-12-27, 10:07 pm

Roger Stafford wrote:
> In article <1167255020.411893.113530@f1g2000cwa.googlegroups.com>,
> "junoexpress" <MTBrenneman@gmail.com> wrote:
>
> --------------------------------
> You are probably using 'quadl' with a very small value of sigma. If so,
> you are dealing with an integrand that is almost identically zero for all
> but values of x very near zero. Numerical integration routines have
> trouble with that kind of situation. In their haste to provide rapid
> answers they start out initially by sampling the integrand at various, but
> not too numerous, points throughout its range. If they encounter what
> looks like all zeros, they often give up quickly without doing further
> testing, under the assumption that the integrand actually is identically
> zero. You need to bring the limits of integration in sufficiently close
> to zero on either side so that 'quadl' won't be fooled in this manner and
> yet far enough away to provide an accurate answer.
>
You were right on the money. Just to get an idea of where the
integration bounds "should" end, I went back and solved for the x that
approximately makes the integrand = 0:

exp(-x^2/(2*Sigma^2)) = realmin = 2.2x10^(-308)

I found that (with the Sigma I was working with) x was approx 1.8,
which puts me about smack of the range I quoted earlier (1-3 for the
upper integration range). Due to the sqrt(1+x), which increases the
integrand, a slightly larger value for the upper integration bounds
works also.

Thank you very much,

Matt Brenneman

Sponsored Links