使用ぎざぎざまたはSTAN観測ノードが潜在ノードの最大値である

Question

I有し以下潜在変数モデル：人のJは、2つの潜在変数、X _J1とX _J2を有しています。私たちが観察できるのは、最大値Y _j = max（X _j1、X _j2）です。潜在変数は二変量法線である。それらはそれぞれ平均μ、分散σを有し、それらの相関はρである。私は、Y患者、j = 1、...、nからのデータを用いて、Y _jのみを使用して3つのパラメータ（μ、Σ、rho）を推定したいと考えています。

私はこのモデルをJAGSに適合させようとしています（パラメータにpriorsを設定しています）。しかし、コンパイルするコードを取得できません。 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) 

それからぎざぎざモデルを定義し、ぎざぎざサンプラーを実行し、サンプルを返すように少し機能を記述します。

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

は残念ながら、これはエラーメッセージをもたらす：「Y [1]の論理ノードであると観察することができない」

# 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は "< - "を使ってY [i]に値を代入するのを私に好まれません（私は "Ack！"という違反行を表します）。私は苦情を理解していますが、これを修正するためにモデルコードを書き換える方法がわかりません。

また、「Ack！」行以外にも問題はないことを実証するために、モデルを再実行しますが、今回は実際に観測されたXデータをフィードします。これは完璧に動作し、私はパラメータの良い推定値を得る：

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

あなたは私と一緒にいいと思いSTANの代わりに、ぎざぎざにこのモデルをプログラムする方法を、見つけることができれば。

Answer 1

私は、このことをStan言語で、尤度を等しい重みの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 
} 

Answer 2

のようになります。理論的にはあなたが最大とないの総和をモデル化しているとして、この場合にはハックのビットを使用DSUM分布を（使用してぎざぎざで、このようなモデルを実装することができます2つの変数）。

set.seed(2017-02-08) 

# 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) 

model.text <- ' 
model { 
    for (i in 1:n) { 
    Y[i] ~ dsum(max_X[i]) 
    max_X[i] <- max(X[i,1], X[i,2]) 
    X[i,1:2] ~ dmnorm(X_mean, X_prec) 
    ranks[i,1:2] <- rank(X[i,1:2]) 
    chosen[i] <- ranks[i,2] 
    } 

    # 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) 

    #data# n, Y 
    #monitor# mu, sigma2, rho, tau, chosen[1:10] 
    #inits# X 
} 
' 

library('runjags') 

results <- run.jags(model.text) 
results 
plot(results) 

2つのことに注意すること：

1. ぎざぎざが十分にスマートではありません - （後述それはどんな本当の意味での「仕事」をしていませんが）しかし、次のコードはコンパイルして実行しdsum（max（X [i、]））制約を満たしながらXの行列を初期化する必要があります。したがって、わかりやすい値を使ってJAGSをXに初期化する必要があります。この場合、私は不正なシミュレーション値を使用しています。あなたが得る答えは、Xの初期値の選択に大きく依存します。実際の世界では、シミュレートされた値は戻ってこないでしょう。

2. max（）制約は、一般的なフレームワークの中で解を考えることができない問題を引き起こします。一方のパラメータが減少する一方で両方のパラメータがまったく使用されない通常のdsum制約X [i、]のmin（）値は無視されるため、サンプラーは自由に行うことができます。これは、サンプラーが2つのX [i]間で 'スイッチ'するのに必要な条件である、Y [i]と同一であるmin（X [i、]）の値に至ることは非常にまれにしか、]。だから、スイッチングは起こりませんし、初期化時に最大値になるように選択されたX []は、これを説明するトレースパラメータ 'selected'を追加しました。

   限り、私は質問私はあるかもしれないが、（ここでは根本的な問題だと思う基本的に同じ非混合罠に陥るだろう「どのように私はこのコードでください」と他の潜在的な解決策を見ることができるよう

間違っていて、そうでないことを示すBUGS/JAGS/Stanコードを大いに歓迎します）。

モデル選択のためのCarlin & Chibbメソッドに類似したものが動作する可能性があります（切り替えを促すためにmin（pseudo_X）パラメータをYに等しくすることができます）。これはうまくいかないかもしれませんが、合理的な量のバグやJAGSの経験があれば、試してみてください。 Carlin、BP、Chib、S.、1995. Bayesian model choiceマルコフ連鎖モンテカルロ法を介して。 J.R.Stat。 Soc。 Ser。 B 57,473-484。

また、問題を多少異なる方法で考えて、最初の列がすべて不足していて2番目の列がすべてYの行列としてXを直接モデル化することもできます。次に、dinterval（）を使用してそれらが対応する最大値よりも低くなければならない欠損値。私はmu/sigma2/rhoを見積もる上でどれくらいうまくいくのかは分かりませんが、試してみる価値があるかもしれません。

ちなみに、これは必ずしもあなたの質問に答えるわけではないが、私はそれが「コード化可能か」と「それは実行可能か」の違いの有用な例だと思う。

マット

psはるかに賢い解決策は、最大2つの正規分布の分布を直接考慮することです。そのような分布が存在するかどうかはわかりませんが、それがPDFである場合、分布を直接最小値の価値を全く考慮する必要なしに0/1のトリックを使用します。