Message-Id: <3.0.32.19980223150857.007f7210@assist.stt.msu.edu> Date: Mon, 23 Feb 1998 15:08:58 -0500 To: r-help@stat.math.ethz.ch From: "Vincent F. Melfi" <melfi@assist.stt.msu.edu> Subject: R-beta: Hypergeometric Probabilities In both versions of R to which I currently have access (R-0.16.1 and R-0.61.1), "phyper" stops returning correct cumulative probabilities as the parameters of the hypergeometric distribution get large. For example, when N1=1345, N2=1055, and n=1330, phyper returns either 0 or 1, and nothing in between. Looking at phyper.c, it's clear what's happening. First a term (called "term") is computed: term = exp(lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n)); Then in a while loop, this is updated by multiplication, and each subsequent "term" is added: while(xr <= x) { sum += term; xr++; NB++; term *= (NR / xr) * (xb / NB); xb--; NR--; } If the first "term" is 0, the rest will be. This is the case in the example given above. One can get the correct cumulative probabilities by using "cumsum" and "dhyper", of course. Possibly it would be wise to modify "phyper" to fix this problem. (One easy modification would be to add an "if" statement to check on the parameters. If they're suitably small, use the current method. If they're larger, compute the probabilities as in "dhyper" and then add. A related question/comment: I assume similar problems may show up in other algorithms in src/math and elsewhere. Would it be wise to either modify the fucntions to return an error message in such cases, or to add something to the help indicating the range of applicability? Thanks. Vince Melfi PS: I'm by no means a programmer, so please be gentle if I've said anything stupid in the above. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request@stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._