如何指定 Wishart 分布尺度矩阵

机器算法验证 r 事先的 锯齿 愿望分配
2022-04-08 02:31:20

我正在使用 rjags 包在 R 中运行以下贝叶斯混合模型，但我很难为 Wishart 分布指定比例矩阵。本质上，我希望 Sigma.inv 是一个协方差矩阵，它解释了从对照样本的重复分析中确定的仪器测量误差。我之前选择了 Wishart 分布作为 Sigma.inv，尺度矩阵采用控制样本的协方差矩阵的形式。然而，这导致模型无法收敛，即使在 100 万次运行和 500,000 次老化后也是如此和 125 细化长度，并且对 p 的估计产生微小的误差。

所以我的问题是比例矩阵应该采用不同的形式吗？或者在这种情况下是否有更合适的 Wishart 分布替代方案？我应该补充一点，如果使用单位矩阵来缩放 Wishart 分布，该模型效果很好，但这当然并不代表我想要合并的测量误差。

任何帮助将不胜感激。

model
{
 # Likelihood
for(i in 1:N) {
Y[i,1:J] ~ dmnorm(mu[i,],Sigma.inv)
mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]
}

# Priors on s
for(i in 1:N) {
for(k in 1:K) {
  s[k,1:J,i] ~ dmnorm(mu.s[,k],Sigma.s.inv[,,k])
  }
}

# Prior on p
for(i in 1:N) {
p[i,1:K] <- expphiV[i,1:K]/sum(expphiV[i,1:K])
phiV[i,1:K] <- phi[i,1:(K-1)]%*%tV
for(k in 1:K) {
  expphiV[i,k] <- exp(phiV[i,k])
}
}

# Prior on phi
for(i in 1:N) {
  for(k in 1:(K-1)) {
  phi[i,k] ~ dnorm(0,kappa.inv[k])
}
philast[i] <- -sum(phi[i,1:(K-1)])
}

for(k in 1:(K-1)) {
kappa.inv[k] ~ dgamma(2,1)
kappa[k] <- 1/kappa.inv[k]
}

# Prior on Sigma
Sigma.inv ~ dwish(R.Sigma,k.Sigma)
}

#

运行模型的 R 代码

#

# Y values
Y<-matrix(c(5.8565, 18.5979,    0.0045, 5.7579, 1.1401, 0.4945, 0.2122, 0.2613, 0.3984,
          14.0696,  10.6925,    0.0084, 8.5126, 2.1779, 1.1406, 0.4051, 0.3472, 0.5971,
        13.7417,    10.0652,    0.0084, 8.2977, 2.2172, 1.1204, 0.3929, 0.3583, 0.6143,
        11.7867,    15.1036,    0.0072, 8.3222, 1.9592, 0.9697, 0.3454, 0.3067, 0.5695),
      byrow=TRUE, ncol=9)
colnames(Y)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")

# mu.s values
mu.s<-matrix(c(6.97463333,  6.657444444,    10.560500000,   13.944108333,
          35.47273333,  13.221666667,   6.794500000,    4.149395833,
          0.00358000,   0.004533333,    0.008404167,    0.008891667,
          5.04523333,   7.128111111,    6.149708333,    6.932004167,
          1.19620000,   1.121222222,    2.086083333,    2.446995833,
          0.61910000,   0.485333333,    1.010833333,    0.880383333,
          0.19653333,   0.199222222,    0.491125000,    0.406904167,
          0.07113333,   0.312111111,    0.339916667,    0.298254167,
          0.45840000,   0.307222222,    0.616708333,    0.658416667),
         byrow=TRUE, ncol=4)
rownames(mu.s)<-c("Al","Ca","Ce","Fe","K","Mg","Na","P","Ti")
colnames(mu.s)<-c("CB","FD","RV","TS")

