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
__________________________________________________
¸Ï¿ì×¢²áÑÅ»¢³¬´óÈÝÁ¿Ãâ·ÑÓÊÏä?
______________________________________________
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.
