2.3.1.1 Asignment

Provide a Markov chain Monte Carlo algorithm to fit a zero-inflated Poisson distribution.

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=TRUE, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2 # Full conditional for cc

for(i in 1:n){
   v = rep(0,2)
   
   if(x[i]==0){
     v[1] = log(w)
     v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
     v    = exp(v - max(v))/sum(exp(v - max(v)))
     
    }else{
      v[1] = 0
      v[2] = 1
    }
   cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

w

## [1] 0.2

lambda

## Error in eval(expr, envir, enclos): object 'lambda' not found

Provide you estimates of the posterior means $\text{E}(\lambda \mid \mathbf{x})$ and $\text{E}(w\mid \mathbf{x})$.

mean(mu_out[burn:rrr])# = 3.051

## [1] 0

mean(w_out[burn:rrr])# = 0.398

## [1] 0

2.3.1.2 Marking

Are the parameters initialized in a reasonable way?

Correctly initializing the indicators $c_1,\dots ,c_n$ is key. In particular, observations such that $x_i\ne 0$ must have $c_i=2$. As for observations for which $x_i=0$, these can be randomly initialized to either component.

For the other parameters, it is reasonable to initialize $\lambda$ to the mean of the observations and the weight $\omega$ to a relatively small value (as we did in the EM algorithm).

Hence, the code for the initialization of the algorithm might look something like:

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=T, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2

However, be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the indicators $c_1,\dots ,c_n$ correct?

The indicators $c_1,\dots ,c_n$ are conditionally independent from each other, and we have

$Pr(c_i=1\mid \cdots )\propto \{\begin{array}{ll}w& x_i=0\\ 0& \text{otherwise}\end{array}$

while

$Pr(c_i=2\mid \cdots )\propto \{\begin{array}{ll}(1-w)\frac{{\lambda }^{x_i}{e}^{-\lambda }}{x_i!}& x_i=0\\ 1& \text{otherwise}\end{array}$

Hence, the code to sample the indicators might look something like this:

# Full conditional for cc
for(i in 1:n){
  v = rep(0,2)
  if(x[i]==0){
    v[1] = log(w)
    v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
    v    = exp(v - max(v))/sum(exp(v - max(v)))
  }else{
    v[1] = 0
    v[2] = 1
  }
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

Please be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the weight ww correct?

This is a simple one, as its structure is common to all mixture models. Recalling that the prior on $\omega$ is a uniform distribution on [0,1], we have

$\omega \mid \cdots \sim \text{Beta}(m(\mathbf{c})+1,n-m(\mathbf{c})+1)$

where $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1. The associated code might look something like this

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+n-sum(cc==1))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + n - sum(cc == 1)): object 'cc' not found

or, alternatively,

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

Please remember to be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the rate $\lambda$ correct?

Because we use a exponential prior on $\lambda$, the full conditional posterior is a Gamma distribution,

$\lambda \mid \cdots \sim \text{Gamma}(1+\sum _{i:c_i=2}{x}_i,1+n-m(\mathbf{c}))$

where, as before, $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1 (so that $n-m(\mathbf{c})$ is the number of observations assigned to component 2), and $\sum _{c_i=2}^i{x}_i$ is the sum of the observations assigned to component 2.

Hence, the code for this prompt might look something like:

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

• 0点 No
• 1点 Yes

Are the posterior means generated by the algorithm correct?

Recall that the posterior means can be approximated by simple averages of the posterior samples obtained after burn in. While there might a bit of Monte Carlo error in the answers reported, the values should be $E(\lambda \mid x)\approx 3.05$ and $E(\omega \mid x)\approx 0.40$. Note that these values are very close to those generated by the EM algorithm.

• 0点 No
• 1点 Yes

2.3.2.1 Asignment

Provide a Markov chain Monte Carlo algorithm to fit a zero-inflated Poisson distribution.

KK        = 2
w         = 1/2
mu        = mean(x)

## Error in mean(x): object 'x' not found

n         = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

aa        = rep(1,KK)  # Uniform prior on w
maxv      = 1 - 1e-6
# Number of iterations of the sampler
rrr       = 6000
burn      = 1000
w_out     = rep(0, rrr)
mu_out    = rep(0, rrr)

for(s in 1:rrr){
  # Sample the indicators
  cc      = rep(0,n)
  for(i in 1:n){
    v     = rep(0,KK)
    v[1]  = log(w)
    #v[1]  = log(w) + dunif(x[i], max=maxv, log=T) # treat first component as uniform (0, .9999)
    v[2]  = log(1-w) + dpois(x[i], mu, log=T)
    v     = exp(v - max(v))/sum(exp(v - max(v)))
    cc[i] = sample(1:KK, 1, replace=T, prob=v)
  }
  
  # Sample the weights
  w = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))

  # Sample the mean
  for(k in 2:KK){
    nk    = sum(cc==k)
    xsumk = sum(x[cc==k])
    mu    = rgamma(1, xsumk + 1, nk + 1)
  }

  # Store samples
  w_out[s]    = w
  mu_out[s]   = mu
}

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