# Covariance matrix accounting for measurement error of control samples 
Sigma.inv<-structure(c(7.516889e-03,  4.719111e-04, 5.117778e-06,   4.896822e-03,   2.683444e-03,   5.368000e-04,   4.266667e-05,   -8.364444e-05,  4.630667e-04,
                  4.719111e-04, 8.097005e-05,   1.132174e-06,   8.006300e-04,   3.065903e-04,   4.421449e-05,   -4.869565e-06,  -5.677295e-06,  7.178551e-05,
                  5.117778e-06, 1.132174e-06,   9.927053e-08,   1.090671e-05,   4.145217e-06,   2.788406e-07,   -4.459903e-07,  -2.015942e-07,  9.478261e-07,
                  4.896822e-03, 8.006300e-04,   1.090671e-05,   9.350603e-03,   3.260499e-03,   4.947188e-04,   -3.441932e-05,  -4.135604e-05,  7.779478e-04,
                  2.683444e-03, 3.065903e-04,   4.145217e-06,   3.260499e-03,   1.390332e-03,   2.255014e-04,   -9.909179e-06,  -3.122995e-05,  2.906319e-04,
                  5.368000e-04, 4.421449e-05,   2.788406e-07,   4.947188e-04,   2.255014e-04,   7.714783e-05,   1.326377e-05,   -6.255072e-06,  4.525217e-05,
                  4.266667e-05, -4.869565e-06,  -4.459903e-07,  -3.441932e-05,  -9.909179e-06,  1.326377e-05,   1.484058e-05,   2.859903e-06,   -3.663768e-06,
                  -8.364444e-05,    -5.677295e-06,  -2.015942e-07,  -4.135604e-05,  -3.122995e-05,  -6.255072e-06,  2.859903e-06,   4.060386e-06,   -5.611594e-06,
                  4.630667e-04, 7.178551e-05,   9.478261e-07,   7.779478e-04,   2.906319e-04,   4.525217e-05,   -3.663768e-06,  -5.611594e-06,  7.074783e-05), dim=c(9,9))
Sigma.inv<-as.positive.definite(Sigma.inv)


