[R] how to improve the precison of this calculation?

[zhaoxing731]

Hello
	T
	I want to order some calculation "result", there will be lots of "result" that need to calculate and order
	PS: the "result" is just a intermediate varible and ordering them is the very aim

# problem:
# For fixed NT and CT, and some pair (c,n). order the pair by corresponding result
# c and n are both random variable

CT<-6000	#assignment to CT
NT<-29535210	#assignment to NT

# formulae for calculation intermediate variable "result", this is just a picture of the calculation which can't implement due to  the precision problem

i<-0:(c-1)
vec<-choose(n,i)*choose(NT-n,CT-i)	#get the vector which need summation
result<-log(sum(vec))			#get the log of summation

# thanks to Petr, we have a solution
calsta <- function(c, n) { 
i <- seq(from=0, length=c) 
logx <- lchoose(NT-n, CT-i) + lchoose(n, i) 
logmax <- max(logx) 
logmax + log(sum(exp(logx - logmax))) 
} 

# now, new problem arise, in theory, the "result" of different (c,n) pair should most probably differ, so I can order them, but
> calsta(918,100000)-calsta(718,100000)
[1] 0
> calsta(718,90000)-calsta(718,90000)
[1] 0
> calsta(718,90000)-calsta(918,100000)
[1] 0
# As Eik pointed out in another post, "calsta constantly returns 57003.6 for c>38 (the summands in
# sum(exp(logx - logmax)) will become 0 for c>38)" (when n= 10,000)


I think there are two ways to solve this problem:

1.Reducing the calcuation formulae algebraically get a conciser kernel for comparing. I think this is the perfect method and I failed many times, though this is not the topic of mailing-list, I hope if someone with stronge algebraical skills could give me some surprise

2.Through skills of R, which is difficult for me

PS: I don't want a perfect solution but a better one, which could have a higher discriminability than before

Thank you in advance

Yours sincerely
 
ZhaoXing
Department of Health Statistics
West China School of Public Health
Sichuan University
No.17 Section 3, South Renmin Road
Chengdu, Sichuan 610041
P.R.China	

__________________________________________________
¸Ï¿ì×¢²áÑÅ»¢³¬´óÈÝÁ¿Ãâ·ÑÓÊÏä?