473,581 Members | 2,224 Online
Bytes | Software Development & Data Engineering Community
+ Post

Home Posts Topics Members FAQ

When (32-bit) double precision isn't precise enough

I've been trying to perform a calculation that has been running into
an underflow (insufficient precision) problem in Microsoft Excel, which
calculates using at most 15 significant digits. For this purpose, that
isn't enough.

I was reading a book about some of the financial scandals of the
1990s called "Inventing Money: The Story of Long-Term Capital Management
and the Legends Behind it" by Nicholas Dunbar. On page 95, he mentions
that the 1987 stock-market crash was designated in economists' computer
models as a "20-sigma event". That is, their models (which obviously were
badly flawed!) put such an event in the exponential "tails" of a normal
Gaussian distribution, outside the range +/- 20 standard deviations from
the mean. In other words - vastly unlikely.

So I wanted to know: just how unlikely is that? Where did they get
the data to support that idea, anyway?

I put together an Excel spreadsheet, attempting to calculate this
probability. But after about 8 sigma I ran into the limits of Excel's
15-digit
precision and could go no further. Since the odds of an 8-sigma event
are about 819 trillion to 1 against, a 20-sigma event is obviously going to
be so unlikely that it's not even worth talking about; would probably
never happen in the lifetime of the universe. Further evidence that those
computer models had some serious problems.

To perform that calculation, I used Excel's built-in error function
ERF(), and here's the output:

s Confidence Interval Probability Odds against

1.0 68.268949131685 700% 31.731050868314 3000000% 3
1.5 86.638554269553 900% 13.361445730446 1000000% 7
2.0 95.449972950730 900% 4.5500270492691 1000000% 22
2.5 98.758066814134 400% 1.2419331858655 6000000% 81
3.0 99.730020387427 800% 0.2699796125722 0200000% 370
3.5 99.953474183672 100% 0.0465258163279 4690000% 2,149
4.0 99.993665751560 700% 0.0063342484393 2140000% 15,787
4.5 99.999320465373 800% 0.0006795346262 3560100% 147,160
5.0 99.999942669685 300% 0.0000573303147 3997840% 1,744,278
5.5 99.999996202087 500% 0.0000037979124 9560668% 26,330,254
6.0 99.999999802682 500% 0.0000001973175 2898267% 506,797,346
6.5 99.999999991968 000% 0.0000000080319 9728949% 12,450,203,405
7.0 99.999999999744 000% 0.0000000002559 6191833% 390,683,116,666
7.5 99.999999999993 600% 0.0000000000063 8378239% 15,664,694,356, 071
8.0 99.999999999999 900% 0.0000000000000 0000000% 818,836,295,885 ,545
8.5 100.00000000000 000% 0.0000000000000 0000000% #DIV/0!
.. . . .
s =ERF(An/SQRT(2)) =1-Bn =1/Cn

(The 'n' in the cell formulas above represents the row number).

Evidently this calculation hits the limit of Excel's computational
precision at about 8 sigma.

For what it's worth, ERF(z) is defined as 2/pi * INT[(0,z) exp(-t²)
dt], and the area under the "normal" Bell curve from -z to z is just
ERF(z/sqrt(2)). There's a discussion of it here:
http://mathworld.wolfram.com/Erf.html, here:
http://mathworld.wolfram.com/ConfidenceInterval.html, and here:
http://jove.prohosting.com/~skripty/page_295.htm.