# Sigma.s.inv values - a precision matrix
Sigma.s.inv<-structure(c(53.2149530065154, 7.78555807206446, -23455.5492262134, 
                     3.98226035437711, 181.147593576025, -554.144707627592, 810.959913329385, 
                     -222.049892052889, -676.160229178196, 7.78555807206446, 1.66555095380872, 
                     -2875.24335449319, 0.59057208979672, 35.2056764704436, -90.5513686454713, 
                     122.419754836375, -24.1464188840069, -110.985869289511, -23455.5492262134, 
                     -2875.24335449319, 26623099.2220065, -551.965187580738, -109728.586986558, 
                     276663.105767497, -386431.378216831, -82170.2904382779, 259038.414073531, 
                     3.98226035437711, 0.59057208979672, -551.965187580738, 0.918694682125983, 
                     9.31656679803296, -36.0723479510435, 38.501457515572, -39.2951967395694, 
                     -43.3516097843844, 181.147593576025, 35.2056764704436, -109728.586986558, 
                     9.31656679803296, 1046.08776695798, -2268.37491730957, 2936.52705816089, 
                     60.3905421599093, -2791.52606704426, -554.144707627592, -90.5513686454713, 
                     276663.105767497, -36.0723479510435, -2268.37491730957, 6336.11679874626, 
                     -8941.85134505843, 1263.96712724819, 7151.07272932156, 810.959913329385, 
                     122.419754836375, -386431.378216831, 38.501457515572, 2936.52705816089, 
                     -8941.85134505843, 14092.5591797243, -2587.63879760624, -10498.261488095, 
                     -222.049892052889, -24.1464188840069, -82170.2904382779, -39.2951967395694, 
                     60.3905421599093, 1263.96712724819, -2587.63879760624, 5072.16198995712, 
                     2677.23974494659, -676.160229178196, -110.985869289511, 259038.414073531, 
                     -43.3516097843844, -2791.52606704426, 7151.07272932156, -10498.261488095, 
                     2677.23974494659, 10377.5013426694, 39071331365.1695, -857206689.09299, 
                     13250725351142.8, 4293968882.52526, -138526471018.157, -346806050176.474, 
                     -185731096148.676, -40154997287.1756, -13832086879.5978, -857206689.09299, 
                     45035996.2737049, -291049796239.663, -94316253.2767046, 3042708991.03323, 
                     7617532441.71197, 4079550081.74548, 881997284.581011, 303819297.524648, 
                     13250725351142.8, -291049796239.663, 4499056063320832, 1457943322445.14, 
                     -47034282019304.7, -117752032879430, -63061801053651.4, -13633939080436.7, 
                     -4696447331898.41, 4293968882.52526, -94316253.2767046, 1457943322445.14, 
                     513861374.526114, -15241712287.5009, -38158180147.8331, -20435516111.352, 
                     -4418151354.45332, -1521909040.23616, -138526471018.157, 3042708991.03323, 
                     -47034282019304.7, -15241712287.5009, 491753442721.598, 1231009860799.41, 
                     659264193054.762, 142532685331.566, 49097861284.5838, -346806050176.474, 
                     7617532441.71197, -117752032879430, -38158180147.8331, 1231009860799.41, 
                     3081922920585.39, 1650491845606.43, 356835753178.233, 122918278499.906, 
                     -185731096148.676, 4079550081.74548, -63061801053651.4, -20435516111.352, 
                     659264193054.762, 1650491845606.43, 883961734443.138, 191102478025.129, 
                     65828570726.718, -40154997287.1756, 881997284.581011, -13633939080436.7, 
                     -4418151354.45332, 142532685331.566, 356835753178.233, 191102478025.129, 
                     41361275161.6651, 14232113705.0418, -13832086879.5978, 303819297.524648, 
                     -4696447331898.41, -1521909040.23616, 49097861284.5838, 122918278499.906, 
                     65828570726.718, 14232113705.0418, 4947128692.00388, 123.244158717727, 
                     15.1472733328015, -29469.8213424378, -162.933686694468, -830.16396214664, 
                     94.7853096969282, -68.8560246643179, 741.612274164014, 2264.45755928194, 
                     15.1472733328015, 5.35479274130578, -399.822164505876, -26.2908447667995, 
                     -48.3456870303433, -20.9894871349876, -7.0699744151633, 67.0416328725067, 
                     116.600590769204, -29469.8213424378, -399.822164505876, 25454078.4009846, 
                     34253.206150171, 238183.903396863, -33771.2583884048, 38523.5333162132, 
                     -288045.203177399, -1016390.47602203, -162.933686694468, -26.2908447667995, 
                     34253.206150171, 237.878595765458, 986.561801969833, -68.5166440497585, 
                     209.1509212556, -915.217874817308, -2932.38505457965, -830.16396214664, 
                     -48.3456870303433, 238183.903396863, 986.561801969833, 6669.29110276054, 
                     -1078.91990285319, 352.917095152379, -5609.73421176863, -17946.4520240348, 
                     94.7853096969282, -20.9894871349876, -33771.2583884048, -68.5166440497585, 
                     -1078.91990285319, 1052.57768971324, -913.546632402807, 484.410845136349, 
                     2431.79560931122, -68.8560246643179, -7.0699744151633, 38523.5333162132, 
                     209.1509212556, 352.917095152379, -913.546632402807, 2990.61062774375, 
                     525.038417043516, -3685.88997836577, 741.612274164014, 67.0416328725067, 
                     -288045.203177399, -915.217874817308, -5609.73421176863, 484.410845136349, 
                     525.038417043516, 7256.45767695138, 17207.7716622753, 2264.45755928194, 
                     116.600590769204, -1016390.47602203, -2932.38505457965, -17946.4520240348, 
                     2431.79560931122, -3685.88997836577, 17207.7716622753, 62661.2158837379, 
                     7.0546793436162, -0.0754213178285816, -2357.02875362113, 2.50749046482452, 
                     12.5166502527345, -86.6941792229212, 1.66923366612312, 7.93904545083449, 
                     -126.082128066627, -0.0754213178285816, 0.618360785191319, 1538.6589391227, 
                     -0.31799222497181, -1.26341570245756, 2.43556703274056, 18.6459686649981, 
                     -1.52020047556644, -1.33047266147014, -2357.02875362113, 1538.6589391227, 
                     7450231.40845542, -1467.47518565875, -8484.65408805661, 38482.6779430954, 
                     36677.1639056923, 5037.57654299884, -16285.3259202736, 2.50749046482452, 
                     -0.31799222497181, -1467.47518565875, 1.63035346006115, 1.13929112479727, 
                     -17.7642468677628, 2.98912100156106, 7.74881374361886, -54.902708360705, 
                     12.5166502527345, -1.26341570245756, -8484.65408805661, 1.13929112479727, 
                     106.706180072993, -244.571429668905, 24.7937851688084, -114.000942411739, 
                     -690.553652091766, -86.6941792229212, 2.43556703274056, 38482.6779430954, 
                     -17.7642468677628, -244.571429668905, 1534.27886445762, 549.368288309068, 
                     64.0901199201054, 980.291164452268, 1.66923366612312, 18.6459686649981, 
                     36677.1639056923, 2.98912100156106, 24.7937851688084, 549.368288309068, 
                     1781.88219182038, -107.322815612669, -1733.09607995999, 7.93904545083449, 
                     -1.52020047556644, 5037.57654299884, 7.74881374361886, -114.000942411739, 
                     64.0901199201054, -107.322815612669, 507.111649502059, 661.404359682356, 
                     -126.082128066627, -1.33047266147014, -16285.3259202736, -54.902708360705, 
                     -690.553652091766, 980.291164452268, -1733.09607995999, 661.404359682356, 
                     9442.32984843676), .Dim = c(9L, 9L, 4L))

