Prev: NI PCI 6220 ADC board
Next: finding remainder and quotient of transfer function
From: Pr B on 12 Apr 2010 17:27
i'm trying to estimate an integration interval for the indefinite integral of a standard normal distribution such that the result has 10 decimals of precision. what's the best way to do this for this integral (or for any integral in general)? is it basically a trial and error process using a quadrature rule (quad, trap, etc.) and checking the result for 10 decimals of precision given an interval, or is there a better way to do this?
From: Walter Roberson on 12 Apr 2010 18:04
Pr B wrote:
> i'm trying to estimate an integration interval for the indefinite
> integral of a standard normal distribution such that the result has 10
> decimals of precision. what's the best way to do this for this integral
> (or for any integral in general)?

Where there is a relatively efficient method to evaluate the integral
(especially if it has a closed form), evaluate from N to infinity and solve
for the result being the required accuracy, if the integral is one-sided. If
the integral is two-sided (as is the case for the normal distribution), use a
similar technique but solve for the result being one half of the required
accuracy; then the integral from -N to +N will have the full accuracy (half if
it being lost to each of the tails.)

> is it basically a trial and error
> process using a quadrature rule (quad, trap, etc.) and checking the
> result for 10 decimals of precision given an interval, or is there a
> better way to do this?

Only if integrating to infinity is not viable. Due to limited precision and
numeric cancellation and other such issues, expanding your range outwards
until you find a point of 10 decimals precision is relatively likely to give
you an inaccurate answer.

I was working with an integral a day or two ago in which the "expand outwards"
numeric approach gave wildly inaccurate answers -- answers that the numeric
calculations were certain that the solution was around 2.92+4.73i, whereas the
true answer was provably real-valued and on the order of 1E-10. It is much
better to solve analytically for the error when you can, than to work
numerically until your system thinks the answer is within the required tolerance.

Which reminds me, that I was working with a different system about 2 weeks
ago, in which to 5000 or 10000 decimal places, it was "obvious" that the real
part of the answer was converging to about Pi*10^50 and that the imaginary
part of the answer was converging to 0 -- but the closed form of the summation
said that the summation value was infinite, not the value that was obviously
staring me in the face. I worked on it for a couple of hours, trying to figure
out what was going on, until I finally thought to take the limit of the
summand as the parameter approached infinity, and saw that the summand
converged to 10000. Well, if you add up something that is at least as big as
10000 an infinite number of times, the answer is of course going to be
infinite -- but my perception bias that it was converging to Pi (which I
cross-checked using 10000 working digits, remember, so I had Evidence) allowed
me to miss that the answer would be ever-increasing. If I had stopped at the
first point at which the difference in adjacent summands was less than 1 part
in 10^10, I would never have discovered the fact the summation did not
converge. Numeric evidence of convergence is not to be trusted without
analysis proving convergence!
From: John D'Errico on 12 Apr 2010 18:32
"Pr B" <pb2297(a)columbia.edu> wrote in message <hq0378$sci$1(a)fred.mathworks.com>...
> i'm trying to estimate an integration interval for the indefinite integral of a standard normal distribution such that the result has 10 decimals of precision. what's the best way to do this for this integral (or for any integral in general)? is it basically a trial and error process using a quadrature rule (quad, trap, etc.) and checking the result for 10 decimals of precision given an interval, or is there a better way to do this?

The better way to do it is by using the stats toolbox,
which ahs functions to do it for you.

If you don't have that toolbox, then use a transformation
of the functions erf, erfc, erfinv and erfcinv. Don't try to
dop this using quad, trapz, etc. Use tools that are
already carefully written to be efficient and accurate
on this exact problem.

John
 | 
Pages: 1
Prev: NI PCI 6220 ADC board
Next: finding remainder and quotient of transfer function