[R] recursive term
Hi r-users; I have this values: eign_val <- c(137.810447,3.538721,2.995161,1.685670) alp <- 1.6549 ; lamda <- eign_val lamda_m <- min(lamda) First I calculated manually: delta0 <- 1 delta1 <- alp*delta0*(4-lamda_m*(1/lamda[1]+1/lamda[2]+1/lamda[3]+1/lamda[4])) delta1 delta2 <- (alp/2)*(delta1*(delta1/alp) + delta0*((1-lamda_m/lamda[1])^2+ (1-lamda_m/lamda[2])^2+(1-lamda_m/lamda[3])^2+(1-lamda_m/lamda[4])^2)) delta2 delta3 <- (alp/3)*(delta2*(delta1/alp) + delta1*((1-lamda_m/lamda[1])^2+ (1-lamda_m/lamda[2])^2+(1-lamda_m/lamda[3])^2+(1-lamda_m/lamda[4])^2) + delta0*((1-lamda_m/lamda[1])^3+ (1-lamda_m/lamda[2])^3+(1-lamda_m/lamda[3])^3+(1-lamda_m/lamda[4])^3)) delta3 c(delta1,delta2,delta3) > c(delta1,delta2,delta3) [1] 3.224772 6.391966 10.091279 Then I wrote this the following code: term <- function(i, lamda, lamda_m) {sum(sapply(lamda, function(lamda, lamda_m,i) ((1-lamda_m/lamda)^i), lamda_m, i))} sm <- sapply(1:N, term, lamda, lamda_m);sm #now calculate the deltas k1 <- 3+1 delta <- rep(1, k1);delta delta_calc <- function(k1, delta, sm, alp) { k <- k1-1 alp/k*sum(sapply(1:k, function(i, delta, sm, k1) (sm[i]*delta[k1-i]), delta, sm, k1)) } delta[2:k1] <- sapply(2:k1, delta_calc, delta, sm, alp);delta[2:k1] > delta[2:k1] <- sapply(2:k1, delta_calc, delta, sm, alp);delta[2:k1] [1] 3.224772 2.804775 2.526796 I could not trace why I the answers are different. I hope somebody can help me with this. Thank you so much for your help. [[alternative HTML version deleted]] __ 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.
Re: [R] recursive term
you might try the following. Pochhammer_n <- function (a,b,c,n) { if(n==0) return(1) return(a*b*Pochhammer_n(a+1, b+1, c+1, n-1)/(c*n)) } hypergeo_sum <- function (a,b,c,z,n) { comb_sum <- 0 for (i in 0:n) { comb_sum <- comb_sum + Pochhammer_n(a,b,c,i)*z^i } return(comb_sum) } hypergeo_sum2 <- function (a,b,c,z,n) { comb <- rep(1,n+1) for (i in 2:(n+1)){ comb[i] <- comb[i-1]*(a+i-2)*(b+i-2)*z/((i-1)*(c+i-2)) } return(sum(comb)) } system.time(hypergeo_sum (1.25,1.75,1.25,0.5,300)) system.time(hypergeo_sum2(1.25,1.75,1.25,0.5,300)) On Dec 10, 4:28 pm, Roslina Zakaria wrote: > Hi, > > I would like to write a function for this formula: > F(a,b,c,z)=1+(ab)/(1!c)+ +(a(a+1)b(b+1))/(2!c(c+1))+ > +(a(a+1)(a+2)b(b+1)(b+2))/(3!c(c+1)(c+2))+… > > I wrote this function but not sure what is not right: > > hypergeo_sum <- function (a,b,c,z,n) > { for (i in 1:n) > { aa <- 1+(a*b*z)/c > aa[i] <- (aa-1)*(a+i)*(b+i)*z/((i+1)*(c+i)) > comb_sum <- aa + aa[i]+ aa[i+2] > } > comb_sum} > > > hypergeo_sum (1.25,1.75,1.25,0.5,3) > > output:> hypergeo_sum (1.25,1.75,1.25,0.5,3) > > [1] 2.394531 NA 1.039062 > > > The answer should be 2.852539. > > Or can you suggest any R book that I can refer to. Thank you so much for your > help. > > __ > r-h...@r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guidehttp://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. __ 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.
[R] recursive term
Hi, I would like to write a function for this formula: F(a,b,c,z)=1+(ab)/(1!c)+ +(a(a+1)b(b+1))/(2!c(c+1))+ +(a(a+1)(a+2)b(b+1)(b+2))/(3!c(c+1)(c+2))+… I wrote this function but not sure what is not right: hypergeo_sum <- function (a,b,c,z,n) { for (i in 1:n) { aa <- 1+(a*b*z)/c aa[i] <- (aa-1)*(a+i)*(b+i)*z/((i+1)*(c+i)) comb_sum <- aa + aa[i]+ aa[i+2] } comb_sum } hypergeo_sum (1.25,1.75,1.25,0.5,3) output: > hypergeo_sum (1.25,1.75,1.25,0.5,3) [1] 2.394531 NA 1.039062 The answer should be 2.852539. Or can you suggest any R book that I can refer to. Thank you so much for your help. __ 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.