Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
516 views
in Technique[技术] by (71.8m points)

r - integrate() gives totally wrong number

integrate() gives horribly wrong answer:

integrate(function (x) dnorm(x, -5, 0.07), -Inf, Inf, subdivisions = 10000L)
# 2.127372e-23 with absolute error < 3.8e-23

The return value should obviously be 1 (normal distrubution integrates to 1), but integrate() returns ridiculously small number, with wrong error reporting, and no warning...

Any ideas?

This seems the default integrate() is horribly buggy... and I just found this by chance! Is there any reliable R package to compute numerical integration?

EDIT: I tried package pracma and I see the same problem! :

require(pracma)
integral(function (x) dnorm(x, -5, 0.07), -Inf, Inf)
# For infinite domains Gauss integration is applied!
# [1] 0

EDIT: Hmm... digging deeper, it seems that he has trouble to find the very narrow domain for the function which is numerically > 0. When I set the limits to certain (very close to 0, 1) quantiles, it starts to work:

integral(function (x) dnorm(x, -5, 0.07), qnorm(1e-10, -5, 0.07), qnorm(1 - 1e-10, -5, 0.07))

But anyway, this is quite horrible gotcha... wonder if there is any remedy for this.

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

From the online documentation: "Like all numerical integration routines, these evaluate the function on a finite set of points. 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."

I take this to mean "caveat emptor". I notice that in your example, the absolute error is greater than value of the integral. Given that you know f(x) > 0 for all x, at least it's giving you the chance to spot that something has gone wrong. It's down to you to take the opportunity.

integrate( function(x) dnorm(x, -5, 0.07), -20, 10, subdivisions=1000L) 

Gives

1 with absolute error < 9.8e-07

The warning in the online doc says to me that, given your apparent definition of buggy, the answer to your question is "no, there is no reliable numerical intergration method. Not in R or any other language". No numerical integration technique should be used blindly. The user needs to check their inputs are sensible and the output is reasonable. It's no good believing an answer just because the computer gave it to you.

See also this post.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
...