Re: binocdf accuracy

[Julien Bect]

On 06/07/2013 17:22, Dr. Alexander Klein wrote:
> Dear list,
>
> the other day I was writing up some statistical functions in Erlang, and kept 
> cross-checking the results with what Octave produced.
>
> One of my experiments included binocdf (0:50,50,1-1e-3 )', and while I admit 
> that the example is somewhat extreme, I suspect that almost all values returned 
> by Octave are inaccurate, while binopdf (0:50,50,1-1e-3 )' produces results 
> nearly identical to my implementation in Erlang, see attachments.
>
> I think this would be easy to patch if only binocdf was afflicted, but binocdf 
> uses the incomplete beta function, and so maybe other functions might have 
> problems, too.
>
> Since I'll need the incomplete beta function for my Erlang-library anyway, I 
> might as well try to come up with an improved version for Octave, but I'd like 
> to be sure that at least my reference values for binocdf are accurate.
>
> Is anyone on the list perhaps able to check my results against a reliable 
> implementation of binocdf?

Hi Alexander,

Here is what I get with R 3.0.0

> options(digits=16)
> as.matrix(pbinom(0:50, 50, 1-1e-3), 51, 1)
                        [,1]
 [1,] 1.000000000000048e-150
 [2,] 4.995100000000185e-146
 [3,] 1.222601176000084e-141
 [4,] 1.954248138157635e-137
 [5,] 2.293997233604237e-133
 [6,] 2.108416766146391e-129
 [7,] 1.579770225969158e-125
 [8,] 9.920310098220520e-122
 [9,] 5.326978264755308e-118
[10,] 2.483507479998758e-114
[11,] 1.017249991777891e-110
[12,] 3.695506847769541e-107
[13,] 1.199877960843249e-103
[14,] 3.503948975379120e-100
[15,]  9.251512155251822e-97
[16,]  2.218227598752389e-93
[17,]  4.847716982481642e-90
[18,]  9.686154229238391e-87
[19,]  1.774099903142615e-83
[20,]  2.985114152841139e-80
[21,]  4.622538019124792e-77
[22,]  6.597384469543245e-74
[23,]  8.688365835296284e-71
[24,]  1.056722867481452e-67
[25,]  1.187704671842907e-64
[26,]  1.234067457739350e-61
[27,]  1.185509907183139e-58
[28,]  1.052822455583338e-55
[29,]  8.640336117699031e-53
[30,]  6.548844400790006e-50
[31,]  4.580113546826153e-47
[32,]  2.952316331432568e-44
[33,]  1.751420477321750e-41
[34,]  9.545076165660551e-39
[35,]  4.768562275628828e-36
[36,]  2.178143731310471e-33
[37,]  9.068455578420740e-31
[38,]  3.428711280917076e-28
[39,]  1.172136642140733e-25
[40,]  3.604149744664342e-23
[41,]  9.905340359574533e-21
[42,]  2.414644258229275e-18
[43,]  5.172005397479049e-16
[44,]  9.619554651075891e-14
[45,]  1.530252903231261e-11
[46,]  2.040786149341117e-09
[47,]  2.219818879633079e-07
[48,]  1.892165514006139e-05
[49,]  1.186482503817934e-03
[50,]  4.879437180296870e-02
[51,]  1.000000000000000e+00

and here is what I get with Matlab R2012a + Statistics toolbox 8.0 :

>> fprintf('%.15e\n', binocdf (0:50,50,1-1e-3 )')
1.000000000000069e-150
4.995100000000349e-146
1.222601176000022e-141
1.954248138157661e-137
2.293997233604027e-133
2.108416766146376e-129
1.579770225969158e-125
9.920310098220262e-122
5.326978264755390e-118
2.483507479998575e-114
1.017249991777909e-110
3.695506847769498e-107
1.199877960843271e-103
3.503948975379016e-100
9.251512155251722e-97
2.218227598752328e-93
4.847716982481582e-90
9.686154229238118e-87
1.774099903142560e-83
2.985114152841078e-80
4.622538019124635e-77
6.597384469543131e-74
8.688365835296052e-71
1.056722867481404e-67
1.187704671842860e-64
1.234067457739307e-61
1.185509907183136e-58
1.052822455583311e-55
8.640336117698896e-53
6.548844400790120e-50
4.580113546826244e-47
2.952316331432627e-44
1.751420477321754e-41
9.545076165660416e-39
4.768562275628819e-36
2.178143731310471e-33
9.068455578420579e-31
3.428711280917081e-28
1.172136642140745e-25
3.604149744664332e-23
9.905340359574605e-21
2.414644258229257e-18
5.172005397479073e-16
9.619554651075910e-14
1.530252903231260e-11
2.040786149341111e-09
2.219818879633077e-07
1.892165514006138e-05
1.186482503817935e-03
4.879437180296867e-02
1.000000000000000e+00

@++
Julien