So I decided to try to tackle this in C. I downloaded a package
called the GNU Scientific Library (http://www.gnu.org/software/gsl/) -
there's a nice precompiled binary for Microsoft Visual C++ at
http://www.network-theory.co.uk/gsl/freedownloads.html. I wrote a little
piece of code to try it out:

=============== =============== =============== =============== ==============
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_erf.h>

int main(void)
{
double odds;
double index;
double sigma;
double result;

for (index = 1.0; index <= 8.5; index += 0.5)
{
sigma = index / M_SQRT2;
result = 1 - (gsl_sf_erf (sigma));
odds = 1 / result;
printf("P(%2.1f \345) = %.18LE \t %-1.18lE\n", index, result, odds);
}
return 0;
}
=============== =============== =============== =============== ==============

And here's its output:

P(1.0 s) = 3.1731050786291 42600E-001 3.1514871875343 75500E+000
P(1.5 s) = 1.3361440253771 61700E-001 7.4842231152268 46800E+000
P(2.0 s) = 4.5500263896358 41700E-002 2.1977894507992 82900E+001
P(2.5 s) = 1.2419330651552 31800E-002 8.0519637334481 30300E+001
P(3.0 s) = 2.6997960632602 06900E-003 3.7039834734495 63900E+002
P(3.5 s) = 4.6525815807108 01700E-004 2.1493443643114 46000E+003
P(4.0 s) = 6.3342483666239 95700E-005 1.5787192767323 96800E+004
P(4.5 s) = 6.7953462494774 18600E-006 1.4715953584806 70300E+005
P(5.0 s) = 5.7330314373604 80700E-007 1.7442778936869 13900E+006
P(5.5 s) = 3.7979124956066 81200E-008 2.6330253821191 82500E+007
P(6.0 s) = 1.9731752898266 56400E-009 5.0679734596101 19500E+008
P(6.5 s) = 8.0319972894926 65000E-011 1.2450203404677 24800E+010
P(7.0 s) = 2.5596191832732 98400E-012 3.9068311666627 59400E+011
P(7.5 s) = 6.3837823915946 50100E-014 1.5664694356071 29100E+013
P(8.0 s) = 1.2212453270876 72200E-015 8.1883629588554 47500E+014
P(8.5 s) = 0.0000000000000 00000E+000 1.#INF000000000 00000E+000
Same limits as in Excel, it seems. GSL's gsl_sf_erf function takes a
double-precision argument and returns double. My Visual C++ compiler (v.
6.0) specifies that both double and long double use an 8-byte
representation: "The long double contains 80 bits: 1 for sign, 15 for
exponent, and 64 for mantissa. Its range is +/- 1.2E4932 with at least 19
digits of precision."

Strangely enough, I only seem to be getting 17 digits of precision.
Regardless, I estimate that this calculation will require at least 38
additional digits of precision. Interesting: according to IEEE 754, the
condition for "positive underflow" (single precision) shouldn't happen
for positive numbers greater than about 1.4E-045 or so. For double
precision, it should happen only for positive numbers less than about
1.0E-308.

So here's what I'd like to know.

Are there 64-bit implementations of something like GSL which would
produce more precise output on appropriate OS/hardware platforms like
WinXP-64, Solaris v. 7-9, Tru64 Unix, Linux, etc? Has anyone
implemented a 128-bit long long double datatype or equivalent?

Can a calculation like this be performed with some kind of arbitrary-
precision numeric package (something like Michael Ring's MAPM or PHP
BCMath), or evaluated directly by Mathematica or the java.math package?

Or maybe I'm just confused and there's an easier way to do this which
I'm just not seeing.

By the way, I don't know of anyone who's bothered to tabulate the
values of this function nearly this far: most such texts (Zwillinger's CRC
Handbooks, Abramowitz & Stegun, Gradshteyn & Ryzhik, etc) only go to about
4.0 sigma.

Any ideas?

-- Dave Schulman (ca******@gate. net)
Nov 13 '05 #1
5 6079
DAVID SCHULMAN wrote:
....
=============== =============== =============== =============== ============== [code] =============== =============== =============== =============== ==============

And here's its output:

P(8.5 s) = 0.0000000000000 00000E+000 1.#INF000000000 00000E+000 .... Can a calculation like this be performed with some kind of arbitrary-
precision numeric package (something like Michael Ring's MAPM or PHP
BCMath), or evaluated directly by Mathematica or the java.math package?
Mathematica, MathCAD, MathLab, Maple all can do such stuff.
Here is what Maple spits out:
Digits := 190: 1 - erf(20.); 1 / %;

-175
.53958656116079 0 10
.18532707668788 893231785722434 367993256559201 748313827910344 5269\
280820793370976 532964246089265 886796811713147 8014046442587\
589709255826797 081322433635062 615830823633134 4839925105294\
176
06688155039 10

Jirka

Nov 13 '05 #2
"DAVID SCHULMAN" <da************ *@verizon.net> writes:

[Long explanation snipped. The poster wants to calculate the exclusion
probability corresponding to 20 standard deviations of a normal
distribution.]

In C99, you can use the `erfc' function. Many C89 implementations provide
it as an extension.

#include <stdio.h>
#include <math.h>

int main (void)
{
printf ("%e\n", erfc (20.0 / sqrt (2.0)));
return 0;
}

This program prints `5.507248e-89'.

Note that although the /mathematical/ definition of `erfc(x)' is
`1.0 - erf(x)', this is not a good way to actually calculate the function.
The `erf' function would have to return a value which is 5.507248e-89 less
than 1.0, which would require at least 294 significant bits to represent.
Usually, floating point numbers are much less precise in C, so `erf' just
returns the closest value it can represent, 1.0.
So I decided to try to tackle this in C. I downloaded a package
called the GNU Scientific Library (http://www.gnu.org/software/gsl/) -
there's a nice precompiled binary for Microsoft Visual C++ at
http://www.network-theory.co.uk/gsl/freedownloads.html. I wrote a little
piece of code to try it out:

=============== =============== =============== =============== ==============
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_erf.h>

int main(void)
{
double odds;
double index;
double sigma;
double result;

for (index = 1.0; index <= 8.5; index += 0.5)
{
sigma = index / M_SQRT2;
result = 1 - (gsl_sf_erf (sigma));
odds = 1 / result;
printf("P(%2.1f \345) = %.18LE \t %-1.18lE\n", index, result, odds);
}
return 0;
}
=============== =============== =============== =============== ==============
While I'm not familiar with the GNU Scientific Library, I would guess that
it also provides an `erfc' function if it provides an `erf' function.
My Visual C++ compiler (v. 6.0) specifies that both double and long
double use an 8-byte representation: "The long double contains 80 bits:
1 for sign, 15 for exponent, and 64 for mantissa. Its range is +/-
1.2E4932 with at least 19 digits of precision."


As explained above `erf(20.0/sqrt(2.0))' can only be represented as a
value different from 1.0 with at least 294 bits for the mantissa.

Martin
Nov 13 '05 #3
Use lcc-win32 C compiler.
Here is a small program that calculates using 350 bits precision. Math
functions are suffixed with the 'q' (qfloat) suffix. lcc-win32 supports standard
350 bits extended floats.

#include <math.h>
#include <stdio.h>
#include <qfloat.h>
int main(void)
{
printf("%.100qe \n",erfcq(20/sqrtq(2.0q)));
return 0;
}

Output
---------
5.5072482372124 673901512455617 149306656149954 695186611353
987433091698372 417656070673834 468842358561507 e-89

Lcc-win32 can be downloaded at no charge from
http://www.cs.virginia.edu/~lcc-win32
"DAVID SCHULMAN" <da************ *@verizon.net> wrote in message
news:t4******** *********@nwrdd c03.gnilink.net ...
I've been trying to perform a calculation that has been running into
an underflow (insufficient precision) problem in Microsoft Excel, which
calculates using at most 15 significant digits. For this purpose, that
isn't enough.

I was reading a book about some of the financial scandals of the
1990s called "Inventing Money: The Story of Long-Term Capital Management
and the Legends Behind it" by Nicholas Dunbar. On page 95, he mentions
that the 1987 stock-market crash was designated in economists' computer
models as a "20-sigma event". That is, their models (which obviously were
badly flawed!) put such an event in the exponential "tails" of a normal
Gaussian distribution, outside the range +/- 20 standard deviations from
the mean. In other words - vastly unlikely.

So I wanted to know: just how unlikely is that? Where did they get
the data to support that idea, anyway?

I put together an Excel spreadsheet, attempting to calculate this
probability. But after about 8 sigma I ran into the limits of Excel's
15-digit
precision and could go no further. Since the odds of an 8-sigma event
are about 819 trillion to 1 against, a 20-sigma event is obviously going to
be so unlikely that it's not even worth talking about; would probably
never happen in the lifetime of the universe. Further evidence that those
computer models had some serious problems.

To perform that calculation, I used Excel's built-in error function
ERF(), and here's the output:

s Confidence Interval Probability Odds against

1.0 68.268949131685 700% 31.731050868314 3000000% 3
1.5 86.638554269553 900% 13.361445730446 1000000% 7
2.0 95.449972950730 900% 4.5500270492691 1000000% 22
2.5 98.758066814134 400% 1.2419331858655 6000000% 81
3.0 99.730020387427 800% 0.2699796125722 0200000% 370
3.5 99.953474183672 100% 0.0465258163279 4690000% 2,149
4.0 99.993665751560 700% 0.0063342484393 2140000% 15,787
4.5 99.999320465373 800% 0.0006795346262 3560100% 147,160
5.0 99.999942669685 300% 0.0000573303147 3997840% 1,744,278
5.5 99.999996202087 500% 0.0000037979124 9560668% 26,330,254
6.0 99.999999802682 500% 0.0000001973175 2898267% 506,797,346
6.5 99.999999991968 000% 0.0000000080319 9728949% 12,450,203,405
7.0 99.999999999744 000% 0.0000000002559 6191833% 390,683,116,666
7.5 99.999999999993 600% 0.0000000000063 8378239% 15,664,694,356, 071
8.0 99.999999999999 900% 0.0000000000000 0000000% 818,836,295,885 ,545
8.5 100.00000000000 000% 0.0000000000000 0000000% #DIV/0!
. . . .
s =ERF(An/SQRT(2)) =1-Bn =1/Cn

(The 'n' in the cell formulas above represents the row number).

Evidently this calculation hits the limit of Excel's computational
precision at about 8 sigma.

For what it's worth, ERF(z) is defined as 2/pi * INT[(0,z) exp(-t²)
dt], and the area under the "normal" Bell curve from -z to z is just
ERF(z/sqrt(2)). There's a discussion of it here:
http://mathworld.wolfram.com/Erf.html, here:
http://mathworld.wolfram.com/ConfidenceInterval.html, and here:
http://jove.prohosting.com/~skripty/page_295.htm.

So I decided to try to tackle this in C. I downloaded a package
called the GNU Scientific Library (http://www.gnu.org/software/gsl/) -
there's a nice precompiled binary for Microsoft Visual C++ at
http://www.network-theory.co.uk/gsl/freedownloads.html. I wrote a little
piece of code to try it out:

=============== =============== =============== =============== ==============
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_erf.h>

int main(void)
{
double odds;
double index;
double sigma;
double result;

for (index = 1.0; index <= 8.5; index += 0.5)
{
sigma = index / M_SQRT2;
result = 1 - (gsl_sf_erf (sigma));
odds = 1 / result;
printf("P(%2.1f \345) = %.18LE \t %-1.18lE\n", index, result, odds);
}
return 0;
}
=============== =============== =============== =============== ==============

And here's its output:

P(1.0 s) = 3.1731050786291 42600E-001 3.1514871875343 75500E+000
P(1.5 s) = 1.3361440253771 61700E-001 7.4842231152268 46800E+000
P(2.0 s) = 4.5500263896358 41700E-002 2.1977894507992 82900E+001
P(2.5 s) = 1.2419330651552 31800E-002 8.0519637334481 30300E+001
P(3.0 s) = 2.6997960632602 06900E-003 3.7039834734495 63900E+002
P(3.5 s) = 4.6525815807108 01700E-004 2.1493443643114 46000E+003
P(4.0 s) = 6.3342483666239 95700E-005 1.5787192767323 96800E+004
P(4.5 s) = 6.7953462494774 18600E-006 1.4715953584806 70300E+005
P(5.0 s) = 5.7330314373604 80700E-007 1.7442778936869 13900E+006
P(5.5 s) = 3.7979124956066 81200E-008 2.6330253821191 82500E+007
P(6.0 s) = 1.9731752898266 56400E-009 5.0679734596101 19500E+008
P(6.5 s) = 8.0319972894926 65000E-011 1.2450203404677 24800E+010
P(7.0 s) = 2.5596191832732 98400E-012 3.9068311666627 59400E+011
P(7.5 s) = 6.3837823915946 50100E-014 1.5664694356071 29100E+013
P(8.0 s) = 1.2212453270876 72200E-015 8.1883629588554 47500E+014
P(8.5 s) = 0.0000000000000 00000E+000 1.#INF000000000 00000E+000
Same limits as in Excel, it seems. GSL's gsl_sf_erf function takes a
double-precision argument and returns double. My Visual C++ compiler (v.
6.0) specifies that both double and long double use an 8-byte
representation: "The long double contains 80 bits: 1 for sign, 15 for
exponent, and 64 for mantissa. Its range is +/- 1.2E4932 with at least 19
digits of precision."

Strangely enough, I only seem to be getting 17 digits of precision.
Regardless, I estimate that this calculation will require at least 38
additional digits of precision. Interesting: according to IEEE 754, the
condition for "positive underflow" (single precision) shouldn't happen
for positive numbers greater than about 1.4E-045 or so. For double
precision, it should happen only for positive numbers less than about
1.0E-308.

So here's what I'd like to know.

Are there 64-bit implementations of something like GSL which would
produce more precise output on appropriate OS/hardware platforms like
WinXP-64, Solaris v. 7-9, Tru64 Unix, Linux, etc? Has anyone
implemented a 128-bit long long double datatype or equivalent?

Can a calculation like this be performed with some kind of arbitrary-
precision numeric package (something like Michael Ring's MAPM or PHP
BCMath), or evaluated directly by Mathematica or the java.math package?

Or maybe I'm just confused and there's an easier way to do this which
I'm just not seeing.

By the way, I don't know of anyone who's bothered to tabulate the
values of this function nearly this far: most such texts (Zwillinger's CRC
Handbooks, Abramowitz & Stegun, Gradshteyn & Ryzhik, etc) only go to about
4.0 sigma.

Any ideas?

-- Dave Schulman (ca******@gate. net)

Nov 13 '05 #4

"DAVID SCHULMAN" <da************ *@verizon.net> wrote in message
news:t4******** *********@nwrdd c03.gnilink.net ...
I've been trying to perform a calculation that has been running into
an underflow (insufficient precision) problem in Microsoft Excel, which
calculates using at most 15 significant digits. For this purpose, that
isn't enough.

I was reading a book about some of the financial scandals of the
1990s called "Inventing Money: The Story of Long-Term Capital Management
and the Legends Behind it" by Nicholas Dunbar. On page 95, he mentions
that the 1987 stock-market crash was designated in economists' computer
models as a "20-sigma event". That is, their models (which obviously were
badly flawed!) put such an event in the exponential "tails" of a normal
Gaussian distribution, outside the range +/- 20 standard deviations from
the mean. In other words - vastly unlikely.

So I wanted to know: just how unlikely is that? Where did they get
the data to support that idea, anyway?

I put together an Excel spreadsheet, attempting to calculate this
probability. But after about 8 sigma I ran into the limits of Excel's
15-digit
precision and could go no further. Since the odds of an 8-sigma event
are about 819 trillion to 1 against, a 20-sigma event is obviously going to be so unlikely that it's not even worth talking about; would probably
never happen in the lifetime of the universe. Further evidence that those
computer models had some serious problems.

To perform that calculation, I used Excel's built-in error function
ERF(), and here's the output:

s Confidence Interval Probability Odds against

1.0 68.268949131685 700% 31.731050868314 3000000% 3
1.5 86.638554269553 900% 13.361445730446 1000000% 7
2.0 95.449972950730 900% 4.5500270492691 1000000% 22
2.5 98.758066814134 400% 1.2419331858655 6000000% 81
3.0 99.730020387427 800% 0.2699796125722 0200000% 370
3.5 99.953474183672 100% 0.0465258163279 4690000% 2,149
4.0 99.993665751560 700% 0.0063342484393 2140000% 15,787
4.5 99.999320465373 800% 0.0006795346262 3560100% 147,160
5.0 99.999942669685 300% 0.0000573303147 3997840% 1,744,278
5.5 99.999996202087 500% 0.0000037979124 9560668% 26,330,254
6.0 99.999999802682 500% 0.0000001973175 2898267% 506,797,346
6.5 99.999999991968 000% 0.0000000080319 9728949% 12,450,203,405
7.0 99.999999999744 000% 0.0000000002559 6191833% 390,683,116,666
7.5 99.999999999993 600% 0.0000000000063 8378239% 15,664,694,356, 071
8.0 99.999999999999 900% 0.0000000000000 0000000% 818,836,295,885 ,545
8.5 100.00000000000 000% 0.0000000000000 0000000% #DIV/0!
. . . .
s =ERF(An/SQRT(2)) =1-Bn =1/Cn

(The 'n' in the cell formulas above represents the row number).

Evidently this calculation hits the limit of Excel's computational
precision at about 8 sigma.

For what it's worth, ERF(z) is defined as 2/pi * INT[(0,z) exp(-t²)
dt], and the area under the "normal" Bell curve from -z to z is just
ERF(z/sqrt(2)). There's a discussion of it here:
http://mathworld.wolfram.com/Erf.html, here:
http://mathworld.wolfram.com/ConfidenceInterval.html, and here:
http://jove.prohosting.com/~skripty/page_295.htm.

So I decided to try to tackle this in C. I downloaded a package
called the GNU Scientific Library (http://www.gnu.org/software/gsl/) -
there's a nice precompiled binary for Microsoft Visual C++ at
http://www.network-theory.co.uk/gsl/freedownloads.html. I wrote a little
piece of code to try it out:

=============== =============== =============== =============== ==============
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_erf.h>

int main(void)
{
double odds;
double index;
double sigma;
double result;

for (index = 1.0; index <= 8.5; index += 0.5)
{
sigma = index / M_SQRT2;
result = 1 - (gsl_sf_erf (sigma));
odds = 1 / result;
printf("P(%2.1f \345) = %.18LE \t %-1.18lE\n", index, result, odds);
}
return 0;
}
=============== =============== =============== =============== ==============

And here's its output:

P(1.0 s) = 3.1731050786291 42600E-001 3.1514871875343 75500E+000
P(1.5 s) = 1.3361440253771 61700E-001 7.4842231152268 46800E+000
P(2.0 s) = 4.5500263896358 41700E-002 2.1977894507992 82900E+001
P(2.5 s) = 1.2419330651552 31800E-002 8.0519637334481 30300E+001
P(3.0 s) = 2.6997960632602 06900E-003 3.7039834734495 63900E+002
P(3.5 s) = 4.6525815807108 01700E-004 2.1493443643114 46000E+003
P(4.0 s) = 6.3342483666239 95700E-005 1.5787192767323 96800E+004
P(4.5 s) = 6.7953462494774 18600E-006 1.4715953584806 70300E+005
P(5.0 s) = 5.7330314373604 80700E-007 1.7442778936869 13900E+006
P(5.5 s) = 3.7979124956066 81200E-008 2.6330253821191 82500E+007
P(6.0 s) = 1.9731752898266 56400E-009 5.0679734596101 19500E+008
P(6.5 s) = 8.0319972894926 65000E-011 1.2450203404677 24800E+010
P(7.0 s) = 2.5596191832732 98400E-012 3.9068311666627 59400E+011
P(7.5 s) = 6.3837823915946 50100E-014 1.5664694356071 29100E+013
P(8.0 s) = 1.2212453270876 72200E-015 8.1883629588554 47500E+014
P(8.5 s) = 0.0000000000000 00000E+000 1.#INF000000000 00000E+000
Same limits as in Excel, it seems. GSL's gsl_sf_erf function takes a
double-precision argument and returns double. My Visual C++ compiler (v.
6.0) specifies that both double and long double use an 8-byte
representation: "The long double contains 80 bits: 1 for sign, 15 for
exponent, and 64 for mantissa. Its range is +/- 1.2E4932 with at least 19
digits of precision."

Strangely enough, I only seem to be getting 17 digits of precision.
Regardless, I estimate that this calculation will require at least 38
additional digits of precision. Interesting: according to IEEE 754, the
condition for "positive underflow" (single precision) shouldn't happen
for positive numbers greater than about 1.4E-045 or so. For double
precision, it should happen only for positive numbers less than about
1.0E-308.

So here's what I'd like to know.

Are there 64-bit implementations of something like GSL which would
produce more precise output on appropriate OS/hardware platforms like
WinXP-64, Solaris v. 7-9, Tru64 Unix, Linux, etc? Has anyone
implemented a 128-bit long long double datatype or equivalent?

Can a calculation like this be performed with some kind of arbitrary-
precision numeric package (something like Michael Ring's MAPM or PHP
BCMath), or evaluated directly by Mathematica or the java.math package?

Or maybe I'm just confused and there's an easier way to do this which
I'm just not seeing.

By the way, I don't know of anyone who's bothered to tabulate the
values of this function nearly this far: most such texts (Zwillinger's CRC
Handbooks, Abramowitz & Stegun, Gradshteyn & Ryzhik, etc) only go to about
4.0 sigma.

Any ideas?

-- Dave Schulman (ca******@gate. net)

&&&&&&&&&&&&&&& &&&&&&&&&&&&&&& &&&&&

Whatever method you used to compute those values of the normal
integral, they are not very accurate.

Here are the true values to 15 places, listed with those you posted:

True, to 15 places Your posted values
1.0 3.1731050786291 41e-01 3.1731050786291 42600E-001
1.5 1.3361440253771 61e-01 1.3361440253771 61700E-001
2.0 4.5500263896358 41e-02 4.5500263896358 41700E-002
2.5 1.2419330651552 27e-02 1.2419330651552 31800E-002
3.0 2.6997960632601 89e-03 2.6997960632602 06900E-003
3.5 4.6525815807105 01e-04 4.6525815807108 01700E-004
4.0 6.3342483666239 84e-05 6.3342483666239 95700E-005
4.5 6.7953462494601 21e-06 6.7953462494774 18600E-006
5.0 5.7330314375838 78e-07 5.7330314373604 80700E-007
5.5 3.7979124931775 44e-08 3.7979124956066 81200E-008
6.0 1.9731752900753 96e-09 1.9731752898266 56400E-009
6.5 8.0320011677182 36e-11 8.0319972894926 65000E-011
7.0 2.5596250877716 70e-12 2.5596191832732 98400E-012
7.5 6.3817833458217 92e-14 6.3837823915946 50100E-014
8.0 1.2441921148543 57e-15 1.2212453270876 72200E-015
8.5 1.8959069644406 64e-17 0.0000000000000 00000E+000
9.0 2.2571768119076 81e-19
9.5 2.0989030150725 21e-21
10.0 1.5239706048321 05e-23
10.5 8.6380126356184 61e-26
11.0 3.8213191489973 50e-28
11.5 1.3191542892226 00e-30
12.0 3.5529642242000 00e-33

Why most of those who deal with the normal integral in probability
theory are still stuck with the historical baggage of the error function
is a puzzle to me, as is the poor quality of the results one gets from
standard library implementations of erf(). (One of the most common
is based on ALGORITHM AS66, APPL. STATIST.(1973) Vol.22, .424 by HILL,
which gives only 6-8 digit accuracy).

Here is a listing of my method:

/*
Marsaglia Complementary Normal Distribution Function
cPhi(x) = integral from x to infinity of exp(-.5*t^2)/sqrt(2*pi), x<15
15-digit accuracy for x<15, returns 0 for x>15.
#include <math.h>
*/

double cPhi(double x){
long double v[]={0.,.655679542 41879847154L,
..4213692292880 5447322L,.30459 029871010329573 L,
..2366523829135 6067062L,.19280 810471531576488 L,
..1623776608968 6746182L,.14010 418345305024160 L,
..1231319632579 3229628L,.10978 728257830829123 L,
..9902859647173 1921395e-1L,.90175675501 064682280e-1L,
..8276628650136 9177252e-1L,.76475761016 248502993e-1L,
..7106958053885 2107091e-1L,.66374235823 250173591e-1L};
long double h,a,b,z,t,sum,p wr;
int i,j;
if(x>15.) return (0.);
if(x<-15.) return (1.);
j=fabs(x)+1.;
z=j;
h=fabs(x)-z;
a=v[j];
b=z*a-1.;
pwr=1.;
sum=a+h*b;
for(i=2;i<60;i+ =2){
a=(a+z*b)/i;
b=(b+z*a)/(i+1);
pwr=pwr*h*h;
t=sum;
sum=sum+pwr*(a+ h*b);
if(sum==t) break; }
sum=sum*exp(-.5*x*x-.91893853320467 274178L);
if(x<0.) sum=1.-sum;
return ((double) sum);
}
*/
end of listing
*/

The method is based on defining phi(x)=exp(-x^2)/sqrt(2pi) and

R(x)=cPhi(x)/phi(x).

The function R(x) is well-behaved and terms of its Taylor
series are readily obtained by a two-term recursion. With an accurate
representation of R(x) at ,say, x=0,1,2,...,15, a simple evaluation
of the Taylor series at intermediate points provides up to
15 digits of accuracy.
An article describing the method will be in the new version of
my Diehard CDROM. A new version of the Diehard tests
of randomness (but not yet the new DVDROM) is at
http://www.csis.hku.hk/~diehard/

One other point about your posting:

As is often done, you mistake odds for probabilities.
The odds for an event should be represented as the ratio
of p to 1-p, not 1/p. Thus if a bookie estimates the
probability that a certain horse will win as .2, then
the odds are 2 to 3 for and 3 to 2 against.
Of course when p is close to zero, the ratio 1/p is close
to (1-p)/p, but it is probably a sound practice to maintain
the distinction between odds and probabilities.

George Marsaglia
Nov 13 '05 #5

"George Marsaglia" <ge*@stat.fsu.e du> wrote in message news:wq******** ************@co mcast.com...


As is often done, you mistake odds for probabilities.
The odds for an event should be represented as the ratio
of p to 1-p, not 1/p. Thus if a bookie estimates the
probability that a certain horse will win as .2, then
the odds are 2 to 3 for and 3 to 2 against.
Of course when p is close to zero, the ratio 1/p is close
to (1-p)/p, but it is probably a sound practice to maintain
the distinction between odds and probabilities.


Typo? P = 0.2 = 1/5 <=> odds of 4-1 against.

John.
Nov 13 '05 #6

This thread has been closed and replies have been disabled. Please start a new discussion.

Similar topics

2
5434
by: Ronny Mandal | last post by:
Is there a function that will do this task properly? -- Thanks Ronny Mandal
19
18368
by: Allagappan | last post by:
Hello everyone, Just check this code: #include<stdio.h> int main() { double a,b; scanf("%lf",&a); scanf("%lf",&b); printf("\nNumber a: %0.16lf",a);
2
5875
by: mbelew | last post by:
I'm seeing a very strange behavior with double precision subtraction. I'm using csUnit for testing. If I run the test by itself, the test passes. When I run the batch of tests, the test fails. Here's the test: public void GetAcuteAngleDifference_PI() { double a = 0; double b = Math.PI;
2
1592
by: Sanjeev Azad | last post by:
I'm porting my application from VC++ 6.0 to VS .Net and experiencing a difference in the double precisions. I have the following code in my application double val = 16e-6; When I debug it, the debugger shows the value of val as "1.599999999999999e-5", not "1.6000000000000e-5", which is what the debugger shows in VC++6.0.
21
2677
by: Karl O. Pinc | last post by:
FYI, It'd be nice if the error message from a REFERENCES constraint mentioned the column name into which the bad data was attempted to be inserted. In PostgreSQL 7.3: sandbox=> insert into foo (id, b) values (3, 2); ERROR: b_is_fkey referential integrity violation - key referenced from
0
1966
by: Marcel Boscher | last post by:
Hello everybody, i get strange error messages when trying to call up my function with a SELECT functionname(with or without int); varying from: ERROR: function chr(double precision) does not exist ERROR: function FUNCTIONNAME() does not exist
6
5707
by: R.Biloti | last post by:
Hi folks I wrote the naive program to show up the unit roundoff (machine precision) for single and double precision: #include <stdio.h> int main (void) { double x;
4
3952
by: marathoner | last post by:
After reading the string representation of a double precision number, -0.417597000000000D+06, from a text file, I am unable to convert it to a double precision number. Here is the code I am using: string result = block.Split(charSeparators, StringSplitOptions.RemoveEmptyEntries); double temp1 = double.Parse(result); After the last...
0
2551
by: Charles Coldwell | last post by:
James Kanze <james.kanze@gmail.comwrites: True, with some additional considerations. The commonly used IEEE 754 floating point formats are single precision: 32 bits including 1 sign bit, 23 significand bits (with an implicit leading 1, for 24 total), and 8 exponent bits double precision: 64 bits including 1 sign bit, 52 significand...
0
7862
marktang
by: marktang | last post by:
ONU (Optical Network Unit) is one of the key components for providing high-speed Internet services. Its primary function is to act as an endpoint device located at the user's premises. However, people are often confused as to whether an ONU can Work As a Router. In this blog post, we’ll explore What is ONU, What Is Router, ONU & Router’s main...
0
7789
by: Hystou | last post by:
Most computers default to English, but sometimes we require a different language, especially when relocating. Forgot to request a specific language before your computer shipped? No problem! You can effortlessly switch the default language on Windows 10 without reinstalling. I'll walk you through it. First, let's disable language...
0
8144
Oralloy
by: Oralloy | last post by:
Hello folks, I am unable to find appropriate documentation on the type promotion of bit-fields when using the generalised comparison operator "<=>". The problem is that using the GNU compilers, it seems that the internal comparison operator "<=>" tries to promote arguments from unsigned to signed. This is as boiled down as I can make it. ...
0
8169
tracyyun
by: tracyyun | last post by:
Dear forum friends, With the development of smart home technology, a variety of wireless communication protocols have appeared on the market, such as Zigbee, Z-Wave, Wi-Fi, Bluetooth, etc. Each protocol has its own unique characteristics and advantages, but as a user who is planning to build a smart home system, I am a bit confused by the...
1
5670
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 1 May 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome a new presenter, Adolph Dupré who will be discussing some powerful techniques for using class modules. He will explain when you may want to use classes...
0
3803
by: TSSRALBI | last post by:
Hello I'm a network technician in training and I need your help. I am currently learning how to create and manage the different types of VPNs and I have a question about LAN-to-LAN VPNs. The last exercise I practiced was to create a LAN-to-LAN VPN between two Pfsense firewalls, by using IPSEC protocols. I succeeded, with both firewalls in...
1
2300
by: 6302768590 | last post by:
Hai team i want code for transfer the data from one system to another through IP address by using C# our system has to for every 5mins then we have to update the data what the data is updated we have to send another system
1
1400
muto222
by: muto222 | last post by:
How can i add a mobile payment intergratation into php mysql website.
0
1132
bsmnconsultancy
by: bsmnconsultancy | last post by:
In today's digital era, a well-designed website is crucial for businesses looking to succeed. Whether you're a small business owner or a large corporation in Toronto, having a strong online presence can significantly impact your brand's success. BSMN Consultancy, a leader in Website Development in Toronto offers valuable insights into creating...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.