From maechler at stat.math.ethz.ch  Tue Dec  7 09:28:16 2010
From: maechler at stat.math.ethz.ch (Martin Maechler)
Date: Tue, 7 Dec 2010 09:28:16 +0100
Subject: [Rd] 0.5 != integrate(dnorm,0,20000) = 0
In-Reply-To: <alpine.LFD.2.00.1012070734070.10594@gannet.stats.ox.ac.uk>
References: <4CFD729C.6050704@structuremonitoring.com>
	<alpine.LFD.2.00.1012070734070.10594@gannet.stats.ox.ac.uk>
Message-ID: <19709.61472.158760.44576@lynne.math.ethz.ch>

>>>>> Prof Brian Ripley <ripley at stats.ox.ac.uk>
>>>>>     on Tue, 7 Dec 2010 07:41:16 +0000 (GMT) writes:

    > On Mon, 6 Dec 2010, Spencer Graves wrote:
    >> Hello:
    >> 
    >> 
    >> The example "integrate(dnorm,0,20000)" says it "fails on many systems". 
    >> I just got 0 from it, when I should have gotten either an error or something 
    >> close to 0.5.  I got this with R 2.12.0 under both Windows Vista_x64 and 
    >> Linux (Fedora 13);  see the results from Windows below.  I thought you might 
    >> want to know.

    > Well, isn't that exactly what the help page says happens?  That 
    > example is part of a section entitled

    > ## integrate can fail if misused

    > and is part of the illustration of

    > If the function is
    > approximately constant (in particular, zero) over nearly all its
    > range it is possible that the result and error estimate may be
    > seriously wrong.

yes, of course, 
and the issue has been known for ``ages''  ..
..........
..........
but it seems that too many useRs are not reading the help
page carefully, but only browse it quickly.
I think we (R developers) have to live with this fact
and should consider adapting to it a bit more, particularly in
this case (see below)

    >> 
    >> Thanks for all your work in creating and maintaining R.
    >> 
    >> 
    >> Best Wishes,
    >> Spencer Graves
    >> ###############################

    >> 
    >> integrate(dnorm,0,20000) ## fails on many systems

    >> 0 with absolute error < 0

and this is particularly unsatisfactory for another reason:

"absolute error < 0"   
is *always* incorrect, so I think we should change *some*thing.

We could just use "<=" (and probably should in any case, or  
"< ~= x" which would convey ``is less than about x'' which I
think is all we can say),
but could consider giving a different message when the integral
evaluates to 0 or, rather actually,
only when the error bound ('abs.error') is 0 *and* 'subdivisions == 1'
as the latter indicates that the algorithm treated the integrand
f(.) as if f() was a linear function.

But in my quick experiments, even for linear (incl. constant)
functions, the 'abs.error' returned is never 0.

If we want to be cautious,
such a warning could be made explicitly suppressable by an argument
      .warn.if.doubtful = TRUE

An additional possibility I'd like to try, is a new argument
   'min.subdivisions = 3' which specifies the *minimal* number
of subdivisions to be used in addition to the already present
   'subdivisions = 100' (= the maximum number of subintervals.)

Martin Maechler,
ETH Zurich



[Rd] 0.5 != integrate(dnorm,0,20000) = 0
Martin Maechler