Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: large numbers in comb(n,k) function: no success


From   Maarten buis <maartenbuis@yahoo.co.uk>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: large numbers in comb(n,k) function: no success
Date   Thu, 5 Feb 2009 17:08:20 +0000 (GMT)

--- On Thu, 5/2/09, Inna Becher wrote:
> the formula I need is actually:
> 
> probability=1- [comb(n-m,k)/comb(n, k)]
> so I have useful numbers after the computation. (n=180000,
> m=2, k=2000)

                 (n-m)!
comb(n-m,k) = -----------
              k! (n-m-k)!

               n!
com(n,k) = ---------
           k! (n-k)!

so:
                            (n-m)!
                         -----------
                          k! (n-m-k)!
comb(n-m,k)/comb(n,k) = ---------------
                              n!
                         ---------
                          k! (n-k)!

                        (n-m)! k! (n-k)!
                     =  ----------------
                        k! (n-m-k)! n!

                        (n-m)! (n-k)!
                     =  -------------
                        (n-m-k)! n!

ln(comb(n-m,k)/comb(n,k)) = 
   ln((n-m)!) + ln((n-k)!) - ln((n-m-k)!) - ln(n!)

As comb(n-m,k)/comb(n,k) = exp(ln(comb(n-m,k)/comb(n,k))) you can compute your probability by typing:

scalar n = 180000
scalar m = 2
scalar k = 2000

di 1- exp(lnfactorial(n-m) + lnfactorial(n-k) - lnfactorial(n-m-k) - lnfactorial(n))

-- Maarten

-----------------------------------------
Maarten L. Buis
Department of Social Research Methodology
Vrije Universiteit Amsterdam
Boelelaan 1081
1081 HV Amsterdam
The Netherlands

visiting address:
Buitenveldertselaan 3 (Metropolitan), room N515

+31 20 5986715

http://home.fsw.vu.nl/m.buis/
-----------------------------------------







      

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index