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]]
______________________________________________
[email protected] 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.
