Dear Eik

What a great idea!!! Thank you so much for your colossal improvment
Yes, you have a unique eye on the numerical problem, I am worrying about this problem right now, hope you could give me new idea again

Hi,
you compute the same results for logx many times. So it is easier and
time saving tabulating all intermediate results.
smth. like
 n<-100000
 CT=6000 #assignment to CT
 NT=29535210 #assignment to NT
 i <- 0:(n-1)
 lookup<- lchoose(NT-n, CT-i) + lchoose(n, i)
 lgmax<-cummax(lookup)
 calsta2<-function(c) lgmax[c] + log(sum(exp(lookup[1:c] - lgmax[c])))
should help for a start, but I think, you are running into numerical
troubles, since you are dealing with very high and low (on log scale)
numbers and calsta constantly returns 57003.6 for c>38 (the summands in
sum(exp(logx - logmax)) will become 0 for c>38).
#check
sapply(1:50,calsta2)
sapply(1:50,calsta)
hth
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

	[[alternative HTML version deleted]]

