4

私は次の潜在変数モデルを持っています: 人 j は 2 つの潜在変数 X j1と X j2を持っています。観察できるのは、それらの最大値 Y j = max(X j1 , X j2 ) だけです。潜在変数は二変量正規です。それらはそれぞれ平均 mu、分散 sigma 2を持ち、それらの相関は rho です。n 人の患者 (j = 1,...,n) からのデータを使用して、Y jのみを使用して 3 つのパラメーター (mu、sigma 2、rho) を推定したいと考えています。

このモデルを JAGS に適合させようとしました (そのため、パラメーターに優先順位を付けています) が、コードをコンパイルすることができません。JAGS を呼び出すために使用している R コードを次に示します。最初に、パラメーターの真の値を指定して、データ (潜在変数と観測変数の両方) を生成します。

# true parameter values
mu <- 3
sigma2 <- 2
rho <- 0.7

# generate data
n <- 100
Sigma <- sigma2 * matrix(c(1, rho, rho, 1), ncol=2)
X <- MASS::mvrnorm(n, c(mu,mu), Sigma) # n-by-2 matrix
Y <- apply(X, 1, max)

次に、JAGS モデルを定義し、JAGS サンプラーを実行してサンプルを返す小さな関数を作成します。

# JAGS model code
model.text <- '
model {
  for (i in 1:n) {
    Y[i] <- max(X[i,1], X[i,2]) # Ack!
    X[i,1:2] ~ dmnorm(X_mean, X_prec)
  }

  # mean vector and precision matrix for X[i,1:2]
  X_mean <- c(mu, mu)
  X_prec[1,1] <- 1 / (sigma2*(1-rho^2))
  X_prec[2,1] <- -rho / (sigma2*(1-rho^2))
  X_prec[1,2] <- X_prec[2,1]
  X_prec[2,2] <- X_prec[1,1]

  mu ~ dnorm(0, 1)
  sigma2 <- 1 / tau
  tau ~ dgamma(2, 1)
  rho ~ dbeta(2, 2)
}
'

# run JAGS code. If latent=FALSE, remove the line defining Y[i] from the JAGS model
fit.jags <- function(latent=TRUE, data, n.adapt=1000, n.burnin, n.samp) {
  require(rjags)
  if (!latent)
    model.text <- sub('\n *Y.*?\n', '\n', model.text)
  textCon <- textConnection(model.text)
  fit <- jags.model(textCon, data, n.adapt=n.adapt)
  close(textCon)
  update(fit, n.iter=n.burnin)
  coda.samples(fit, variable.names=c("mu","sigma2","rho"), n.iter=n.samp)[[1]]
}

最後に、JAGS を呼び出して、観測されたデータのみを供給します。

samp1 <- fit.jags(latent=TRUE, data=list(n=n, Y=Y), n.burnin=1000, n.samp=2000)

悲しいことに、これによりエラー メッセージが表示されます。「Y[1] は論理ノードであり、監視できません」。JAGS は、私が "<-" を使用して Y[i] に値を代入することを好みません (問題のある行を "Ack!" で示します)。苦情は理解できますが、これを修正するためにモデル コードを書き直す方法がわかりません。

また、他のすべて (「Ack!」行以外) が正常であることを示すために、モデルを再度実行しますが、今回は、実際に観測されているふりをして X データをフィードします。これは完全に実行され、パラメーターの適切な推定値が得られます。

samp2 <- fit.jags(latent=FALSE, data=list(n=n, X=X), n.burnin=1000, n.samp=2000)
colMeans(samp2)

このモデルを JAGS の代わりに STAN でプログラムする方法を見つけることができれば、それで問題ありません。

4

2 に答える 2

0

確率を同じ重みを持つ 2 つの成分の混合物として扱うスタン言語でこれをモデル化できると思います。スタンコードは次のようになります

data {
  int<lower=1> N;
  vector[N] Y;
}
parameters {
  vector<upper=0>[2] diff[N];
  real mu;
  real<lower=0> sigma;
  real<lower=-1,upper=1> rho;
}
model {
  vector[2] case_1[N];
  vector[2] case_2[N];
  vector[2] mu_vec;
  matrix[2,2] Sigma;
  for (n in 1:N) {
    case_1[n][1] = Y[n]; case_1[n][2] = Y[n] + diff[n][1];
    case_2[n][2] = Y[n]; case_2[n][1] = Y[n] + diff[n][2];
  }
  mu_vec[1] = mu; mu_vec[2] = mu;
  Sigma[1,1] = square(sigma);
  Sigma[2,2] = Sigma[1,1];
  Sigma[1,2] = Sigma[1,1] * rho;
  Sigma[2,1] = Sigma[1,2];
  // log-likelihood
  target += log_mix(0.5, multi_normal_lpdf(case_1 | mu_vec, Sigma), 
                         multi_normal_lpdf(case_2 | mu_vec, Sigma));
  // insert priors on mu, sigma, and rho
}
于 2016-11-12T05:04:03.487 に答える