求助：Gibbs算法的求解

如上图片所示。不知什么原因，用这些参数的满条件分布得到的模拟结果只有mu是正确的，sigma0和sigma不对，我想是否在sigma0和sigma中推导的满条件分布不对，但我又推导了一次还是这样，求大神解答。我附上我的代码：

1. Gib <- function(y, sigma0, sigma, mu, simnum){
   
2.   n = length(y)
   
3.   gibchain = matrix(NA, simnum, 3)
   
4.   gibchain[1,] = c(sigma0, sigma, mu)
   
5.   for (i in 2:simnum){
   
6.     currsigma0 = gibchain[i-1, 1]
   
7.     currsigma = gibchain[i-1, 2]
   
8.     currmu = gibchain[i-1, 3]
   
9.     ## sample from full conditional distribution of sigma0 (Gibbs)
   
10.     currsigma0 = 1/rgamma(1, shape = 2, scale = 1+currmu^2/2)
    
11.     currsigma = 1/rgamma(1, shape = 1+(n+6)/2, scale = 1+sum(y-currmu)^2/2)
    
12.     currmu = rnorm(1, mean = currsigma0*sum(y)/(n*currsigma0+2*currsigma),
    
13.                    sd = sqrt(currsigma0*currsigma))
    
14.     gibchain[i,] = c(currsigma0, currsigma, currmu)
    
15.   }
    
16.   return(gibchain)
    
17. }

复制代码

我用y1模拟数据作为实际数据

1. y1 <- rnorm(100, 0, sd = 4)
   
2. gibbs <- Gib(y = y1, sigma0 = 6, sigma = 2, mu = 0, simnum = 10000)

复制代码

除了mu正确，其他两个都不对。