1
$\begingroup$

I'm trying to implement the Metropolis Hastings algorithm in this problem but I'm having problems with the convergence.

$$Y_i|\beta_0,\beta_1 \sim \text{Binomial}(m_i,\theta_i)$$

where $logit(\theta) = \beta_0 + \beta_1x_i$.

$Y_i$ is the number of damage incidents of 6 possible (i.e number of successes). This also means $m_i = 6$.

I assume a normal distribution for the regression coefficients and I'm using random walk as a candidate distribution, so

$$\beta_0 \sim N(0,10)$$ $$\beta_0 \sim N(0,10)$$

$$\beta_j^*|\beta_{j}^{t-1} \sim N(\beta_{j}^{t-1},\sigma^2)$$

For the data I'm using

temp <- c(53,57,58,63,66,67,67,67,68,69,70,70,
          70,70,72,73,75,75,76,76,78,79,81)

success <- c(5,1,1,1,0,0,0,0,0,0,1,
             0,1,0,0,0,0,1,0,0,0,0,0)

failures <- c(1,5,5,5,6,6,6,6,6,6,
              5,6,5,6,6,6,6,5,6,6,6,6,6)

dat <- data.frame(success=success,failures=failures,temp=temp)
# Posterior distribution
post_beta <- function(Y,x,beta,m=6,mu_beta0=0,s2_beta0=10,mu_beta1=0,s2_beta1=10){
  b0 <- beta[1]
  b1 <- beta[2]
  pred <- b0+b1*x
  theta <- exp(pred)/(1+exp(pred))
  like <- dbinom(Y,m,prob = theta)
  prior <- dnorm(x=b0,mean = mu_beta0,sd = s2_beta0) *
    dnorm(x=b1,mean = mu_beta1,sd = s2_beta1)
  return(like*prior)
}

# Metropolis-Hastings

mh <- function(S,Y,x,init,fixed_sd){
  samples <- matrix(NA,S,2)
  colnames(samples) <- c("b0","b1")
  beta <- init
  for(s in 1:S){
    for(i in 1:length(beta)){
      can <- beta
      can[i] <- rnorm(1,mean = beta[i],fixed_sd)
      R <- min(1,post_beta(Y,x,can)/post_beta(Y,x,beta))
      if(runif(1) <= R){
        beta <- can
      }
    }
    samples[s,] <- c(beta)
  }
  return(samples)
}


S <- 50000 # Simulations
Y <- dat$success
x <- dat$temp
init <- c(11.66299,-0.21623) # Initial values (values from glm output)
fixed_sd <- 0.01 # Fixed sd for the M-H

res <- mh(S = S,Y = Y,x = x,init = init,fixed_sd = fixed_sd)
plot(res[,1],type = "l",main="b0")
plot(res[,2],type = "l",main="b1")

The results using a frequentist framework using glm are the following.

res_freq = glm(cbind(data_p2$success,data_p2$failures)~x,family = binomial())
summary(res_freq )

And I tried to use this result as initial values but that doesn't seem to work to either. I have the feeling that the mistake could be in post_beta function.

Thanks in advance!

$\endgroup$
3
  • $\begingroup$ this line R <- min(1,post_beta(Y,x,can)/post_beta(Y,x,beta)) is missing a correction factor I think $\endgroup$
    – Jellyfish
    Jun 24 at 2:40
  • $\begingroup$ post betta doesn't need the prior involved either because it cancels out in the acceptance probability. that being said, convergence is tough for these kinds of problems $\endgroup$
    – Jellyfish
    Jun 24 at 3:26
  • $\begingroup$ @Jellyfish thanks for the reply! What I don't know is if the likelihood is well defined in my post_beta function and the other thing is perhaps the standard deviation I'm using for the candidate is not the right one. $\endgroup$
    – Seb
    Jun 24 at 5:45

1 Answer 1

0
$\begingroup$

You can get nice grassy oscillation in the time series plots by changing the prior to

$$\begin{split}\beta_0&\sim \mathcal N(11.66299, .1^2)\\ \beta_1 &\sim \mathcal N(-0.21623, .1^2)\end{split}$$

and changing the standard deviation of your proposal/candidate to $.1$ for $\beta_0$ (keeping it at $.01$ for $\beta_1)$.

Lastly, the acceptance probability is missing the factor involving the candidate distribution.

temp <- c(53,57,58,63,66,67,67,67,68,69,70,70,
          70,70,72,73,75,75,76,76,78,79,81)

success <- c(5,1,1,1,0,0,0,0,0,0,1,
             0,1,0,0,0,0,1,0,0,0,0,0)

failures <- c(1,5,5,5,6,6,6,6,6,6,
              5,6,5,6,6,6,6,5,6,6,6,6,6)

dat <- data.frame(success=success,failures=failures,temp=temp)
# Posterior distribution
post_beta <- function(Y,x,beta,m=6,mu_beta0=11.66299,s2_beta0=.1,mu_beta1=-0.21623,s2_beta1=.1){
  b0 <- beta[1]
  b1 <- beta[2]
  pred <- b0+b1*x
  theta <- exp(pred)/(1+exp(pred))
  like <- dbinom(Y,m,prob = theta)
  prior <- dnorm(x=b0,mean = mu_beta0,sd = s2_beta0) *
    dnorm(x=b1,mean = mu_beta1,sd = s2_beta1)
  return(like*prior)
}

# Metropolis-Hastings

mh <- function(S,Y,x,init,fixed_sd){
  samples <- matrix(NA,S,2)
  colnames(samples) <- c("b0","b1")
  beta <- init
  for(s in 1:S){
    for(i in 1:length(beta)){
      can <- beta
      can[i] <- rnorm(1,mean = beta[i],fixed_sd[i])
      R <- min(1,post_beta(Y,x,can)/post_beta(Y,x,beta)*dnorm(beta[i],can[i],fixed_sd[i])/dnorm(can[i],beta[i],fixed_sd[i]))
      if(runif(1) <= R){
        beta <- can
      }
    }
    samples[s,] <- c(beta)
  }
  return(samples)
}


S <- 50000 # Simulations
Y <- dat$success
x <- dat$temp
init <- c(11.66299,-0.21623) # Initial values (values from glm output)
fixed_sd <- c(.1, 0.01) # Fixed sd for the M-H

res <- mh(S = S,Y = Y,x = x,init = init,fixed_sd = fixed_sd)
par(mfrow=c(2,2))
plot(res[,1],type = "l",main="b0")
plot(res[,2],type = "l",main="b1")
acf(res[,1])
acf(res[,2])

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks! On the factor involving the candidate distribution, I didn't consider it since I'm using a symmetric distribution for the candidate this line cancels out dnorm(beta[i],can[i],fixed_sd[i])/dnorm(can[i],beta[i],fixed_sd[i]). One more thing, how did you come up with those values for the standard deviation? $\endgroup$
    – Seb
    Jun 24 at 17:05
  • $\begingroup$ @Seb By trial and error. I chose priors close to the frequentist framework with small standard deviation. Then based on the magnitudes of the mean, I chose the standard deviation of the proposal. I think convergence might be easier if the model fit the data better, but I’m not sure, because using Bayesian software like JAGS would probably not require such a specific prior and still have good results. So, I’m not entirely sure either $\endgroup$
    – Jellyfish
    Jun 24 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.