Hi,

I am university student learning Bayesian statistics currently, and I am a first time user here. I would like to ask some support, please.

I intend to create a model for variable selection using Bayesian Lasso approach on the publicly available ‘wine quality data’. The data is attached, and it is also available by this link: https://archive.ics.uci.edu/ml/datasets/wine+quality.

I am using Metropolis-Hastings sampling implementing it with base R codes. I am aware there are other methods, but this is what I would like to use now.

My current aim is to have converging chains, and it is monitored by Gelman-Rubin statistic.

I used the following steps to reach convergences for all parameters:

1, I created the first variance-covariance (cov.prop) matrix with random values in all diagonal positions, and scaler (s=1) for this.

2, I run the sampler 3 times to get 1000 samples, and I checked the acceptance ratio in each cycle; I maintained the acceptance ratio by adjusting the ‘s’ scaling factor.

3, After the fourth run, I created a new variance-covariance matrix with the last 1000 samples.

4, I repeated step2 and step3 many times.

I expected the chains to converge; it happened with only the first 12 columns of ‘draws’ until reaching around 120000 samples then even those become unconverging again. I would like to find out why it is happening that (a) only the some parameters were converging, and (b) for what reasons even those became unconverging again.

The code I used:

- for the very first cov matrix:

A<- rnorm(25, 0.1, 0.01) # updated since the first post

cov.prop <- diag(A)

- for the updated cov matrix

int = c((step1-1000):step1)

M <- (draws[int, , 1] + draws[int, , 2] + draws[int, , 3] + draws[int, , 4])/4

cov.prop <- cov(M)

- for sampling:

d <- winequality_red

d <- as.data.frame(d)

y <- d$quality

X <- d[,-12]

X <- matrix(unlist(X), ncol = dim(X)[2], byrow = FALSE)

X <- cbind(rep(1,nrow(X)),scale(X))

n <- length(y)

k <- ncol(X)

npars <- k + (k-1) + 2

logprior <- function(pars) {

b <- as.matrix(pars[1:k])

tau <- pars[(k+1):(2*k-1)]
lambda2 <- pars[2*k]

sigma2 <- pars[npars]

if (lambda2 <= 0)

return(log(0))

if (sigma2 <= 0)

return(log(0))

if (any(tau<0))

return(log(0))

dens <- (sum(dnorm(abs(b[2:k]),0,sqrt(sigma2*tau),log=TRUE))+

dgamma(1/sigma2,4,60,log=TRUE) +

dgamma(lambda2,0.025,0.1,log=TRUE)+

sum(dexp(tau,lambda2/2,log=TRUE)))

return(dens)

}

loglikelihood <- function(y,X,par) {

b <- as.matrix(par[1:k])

sigma2 <- par[npars]

if (sigma2 <= 0)

return(log(0))

dens <- sum(dnorm(y,t(X%*%b),sqrt(sigma2),log=TRUE))

return(dens)

}

draws <- array(0, dim=c(200000, npars, 4))

draws[1,npars, 1:4] <- 1/rgamma(1,4,60)

draws[1,npars-1, 1:4] <- 1

draws[1,1, 1:4] <- mean(y)

draws[1,2:k, 1:4] <- rep(0,k-1)

draws[1,(k+1):(2*k-1), 1:4] <- rep(1,k-1)

converged <- c(1:npars)

converged[1:npars] <- 0

step1 <- 2

proposed <- c(1:npars)

proposed[1:npars] <- -1

accepted.prop <- rep(1, 4)

while (any(converged<1)) {

for (i in 1:4) {

proposed <- rmvnorm(1,draws[step1-1, , i], s*cov.prop)

r <- loglikelihood(y,X,proposed)+

logprior(proposed)-

loglikelihood(y,X,draws[step1-1, , i])-

logprior(draws[step1-1, , i])

u <- runif(1)

if (log(u) < r) {

draws[step1, , i] <- proposed

accepted.prop[i] <- accepted.prop[i]+1

} else {

draws[step1, , i] <- draws[step1-1, , i]

}

}

step1 <- step1 + 1

if (any(step1 == as.vector(seq(1000, 200000, by = 1000))) == TRUE) {

for ( z in 1:npars) {

chainlist <- mcmc.list(Chain1=mcmc(draws[1:(step1-1), z, 1]),

Chain2=mcmc(draws[1:(step1-1), z, 2]),

chain3=mcmc(draws[1:(step1-1), z, 3]),

chain4=mcmc(draws[1:(step1-1), z, 4]))

if ((gelman.diag(chainlist)$psrf[1,2]) < 1.4) {

converged[z] <- converged[z]+1

```
}
}
print(converged)
print((accepted.prop-accepted.prop.print))
accepted.prop.print <- accepted.prop
```

}

}