Hello R community, following is my code and it shows error, can some one
fix this error and explain why this occurs?
gibbs <-function(m,n, theta = 0, lambda = 1){
alpha <- 1.5
beta <- 1.5
gamma <- 1.5
x<- array(0,c(m+1, 3))
x[1,1] <- theta
x[1,2] <- lambda
x[1,3]<- n
for(t in 2:m+1){
x[t,1] <- rbinom(x[t-1,3], 1, x[t-1,1])
x[t,2]<-rbeta(m, x[t-1,1] + alpha, tx[t-1,3] - x[t-1,1] + beta)
x[t,3] <- rpois(x[t-1,3] - x[t-1,1],(1 - x[t-1,2])*gamma)
}
x
}
gibbs(100, 10)
Gyn
[[alternative HTML version deleted]]
Hi, I see two problems right off: On Mon, Nov 7, 2011 at 3:10 PM, Gyanendra Pokharel <gyanendra.pokharel at gmail.com> wrote:> Hello R community, following is my code and it shows error, can some one > fix this error and explain why this occurs? > > gibbs <-function(m,n, theta = 0, lambda = 1){ > ? ?alpha <- 1.5 > ? ?beta <- 1.5 > ? ?gamma <- 1.5 > ? ?x<- array(0,c(m+1, 3)) > ? ?x[1,1] <- theta > ? ?x[1,2] <- lambda > ? ?x[1,3]<- n > ? ?for(t in 2:m+1){ > ? ? ? ?x[t,1] <- rbinom(x[t-1,3], 1, x[t-1,1])The rbinom() command here returns a vector of values, but you're trying to assign it to a single matrix element. You might want to double-check the help for rbinom() again.> ? ? ? ?x[t,2]<-rbeta(m, x[t-1,1] + alpha, tx[t-1,3] - x[t-1,1] + beta)tx isn't defined, but is probably a typo for x> ? ? ? ?x[t,3] <- rpois(x[t-1,3] - x[t-1,1],(1 - x[t-1,2])*gamma) > ? ?} > ? ?x > } > gibbs(100, 10)With a nice short function such as this, it's easy to set your arguments to their default values and actually run the function line by line at the command prompt so you can see whether what is happening is what you expect. Sarah -- Sarah Goslee http://www.functionaldiversity.org
The first argument to rbinom() is how many random samples you want to draw, not whatever you seem to think it is. It's not matching the size of what you mean to assign it to: in particular note that x[t-1, 3] is zero for t=3 which is where you initialize it. (I.e., you are also probably getting tripped up by an order of operations error) Michael On Mon, Nov 7, 2011 at 3:10 PM, Gyanendra Pokharel <gyanendra.pokharel at gmail.com> wrote:> Hello R community, following is my code and it shows error, can some one > fix this error and explain why this occurs? > > gibbs <-function(m,n, theta = 0, lambda = 1){ > ? ?alpha <- 1.5 > ? ?beta <- 1.5 > ? ?gamma <- 1.5 > ? ?x<- array(0,c(m+1, 3)) > ? ?x[1,1] <- theta > ? ?x[1,2] <- lambda > ? ?x[1,3]<- n > ? ?for(t in 2:m+1){ > ? ? ? ?x[t,1] <- rbinom(x[t-1,3], 1, x[t-1,1]) > ? ? ? ?x[t,2]<-rbeta(m, x[t-1,1] + alpha, tx[t-1,3] - x[t-1,1] + beta) > ? ? ? ?x[t,3] <- rpois(x[t-1,3] - x[t-1,1],(1 - x[t-1,2])*gamma) > ? ?} > ? ?x > } > gibbs(100, 10) > > Gyn > > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > 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. >
Hi: In your function call, x[1, 1] = theta = 0. In the first line of the loop, your rbinom() call works out to be x[2, 1] <- rbinom(x[1, 3], 1, x[1, 1]) <=> rbinom(10, 1, 0) That likely accounts for the error message: Error in x[t, 1] <- rbinom(x[t - 1, 3], 1, x[t - 1, 1]) : replacement has length zero HTH, Dennis On Mon, Nov 7, 2011 at 12:10 PM, Gyanendra Pokharel <gyanendra.pokharel at gmail.com> wrote:> Hello R community, following is my code and it shows error, can some one > fix this error and explain why this occurs? > > gibbs <-function(m,n, theta = 0, lambda = 1){ > ? ?alpha <- 1.5 > ? ?beta <- 1.5 > ? ?gamma <- 1.5 > ? ?x<- array(0,c(m+1, 3)) > ? ?x[1,1] <- theta > ? ?x[1,2] <- lambda > ? ?x[1,3]<- n > ? ?for(t in 2:m+1){ > ? ? ? ?x[t,1] <- rbinom(x[t-1,3], 1, x[t-1,1]) > ? ? ? ?x[t,2]<-rbeta(m, x[t-1,1] + alpha, tx[t-1,3] - x[t-1,1] + beta) > ? ? ? ?x[t,3] <- rpois(x[t-1,3] - x[t-1,1],(1 - x[t-1,2])*gamma) > ? ?} > ? ?x > } > gibbs(100, 10) > > Gyn > > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > 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. >