Re: strange behavior of quad ?

In article <e2j8kh$4rc$1@xxxxxxxxxxxxxxxxxxxxxxxxx>, "Gg"
<nc-schuetes@xxxxxxxxxxxxx> wrote:

Dear All,

I am trying to do some analysis involving lognormal (and other)
distributions which often involves indefinite integrals (i.e. integrals with
limits -inf and/or +inf).
I was trying to do this like

quadl(@(x)(lognpdf(x,1,2)),lowlim, hilim)

where lowlim and hilim are very small or very large numbers supposedly
proxies for -inf and +inf.
I noticed the following oddity of quad and quadl:

First Test:
----------------------------------------------
quad(@(x)(lognpdf(x,1,2)),0, 1e40)

ans =

0
----------------------------------------------
This sure looks strange, doesn't it ? But it did not give us any warning or
error.

Second test:
----------------------------------------------
quadl(@(x)(lognpdf(x,1,2)),0.0000001, 1e13)

ans =

1.0000

quadl(@(x)(lognpdf(x,1,2)),0.00000000001, 1e13)

ans =

1.2222e-016

quad(@(x)(lognpdf(x,1,2)),0.00000000001, 1e13)

ans =

9.8571e-017

quadl(@(x)(lognpdf(x,1,2)),0, 1e13)

ans =

1.7873e-039
----------------------------------------------

Question:
- Is it fair to say that the implementation of quad and quadl is somewhat
crummy, since it gives us wrong answers without warning ?
- Is there a better way to compute those indefinite integrals ? Especially
considering the fact that I know beforehand that all integrals will
converge.

Note that it might be possible to come up with more or less ingenious ways
to "preprocess" the integrals, so that quad or quadl might work. But
generally my integrands will be more complicated than logn. In addition I
would like to program something generic, which does not require intelligent
intervention.

Thank you for any hints and help
gg
-------------------------
(This is an addition to John's remarks. It took me a while to prepare
so I can't resist sending it.)

I would say you are putting an unusual strain on 'quadl' for it to be
able to encompass a range of x all the way from 0 to 1e40 with a function
such as lognpdf(x,1,2). Remember, 'quadl' is a numerical integration
algorithm with all the shortcomings that numerical procedures have. When
presented with such a huge range for which all but the initial portion -
minuscule relative to the whole range - is essentially zero as far as
numerical computation is concerned, the algorithm does some searching over
possible discrete sets of x values but in this case it never sees the part
that hasn't been rounded down to zero. You are asking an awful lot in
giving such an extreme range. Even in the shorter ranges up to 1e13, this
is still giving rise to a large amount of inaccuracy because so little
sampling of the important nonzero range is done.

To achieve reasonable accuracy in such extreme cases, you should be
breaking up your range into two or more segments in which the first
segment encompasses the part where the bulk of the integrands' appreciably
nonzero values reside, with other added segments intended only for adding
on a small degree of accuracy. I think 'quadl' is not smart enough to do
all that for you in one call.

In short, to do the kind of problem you envision will require a bit of
careful planning to get the proper answers. It cannot be completely
automatic. The individual characteristics of the distributions you are
dealing with need to be taken into account.

I would not characterize the existing quadrature functions as "crummy".
Perhaps with a lot of hard word and planning, more versatile quadrature
procedures could be designed to work properly in the face of such extreme
skewness of distribution, but very likely it would be at an appreciable
cost in execution time for the more ordinary kinds of integration. They
would have to resort to some very extensive searching patterns.

By the way, it is misleading to refer to integrals with infinities in
their range of integration, as being "indefinite" integrals. That is a
term reserved in calculus for integrals which are the reverse derivatives
of functions and have an arbitrary additive constant. "Improper" is the
accepted term for integrals with infinities in their range of integration
(or in their integrand values.)

(Remove "xyzzy" and ".invalid" to send me email.)
Roger Stafford
.

Relevant Pages

  • strange behavior of quad ?
    ... distributions which often involves indefinite integrals (i.e. integrals with ... Is it fair to say that the implementation of quad and quadl is somewhat ... - Is there a better way to compute those indefinite integrals? ...
    (comp.soft-sys.matlab)
  • Re: strange behavior of quad ?
    ... distributions which often involves indefinite integrals (i.e. integrals ... QUAD is an adaptive quadrature routine. ... integration into two subintervals, ... QUADL works using a similar algorithm to the one QUAD uses, ...
    (comp.soft-sys.matlab)
  • Re: quadl integration plot
    ... like to use quadl. ... Here's an example I just threw together using QUADGK. ... % Y is the array of from integrals from Xto X. ...
    (comp.soft-sys.matlab)
  • Re: Wanted: Lazarus Long rebuttal quote
    ... In rec.arts.sf.written Mike Schilling wrote: ... was the fact that in integrating there were perfectly valid ... Did the exams really give you difficult indefinite integrals? ...
    (rec.arts.sf.written)
  • Re: strange behavior of quad ?
    ... distributions which often involves indefinite integrals (i.e. integrals with ... - Is there a better way to compute those indefinite integrals? ... so that quad or quadl might work. ... The best material model of a cat is another, or preferably the same, cat. ...
    (comp.soft-sys.matlab)
Loading