[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: tolerance in binopdf.m
From: |
Dr. Alexander Klein |
Subject: |
Re: tolerance in binopdf.m |
Date: |
Wed, 21 Sep 2011 09:08:44 +0200 |
Am 21.09.2011 um 07:23 schrieb Rik:
On 09/20/2011 08:50 PM, Ben Abbott wrote:
On Sep 20, 2011, at 9:03 PM, Ben Abbott wrote:
More info on this particular test. After some simplification ...
exp (gammaln (3)) * exp (2 * log (0.5)) - 0.5
ans = 0
exp (gammaln (3) + 2 * log (0.5)) - 0.5
ans = 1.1102e-16
Good morning,
I just tested this with 3.2.3 on OSX 10.5.8 which commes with lgamma
in math.h, and in this case, too, the result of lgamma is off by 0.5ulp.
I did some further investigation, first on the Mac:
##Set up some factorials, Y doesn't have enough digits to warrant any
rounding errors:
octave:37> X=1:15;Y=[1,cumprod(X)(1:end-1)];
##First let's see how well gamma fares
octave:38> Y-gamma(X)
ans =
Columns 1 through 12:
0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00
0.0000e+00 -2.9802e-08
Columns 13 through 15:
2.3842e-06 1.1444e-05 1.5259e-04
##That seems to be ok, even with the last few values the relative
error is always smaller than 1e-14.
##Now let's see how lgamma behaves (assuming that log does)
octave:39> log(Y)-lgamma(X)
ans =
Columns 1 through 12:
0.0000e+00 0.0000e+00 -1.1102e-16 -2.2204e-16 4.4409e-16
8.8818e-16 0.0000e+00 1.7764e-15 3.5527e-15 0.0000e+00
0.0000e+00 0.0000e+00
Columns 13 through 15:
-3.5527e-15 3.5527e-15 0.0000e+00
##No real issues here, either.
##Just to be sure compare log(gamma()) to lgamma()
octave:40> log(gamma(X))-lgamma(X)
ans =
Columns 1 through 12:
0.0000e+00 0.0000e+00 -1.1102e-16 -2.2204e-16 4.4409e-16
8.8818e-16 0.0000e+00 1.7764e-15 3.5527e-15 0.0000e+00
0.0000e+00 0.0000e+00
Columns 13 through 15:
-7.1054e-15 0.0000e+00 0.0000e+00
##One value becomes better, one becomes worse, all others remain the
same.
Now Xubuntu 10.04, also with 3.2.3 installed via Synaptic:
##Same setup as above:
octave:10> X=1:15;Y=[1,cumprod(X)(1:end-1)];
##Again, see how well gamma behaves
octave:11> Y-gamma(X)
ans =
Columns 1 through 9:
0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 -3.5527e-15
2.8422e-14 -1.1369e-13 -1.8190e-12 -2.1828e-11
Columns 10 through 15:
1.7462e-10 -2.3283e-09 -2.9802e-08 7.1526e-07 -1.0490e-05
1.5259e-04
##That's somewhat worse than the Mac version.
##But now for the real surprises:
octave:12> log(Y)-lgamma(X)
ans =
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
octave:13> log(gamma(X))-lgamma(X)
ans =
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Even with gamma having more roundoff error than on the mac, we get all
zeroes here. Very interesting!
After all, I think that none of the above discrepancies warrant any
concern. They're all in a range that is to be expected and may even be
caused my simply moving data in and out of x87 FPU registers that
provide some extra hidden bits of precision unless you force strict
IEEE754 compliance for an implementation.
I think Michaels issues are much more worthwile considering.
Kind regards,
Alex
P.s.: I was just asked to review a paper about roundoff error analysis
where the algorithms presented differed by several orders of magnitude
in their errors, and still none of them deserves to be labelled
"bad". ;)
--
Dr. Alexander Klein, Diplom-Mathematiker
Physiologisches Institut | TransMIT Zentrum
Raum 543 | für Numerische Methoden
Aulweg 129 | Heinrich-Buff-Ring 44
35392 Giessen | 35392 Giessen
- Re: tolerance in binopdf.m, (continued)
- Re: tolerance in binopdf.m, Rik, 2011/09/21
- Re: tolerance in binopdf.m, Marco Atzeri, 2011/09/21
- Re: tolerance in binopdf.m, Jordi Gutiérrez Hermoso, 2011/09/21
- Re: tolerance in binopdf.m, Michael D Godfrey, 2011/09/21
- Re: tolerance in binopdf.m, Michael D Godfrey, 2011/09/21
- Re: tolerance in binopdf.m, Ben Abbott, 2011/09/21
- Re: tolerance in binopdf.m, Ben Abbott, 2011/09/21
- Re: tolerance in binopdf.m, Marco atzeri, 2011/09/21
- Re: tolerance in binopdf.m, Ben Abbott, 2011/09/21
- Re: tolerance in binopdf.m, Jordi Gutiérrez Hermoso, 2011/09/21
- Re: tolerance in binopdf.m,
Dr. Alexander Klein <=
Re: Overhaul of statistical distribution functions, Michael D Godfrey, 2011/09/21
Re: Overhaul of statistical distribution functions, Michael D Godfrey, 2011/09/21