Provide you estimates of the posterior means $\text{E}(\lambda \mid \mathbf{x})$ and $\text{E}(w\mid \mathbf{x})$.

mean(mu_out[burn:rrr])# = 3.051

## [1] 0

mean(w_out[burn:rrr])# = 0.398

## [1] 0

2.3.2.2 Marking

Are the parameters initialized in a reasonable way?

Correctly initializing the indicators $c_1,\dots ,c_n$ is key. In particular, observations such that $x_i\ne 0$ must have $c_i=2$. As for observations for which $x_i=0$, these can be randomly initialized to either component.

For the other parameters, it is reasonable to initialize $\lambda$ to the mean of the observations and the weight $\omega$ to a relatively small value (as we did in the EM algorithm).

Hence, the code for the initialization of the algorithm might look something like:

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=T, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2

However, be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the indicators $c_1,\dots ,c_n$ correct?

The indicators $c_1,\dots ,c_n$ are conditionally independent from each other, and we have

$Pr(c_i=1\mid \cdots )\propto \{\begin{array}{ll}w& x_i=0\\ 0& \text{otherwise}\end{array}$

while

$Pr(c_i=2\mid \cdots )\propto \{\begin{array}{ll}(1-w)\frac{{\lambda }^{x_i}{e}^{-\lambda }}{x_i!}& x_i=0\\ 1& \text{otherwise}\end{array}$

Hence, the code to sample the indicators might look something like this:

# Full conditional for cc
for(i in 1:n){
  v = rep(0,2)
  if(x[i]==0){
    v[1] = log(w)
    v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
    v    = exp(v - max(v))/sum(exp(v - max(v)))
  }else{
    v[1] = 0
    v[2] = 1
  }
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

Please be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the weight ww correct?

This is a simple one, as its structure is common to all mixture models. Recalling that the prior on $\omega$ is a uniform distribution on [0,1], we have

$\omega \mid \cdots \sim \text{Beta}(m(\mathbf{c})+1,n-m(\mathbf{c})+1)$

where $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1. The associated code might look something like this

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+n-sum(cc==1))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + n - sum(cc == 1)): object 'cc' not found

or, alternatively,

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

Please remember to be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the rate $\lambda$ correct?

Because we use a exponential prior on $\lambda$, the full conditional posterior is a Gamma distribution,

$\lambda \mid \cdots \sim \text{Gamma}(1+\sum _{i:c_i=2}{x}_i,1+n-m(\mathbf{c}))$

where, as before, $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1 (so that $n-m(\mathbf{c})$ is the number of observations assigned to component 2), and $\sum _{c_i=2}^i{x}_i$ is the sum of the observations assigned to component 2.

Hence, the code for this prompt might look something like:

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

• 0点 No
• 1点 Yes

Are the posterior means generated by the algorithm correct?

Recall that the posterior means can be approximated by simple averages of the posterior samples obtained after burn in. While there might a bit of Monte Carlo error in the answers reported, the values should be $E(\lambda \mid x)\approx 3.05$ and $E(\omega \mid x)\approx 0.40$. Note that these values are very close to those generated by the EM algorithm.

• 0点 No
• 1点 Yes

2.3.3.1 Asignment

Provide a Markov chain Monte Carlo algorithm to fit a zero-inflated Poisson distribution.

## Initialize the parameters
KK     = 2                               # number of components
n      = length(x)                         # number of samples

## Error in eval(expr, envir, enclos): object 'x' not found

v      = array(0, dim=c(n,KK))

## Error in array(0, dim = c(n, KK)): 'list' object cannot be coerced to type 'integer'

v[,1]  = 0.5*(x==0)                    #Assign half weight for nests with 0 eggs to first component

## Error in eval(expr, envir, enclos): object 'x' not found

v[,2]  = 1-v[,1]                              #Assign all of the remaining nests to the second component

## Error in eval(expr, envir, enclos): object 'v' not found

w      = mean(v[,1])                          #weight of the first component

## Error in mean(v[, 1]): object 'v' not found

lambda = sum(x*v[,2])/sum(v[,2])     #parameter (mean) of the second component

## Error in eval(expr, envir, enclos): object 'x' not found

## The actual MCMC algorithm starts here
# Priors
aa     = rep(1,KK)  # Uniform prior on w
alpha  = 2 #For Gamma distribution prior number of obs = alpha-1
beta   = 1 #Sum of observations = beta

# Number of iterations of the sampler
rrr    = 6000
burn   = 1000

# Storing the samples
cc.out      = array(0, dim=c(rrr, n))

## Error in array(0, dim = c(rrr, n)): 'list' object cannot be coerced to type 'integer'

w.out       = rep(0, rrr)
lambda.out  = array(0, dim=c(rrr, 1))
logpost     = rep(0, rrr)

# MCMC iterations
# MCMC iterations
for(s in 1:rrr){
  # Sample the indicators
  cc = rep(0,n)
  for(i in 1:n){
    v       = rep(0,KK)
    if (x[i]==0) {
      v[1]  = (w/(w+(1-w)*exp(-lambda)))
    } else {
      v[1]  = 0.0
    }
    v[2]    = 1.0 - v[1]
    cc[i]   = sample(1:KK, 1, replace=TRUE, prob=v)
  }
  
  # Sample the weights
  w         = rbeta(1, aa[1] + sum(cc==1), aa[2] + sum(cc==2))
  
  # Sample lambda
  nk        = sum(cc==2)
  xsumk     = sum(x[cc==2])
  alpha.hat = alpha + xsumk
  beta.hat  = beta + nk
  lambda    = rgamma(1, shape = alpha.hat, rate = beta.hat)
  
  # Store samples
  cc.out[s,]     = cc
  w.out[s]       = w
  lambda.out[s,] = lambda
  for(i in 1:n){
    if(cc[i]==1){
      logpost[s] = logpost[s] + log(w)
    }else{
      logpost[s] = logpost[s] + log(1-w) + dpois(x[i], lambda, log=TRUE)
    }
  }
  logpost[s] = logpost[s] + dbeta(w, aa[1], aa[2],log = T)
  logpost[s] = logpost[s] + dgamma(lambda, shape = alpha, rate = beta)
  if(s/500==floor(s/500)){
    print(paste("s =",s, w, lambda, logpost[s]))
  }
}

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

## Plot the logposterior distribution for various samples
par(mfrow=c(1,1))
par(mar=c(4,4,1,1)+0.1)
plot(logpost, type="l", xlab="Iterations", ylab="Log posterior")

w_avg = sum(w.out[(burn+1):rrr])/(rrr-burn)
lambda_avg = sum(lambda.out[(burn+1):rrr])/(rrr-burn)
print(paste(w_avg, lambda_avg))

## [1] "0 0"

w_sum = 0
lambda_sum = 0
counter = 0
for(s in 1:(rrr-burn)){
  w_sum = w_sum + w.out[s+burn]
  lambda_sum = lambda_sum + lambda.out[s+burn]
  counter = counter + 1
}
print(paste(counter, w_sum/counter, lambda_sum/counter))

## [1] "5000 0 0"

Provide you estimates of the posterior means $\text{E}(\lambda \mid \mathbf{x})$ and $\text{E}(w\mid \mathbf{x})$.

#E(lambda|x) = 3.056 and E(w|x) = 0.398
mean(lambda_sum)

## [1] 0

mean(w_sum)

## [1] 0

2.3.3.2 Marking

Are the parameters initialized in a reasonable way?

Correctly initializing the indicators $c_1,\dots ,c_n$ is key. In particular, observations such that $x_i\ne 0$ must have $c_i=2$. As for observations for which $x_i=0$, these can be randomly initialized to either component.

For the other parameters, it is reasonable to initialize $\lambda$ to the mean of the observations and the weight $\omega$ to a relatively small value (as we did in the EM algorithm).

Hence, the code for the initialization of the algorithm might look something like:

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=T, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2

However, be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the indicators $c_1,\dots ,c_n$ correct?

The indicators $c_1,\dots ,c_n$ are conditionally independent from each other, and we have

$Pr(c_i=1\mid \cdots )\propto \{\begin{array}{ll}w& x_i=0\\ 0& \text{otherwise}\end{array}$

while

$Pr(c_i=2\mid \cdots )\propto \{\begin{array}{ll}(1-w)\frac{{\lambda }^{x_i}{e}^{-\lambda }}{x_i!}& x_i=0\\ 1& \text{otherwise}\end{array}$

Hence, the code to sample the indicators might look something like this:

# Full conditional for cc
for(i in 1:n){
  v = rep(0,2)
  if(x[i]==0){
    v[1] = log(w)
    v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
    v    = exp(v - max(v))/sum(exp(v - max(v)))
  }else{
    v[1] = 0
    v[2] = 1
  }
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

Please be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the weight ww correct?

This is a simple one, as its structure is common to all mixture models. Recalling that the prior on $\omega$ is a uniform distribution on [0,1], we have

$\omega \mid \cdots \sim \text{Beta}(m(\mathbf{c})+1,n-m(\mathbf{c})+1)$

where $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1. The associated code might look something like this

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+n-sum(cc==1))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + n - sum(cc == 1)): object 'cc' not found

or, alternatively,

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

Please remember to be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the rate $\lambda$ correct?

Because we use a exponential prior on $\lambda$, the full conditional posterior is a Gamma distribution,

$\lambda \mid \cdots \sim \text{Gamma}(1+\sum _{i:c_i=2}{x}_i,1+n-m(\mathbf{c}))$

where, as before, $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1 (so that $n-m(\mathbf{c})$ is the number of observations assigned to component 2), and $\sum _{c_i=2}^i{x}_i$ is the sum of the observations assigned to component 2.

Hence, the code for this prompt might look something like:

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

• 0点 No
• 1点 Yes

Are the posterior means generated by the algorithm correct?

Recall that the posterior means can be approximated by simple averages of the posterior samples obtained after burn in. While there might a bit of Monte Carlo error in the answers reported, the values should be $E(\lambda \mid x)\approx 3.05$ and $E(\omega \mid x)\approx 0.40$. Note that these values are very close to those generated by the EM algorithm.

• 0点 No
• 1点 Yes

2.3.4.1 Asignment

Provide a Markov chain Monte Carlo algorithm to fit a zero-inflated Poisson distribution.

Provide you estimates of the posterior means $\text{E}(\lambda \mid \mathbf{x})$ and $\text{E}(w\mid \mathbf{x})$.

2.3.4.2 Marking

Are the parameters initialized in a reasonable way?

Correctly initializing the indicators $c_1,\dots ,c_n$ is key. In particular, observations such that $x_i\ne 0$ must have $c_i=2$. As for observations for which $x_i=0$, these can be randomly initialized to either component.

For the other parameters, it is reasonable to initialize $\lambda$ to the mean of the observations and the weight $\omega$ to a relatively small value (as we did in the EM algorithm).

Hence, the code for the initialization of the algorithm might look something like:

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=T, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2

However, be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the indicators $c_1,\dots ,c_n$ correct?

The indicators $c_1,\dots ,c_n$ are conditionally independent from each other, and we have

$Pr(c_i=1\mid \cdots )\propto \{\begin{array}{ll}w& x_i=0\\ 0& \text{otherwise}\end{array}$

while

$Pr(c_i=2\mid \cdots )\propto \{\begin{array}{ll}(1-w)\frac{{\lambda }^{x_i}{e}^{-\lambda }}{x_i!}& x_i=0\\ 1& \text{otherwise}\end{array}$

Hence, the code to sample the indicators might look something like this:

# Full conditional for cc
for(i in 1:n){
  v = rep(0,2)
  if(x[i]==0){
    v[1] = log(w)
    v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
    v    = exp(v - max(v))/sum(exp(v - max(v)))
  }else{
    v[1] = 0
    v[2] = 1
  }
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

Please be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the weight ww correct?

This is a simple one, as its structure is common to all mixture models. Recalling that the prior on $\omega$ is a uniform distribution on [0,1], we have

$\omega \mid \cdots \sim \text{Beta}(m(\mathbf{c})+1,n-m(\mathbf{c})+1)$

where $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1. The associated code might look something like this

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+n-sum(cc==1))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + n - sum(cc == 1)): object 'cc' not found

or, alternatively,

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

Please remember to be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the rate $\lambda$ correct?

Because we use a exponential prior on $\lambda$, the full conditional posterior is a Gamma distribution,

$\lambda \mid \cdots \sim \text{Gamma}(1+\sum _{i:c_i=2}{x}_i,1+n-m(\mathbf{c}))$

where, as before, $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1 (so that $n-m(\mathbf{c})$ is the number of observations assigned to component 2), and $\sum _{c_i=2}^i{x}_i$ is the sum of the observations assigned to component 2.

Hence, the code for this prompt might look something like:

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

• 0点 No
• 1点 Yes

Are the posterior means generated by the algorithm correct?

Recall that the posterior means can be approximated by simple averages of the posterior samples obtained after burn in. While there might a bit of Monte Carlo error in the answers reported, the values should be $E(\lambda \mid x)\approx 3.05$ and $E(\omega \mid x)\approx 0.40$. Note that these values are very close to those generated by the EM algorithm.

• 0点 No
• 1点 Yes

2.3.5.1 Asignment

Provide a Markov chain Monte Carlo algorithm to fit a zero-inflated Poisson distribution.

Provide you estimates of the posterior means $\text{E}(\lambda \mid \mathbf{x})$ and $\text{E}(w\mid \mathbf{x})$.

2.3.5.2 Marking

Are the parameters initialized in a reasonable way?

Correctly initializing the indicators $c_1,\dots ,c_n$ is key. In particular, observations such that $x_i\ne 0$ must have $c_i=2$. As for observations for which $x_i=0$, these can be randomly initialized to either component.

For the other parameters, it is reasonable to initialize $\lambda$ to the mean of the observations and the weight $\omega$ to a relatively small value (as we did in the EM algorithm).

Hence, the code for the initialization of the algorithm might look something like:

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=T, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2

However, be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the indicators $c_1,\dots ,c_n$ correct?

The indicators $c_1,\dots ,c_n$ are conditionally independent from each other, and we have

$Pr(c_i=1\mid \cdots )\propto \{\begin{array}{ll}w& x_i=0\\ 0& \text{otherwise}\end{array}$

while

$Pr(c_i=2\mid \cdots )\propto \{\begin{array}{ll}(1-w)\frac{{\lambda }^{x_i}{e}^{-\lambda }}{x_i!}& x_i=0\\ 1& \text{otherwise}\end{array}$

Hence, the code to sample the indicators might look something like this:

# Full conditional for cc
for(i in 1:n){
  v = rep(0,2)
  if(x[i]==0){
    v[1] = log(w)
    v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
    v    = exp(v - max(v))/sum(exp(v - max(v)))
  }else{
    v[1] = 0
    v[2] = 1
  }
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

Please be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the weight ww correct?

This is a simple one, as its structure is common to all mixture models. Recalling that the prior on $\omega$ is a uniform distribution on [0,1], we have

$\omega \mid \cdots \sim \text{Beta}(m(\mathbf{c})+1,n-m(\mathbf{c})+1)$

where $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1. The associated code might look something like this

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+n-sum(cc==1))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + n - sum(cc == 1)): object 'cc' not found

or, alternatively,

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

Please remember to be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the rate $\lambda$ correct?

Because we use a exponential prior on $\lambda$, the full conditional posterior is a Gamma distribution,

$\lambda \mid \cdots \sim \text{Gamma}(1+\sum _{i:c_i=2}{x}_i,1+n-m(\mathbf{c}))$

where, as before, $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1 (so that $n-m(\mathbf{c})$ is the number of observations assigned to component 2), and $\sum _{c_i=2}^i{x}_i$ is the sum of the observations assigned to component 2.

Hence, the code for this prompt might look something like:

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

• 0点 No
• 1点 Yes

Are the posterior means generated by the algorithm correct?

Recall that the posterior means can be approximated by simple averages of the posterior samples obtained after burn in. While there might a bit of Monte Carlo error in the answers reported, the values should be $E(\lambda \mid x)\approx 3.05$ and $E(\omega \mid x)\approx 0.40$. Note that these values are very close to those generated by the EM algorithm.

• 0点 No
• 1点 Yes

2.3.6.1 Asignment

Provide a Markov chain Monte Carlo algorithm to fit a zero-inflated Poisson distribution.

Provide you estimates of the posterior means $\text{E}(\lambda \mid \mathbf{x})$ and $\text{E}(w\mid \mathbf{x})$.

2.3.6.2 Marking

Are the parameters initialized in a reasonable way?

Correctly initializing the indicators $c_1,\dots ,c_n$ is key. In particular, observations such that $x_i\ne 0$ must have $c_i=2$. As for observations for which $x_i=0$, these can be randomly initialized to either component.

For the other parameters, it is reasonable to initialize $\lambda$ to the mean of the observations and the weight $\omega$ to a relatively small value (as we did in the EM algorithm).

Hence, the code for the initialization of the algorithm might look something like:

n        = length(x)

## Error in eval(expr, envir, enclos): object 'x' not found

cc       = rep(0, n)

## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'

cc[x==0] = sample(1:2, sum(x==0), replace=T, prob=c(1/2, 1/2))

## Error in sample.int(length(x), size, replace, prob): object 'x' not found

cc[x!=0] = 2

## Error in cc[x != 0] = 2: object 'cc' not found

lambda   = mean(x)

## Error in mean(x): object 'x' not found

w        = 0.2

However, be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the indicators $c_1,\dots ,c_n$ correct?

The indicators $c_1,\dots ,c_n$ are conditionally independent from each other, and we have

$Pr(c_i=1\mid \cdots )\propto \{\begin{array}{ll}w& x_i=0\\ 0& \text{otherwise}\end{array}$

while

$Pr(c_i=2\mid \cdots )\propto \{\begin{array}{ll}(1-w)\frac{{\lambda }^{x_i}{e}^{-\lambda }}{x_i!}& x_i=0\\ 1& \text{otherwise}\end{array}$

Hence, the code to sample the indicators might look something like this:

# Full conditional for cc
for(i in 1:n){
  v = rep(0,2)
  if(x[i]==0){
    v[1] = log(w)
    v[2] = log(1-w) + dpois(x[i], lambda, log=TRUE)
    v    = exp(v - max(v))/sum(exp(v - max(v)))
  }else{
    v[1] = 0
    v[2] = 1
  }
  cc[i] = sample(1:2, 1, replace=TRUE, prob=v)
}

## Error in 1:n: NA/NaN argument

Please be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the weight ww correct?

This is a simple one, as its structure is common to all mixture models. Recalling that the prior on $\omega$ is a uniform distribution on [0,1], we have

$\omega \mid \cdots \sim \text{Beta}(m(\mathbf{c})+1,n-m(\mathbf{c})+1)$

where $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1. The associated code might look something like this

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+n-sum(cc==1))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + n - sum(cc == 1)): object 'cc' not found

or, alternatively,

# Full conditional for w
w = rbeta(1, 1+sum(cc==1), 1+sum(cc==2))

## Error in rbeta(1, 1 + sum(cc == 1), 1 + sum(cc == 2)): object 'cc' not found

Please remember to be open-minded when reviewing this code, as there are many ways to accomplish the same thing in R.

• 0点 No
• 1点 Yes

Is the full conditional for the rate $\lambda$ correct?

Because we use a exponential prior on $\lambda$, the full conditional posterior is a Gamma distribution,

$\lambda \mid \cdots \sim \text{Gamma}(1+\sum _{i:c_i=2}{x}_i,1+n-m(\mathbf{c}))$

where, as before, $m(\mathbf{c})$ is the number of observations that $\mathbf{c}$ assigns to component 1 (so that $n-m(\mathbf{c})$ is the number of observations assigned to component 2), and $\sum _{c_i=2}^i{x}_i$ is the sum of the observations assigned to component 2.

Hence, the code for this prompt might look something like:

lambda = rgamma(1, sum(x[cc==2]) + 1, sum(cc==2) + 1)

## Error in rgamma(1, sum(x[cc == 2]) + 1, sum(cc == 2) + 1): object 'x' not found

• 0点 No
• 1点 Yes

Are the posterior means generated by the algorithm correct?

Recall that the posterior means can be approximated by simple averages of the posterior samples obtained after burn in. While there might a bit of Monte Carlo error in the answers reported, the values should be $E(\lambda \mid x)\approx 3.05$ and $E(\omega \mid x)\approx 0.40$. Note that these values are very close to those generated by the EM algorithm.

• 0点 No
• 1点 Yes