Gibbs Sampler for Normal and Inverse Gamma Distribution in R I'm trying to implement a Gibbs sampler for the following conditional distributions using R:

This is the code I have in R so far:

gibbs = function(N, v0, alpha, beta){

# Data
X = c(x1,...,xn);

# Parameters
n = length(X);
xbar = mean(X);
sigma2 = var(X);

# Starting values
mu0 = v0[1];
sigma20 = v0[2];
out = NULL;

for (i in 1:N){
  mu0 = rnorm(1, mean = ((n*xbar/sigma2 + mu0/sigma20) / (n/sigma2) + (1/sigma20)), sd = n/sigma2 + (1/sigma20)^-1);
  sigma20 = rinvgamma(1, n + alpha, 0.5*(sum(X - mu0))^2 + beta);
  out = rbind(out, c(mu0, sigma20)) 
}

return(out)

}


When I run the function using gibbs(1000, v0 = c(1,1), 0.001, 0.01), the code runs for 5 or so iterations before starting to output NaN for everything.
I'm new to R, and to MCMC methods, so any help would be greatly appreciated.
Thanks!
 A: Several conceptual and R-coding errors:

*

*$\mu_0$ and $\sigma_0$ are hyper-parameters for the prior, hence should not be modified along iterations


*the conditionals in the Gibbs sampler are about the sampling model parameters $\mu$ and $\sigma$ which means updating mu=rnorm(... and sigma2=rinvgamma(... in the R code


*there was a factor 2 missing in n/2 in the shape of the rinvgamma call


*The term sum(X - mu0))^2 is doubly wrong as it should be sum((X - mu)^2)


*the power -1 is not correctly set in the standard deviation in the  rnorm call, it should be
(n/sigma2 + 1/sigma20)^(-1/2)


*the rinvgamma function in MCMCpack is parameterised in terms of scale and shape, not rate and shape, hence the second parameter is the inverse of what it should be:
sigma2 = rinvgamma(1, n/2 + alpha, 1/(0.5*sum((X - mu)^2) + beta))
which is the principal reason for the diverging chain in the original version.
The entire R code could thus be
gibbs = function(N, v0, alpha, beta){    
# parameters
n = length(X)
nxbar = sum(X)
sigma2 = var(X)
# hyperparameters
mu0 = v0[1]
sigma20 = v0[2]
out = NULL   
for (i in 1:N){
  mu = rnorm(1, 
    mean = (nxbar/sigma2 + mu0/sigma20) / (n/sigma2) + (1/sigma20), 
    sd = 1/sqrt(n/sigma2 + 1/sigma20))
  sigma2 = rinvgamma(1, n/2 + alpha, 1/(sum((X - mu)^2)/2 + beta))
  out = rbind(out, c(mu, sigma2)) # yuk!!
}
return(out)}

