On Sat, Jan 08, 2011 at 12:20:59AM +0800, zhaoxing731 wrote: > Hello > > I want to calculate natural logarithm of sum of combinations as follow: > (R code) > > { > > com_sum=choose(2000000,482)*choose(1000000,118)+choose(2000000,483)*choose(1000000,117)+...+choose(2000000,i)*choose(1000000,600-i)+...+choose(2000000,600)*choose(1000000,0) > #calculate the sum > result=log(com_sum) #calculate > the log of the sum > }
Let me present a more detailed calculation. The natural logarithms of the products above are i <- 482:600 logx <- lchoose(2000000,i) + lchoose(1000000,600-i) The logarithm of the largest of the products is maxlog <- max(logx) maxlog # [1] 5675.315 The ratios of the products divided by the largest of them may be computed as exp(logx - maxlog) # [1] 1.000000e+00 4.885522e-01 2.361712e-01 1.129583e-01 5.345075e-02 # ... Only a few of the ratios are not negligible, so the sum of the products differs from the largest product by a small factor. This factor may be estimated as sum(exp(logx - maxlog)) # [1] 1.937370 The required logarithm differs from maxlog by the logarithm of this factor. So, the result is options(digits=15) maxlog + log(sum(exp(logx - maxlog))) # [1] 5675.97681976982 Computing the sum of products exactly using long integer arithmetic and computing its logarithm with large precision produced 5675.9768197698224041850302683544102271115821956713424401649... which matches well the result obtained using double precision in R. The above approach may be derived from the numerical part of the algorithm suggested by Ted Harding for computing extremely small p-values in the Fisher Exact test. There, we are computing a sum of small numbers, not large ones. However, the problem is similar, since we cannot represent the summands themselves, but we can represent their logarithms. In my opinion, P(a, b, c, d) is the largest of the summands or close to the largest, so we calculate the ratios of the summands to the largest of them or to a close approximation of the largest. Petr Savicky. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.