#N
N<-nrow(Y)
#K
K<-4
#J
J<-9

myjagsfile<-"JAGSmodel.txt" # the above jags model
tV = t(ilrBase(D=K))
mydata <-list(N=N, J=J, K=K, Y=Y, mu.s=mu.s, Sigma.s.inv=Sigma.s.inv, R.Sigma=Sigma.inv, k.Sigma=J, tV=tV)
mypars <-c("Sigma.inv","p","kappa","phi", "s")
myiterations<-25000
myburnin<-12500
mythin<-175
set.seed(5000)
run.mcmc <- jags(data=mydata, inits=NULL, parameters.to.save=mypars, model.file=myjagsfile, n.chains=2, 
             n.iter=myiterations, n.burnin=myburnin, n.thin=mythin, DIC=TRUE)
1个回答

rjags 的作者提到，Wishhart 分布只能在强烈意义上表现为共轭分布时才能使用：

http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/eab372de

不过有一种方法可以解决这个问题（不过，我必须感谢R包的作者bamdit，这就是我第一次发现以下技巧的地方：

model
{
for( i in 1 : n ) {
    m[i,1:2] ~ dmnorm(mu.0[1:2], sigma.inv[1:2, 1:2]) ## m is not the data
    w[i] ~ dgamma(nu.2, nu.2) T(0.1, 3)
    y[i, 1] <- mu[1] + m[i, 1] / sqrt(w[i]) ## y is the data
    y[i, 2] <- mu[2] + m[i, 2] / sqrt(w[i]) ## the scale enters here 
}
# Priors ...
mu[1] ~ dnorm(m.0[1], pre.mu[1])
mu[2] ~ dnorm(m.0[2], pre.mu[2])
mu.0[1] <- 0
mu.0[2] <- 0
# Weights distribution
nu.2 <- nu / 2
nu ~ dexp(nu.0) # prior for df
sigma.inv[1:2,1:2] ~ dwish(R[1:2,1:2], k)
sigma[1:2, 1:2] <- inverse(sigma.inv[1:2, 1:2])
}
其它你可能感兴趣的问题
上一篇什么是高斯过程中的平滑 下一篇使用对数转换数据进行线性回归的 R 中的发生率比 (IRR)？