[Rd] quantiles of the hypergeometric distribution (PR#502)

[Peter Dalgaard BSA]

Torsten Hothorn <hothorn@ci.tuwien.ac.at> writes:

>  term = exp(lfastchoose(NR, xr) + lfastchoose(NB, xb) - lfastchoose(N, n));
> 
> The expression inside the exp() function gives -1415.411856 and exp()
> returns 0. So term and sum in 
> 
>  while(sum < p && xr < xend) {
>         xr++;
>         NB++;
>         term *= (NR / xr) * (xb / NB);
>         sum += term;
>         xb--;
>         NR--;
>     }
> 
> are 0 all the time and xend = 600 is returned. 
> 
> R-Core: How to fix this?

Ick! That is of course only to be expected with recurrence relations
like this applied to large numbers. Since the term formula (as far as
I remember!) applies to all the point probabilities, I suppose that
one could start the recurrence from a "known positive" value and work
both forwards and backwards from there. Since you still want to do the
sums from the low end (although that is perhaps not all that
important?) you might need an intermediate array to hold the values.
One idea could be to start from the normal approximation of the
desired quantile. Probably a bit more than an afternoons work, but
hardly several weeks...

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel 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-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._