R help - Feb 2011 - How can we make a vector call a function element-wise efficiently?

Hello
	I have a time-comsuming program which need to simplify, I have tested the
annotated program as follow:
> #define function which will be call
> calsta <- function(c, n=100000) 
+ { 
+   i <- seq(from=0, length=c) 
+   logx <- lchoose(NT-n, CT-i) + lchoose(n, i) 
+   logmax <- max(logx) 
+   logmax + log(sum(exp(logx - logmax))) 
+ } 
> CT=6000	#assignment to CT
> NT=29535210	#assignment to NT
> 
> vec<-c(2331,524,918,218,1100,547,289,1167,450,1723)
> vec
 [1] 2331  524  918  218 1100  547  289 1167  450 1723
> vec<-rep(vec,1000)#replicate the vec 1000 times
> length(vec)
[1] 10000 
> #then I'd like to make vector "vec" call function calsta
element-wise
> #and save the output to vector "result"
> system.time(result<-sapply(vec,calsta))
   user  system elapsed 
  26.45    0.03   26.70 
> 
> system.time(for (i in 1:10000) result[i]=calsta(vec[i]))
   user  system elapsed 
  27.74    0.14   28.94 

I have about  300,000 such 26.70/ 28.94 seconds, so the approximate computation
time is 100 days
What a terrible thing to do!!!
Any modification, nomatter how subtle, will be a big help for me

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

__________________________________________________
?????????????????????????????

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

Am 11.02.2011 06:12, schrieb zhaoxing731:
> Hello
> 	I have a time-comsuming program which need to simplify, I have tested the
annotated program as follow:
> 
>> #define function which will be call
> 
>> calsta <- function(c, n=100000) 
> + { 
> +   i <- seq(from=0, length=c) 
> +   logx <- lchoose(NT-n, CT-i) + lchoose(n, i) 
> +   logmax <- max(logx) 
> +   logmax + log(sum(exp(logx - logmax))) 
> + } 
>> CT=6000	#assignment to CT
>> NT=29535210	#assignment to NT
>>
>> vec<-c(2331,524,918,218,1100,547,289,1167,450,1723)
>> vec
>  [1] 2331  524  918  218 1100  547  289 1167  450 1723
>> vec<-rep(vec,1000)#replicate the vec 1000 times
>> length(vec)
> [1] 10000 
> 
>> #then I'd like to make vector "vec" call function calsta
element-wise
>> #and save the output to vector "result"
> 
>> system.time(result<-sapply(vec,calsta))
>    user  system elapsed 
>   26.45    0.03   26.70 
>>
>> system.time(for (i in 1:10000) result[i]=calsta(vec[i]))
>    user  system elapsed 
>   27.74    0.14   28.94 
> 
> I have about  300,000 such 26.70/ 28.94 seconds, so the approximate
computation time is 100 days
> What a terrible thing to do!!!
> Any modification, nomatter how subtle, will be a big help for me
> 
> 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 at 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.

-- 
Eik Vettorazzi
Institut f?r Medizinische Biometrie und Epidemiologie
Universit?tsklinikum Hamburg-Eppendorf

Martinistr. 52
20246 Hamburg

T ++49/40/7410-58243
F ++49/40/7410-57790

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]]

Hi ZhaoXing,
without knowledge about the ultimate purpose of your calculations its
quite difficult to give further hints.

best regards

Am 11.02.2011 14:35, schrieb zhaoxing731:
> 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/

-- 
Eik Vettorazzi
Institut f?r Medizinische Biometrie und Epidemiologie
Universit?tsklinikum Hamburg-Eppendorf

Martinistr. 52
20246 Hamburg

T ++49/40/7410-58243
F ++49/40/7410-57790