観測データ y1 と y2 を考えます。y1 は連続スケールで測定され、y2 はバイナリ スケールで測定されます。連続潜在変数 z は、y2 = I(z > 0) として y2 を生成すると仮定されます。(z が正規の場合、y2 はわずかに 2 進プロビットです)。さらに、y1 と z の間の依存関係をモデル化するためにコピュラが使用されます。このモデルは、次のように階層的に記述できます (表記の乱用があります)。
y2 = I(z > 0)
(y1, z) ~ C(F_y1( |w), F_z( |w) | ファイ)
w, ファイ ~ 事前確率
ここで、w は y1 と z の周辺パラメーターのベクトル、F_y1 と F_z はそれぞれ y1 と z の周辺累積分布関数、phi はコピュラ パラメーターです。
これをスタンでどのようにモデル化できますか? コピュラによって生成された二変量尤度から y1 と z をサンプリングするカスタム確率関数を作成しました。どうすればよいかわからないのは、潜在変数 z を説明する (生成する) ことと、y2 と z の関係を指定する方法です。
stan でデータ拡張を使用したプロビット回帰を既に見ましたが、モデルにコピュラがあるため、これは役に立たないようです。
編集：上記のリンクが役に立たないことについて誤解している可能性があります。私は次のコードを書きました。(理論的に)正しいように見えるかどうかについてコメントをいただければ幸いです。
functions {
real copulapdf_log(real[] y1, real[] z, vector mu1, vector mu2, real sigma1, real phi, int n){
real logl;
real s;
logl <- 0.0;
for (i in 1:n){
s <- log(dCphi_du1du2_s(normal_cdf(y1[i],mu1[i],sigma1), logistic_cdf(z[i],mu2[i],1), phi)) + normal_log(y1[i],mu1[i],sigma1) + logistic_log(z[i],mu2[i],1);
logl <- logl + s;
}
return logl;
}
}
data {
int<lower=0> n; // number of subjects
int<lower=0> k1; // number of predictors for y1
int<lower=0> k2; // number of predictors for y2
real y1[n]; // continuous data
real y2[n]; // 0/1 binary data
matrix[n, k1] x1; // predictor variables for y1
matrix[n, k2] x2; // predictor variables for y2
}
transformed data{
int<lower=-1, upper=1> sign[n];
for (i in 1:n) {
if (y2[i]==1)
sign[i] <- 1;
else
sign[i] <- -1;
}
}
parameters {
real phi; // frank copula param
vector[k1] b1; // beta coefficients for y1
vector[k2] b2; // beta coefficients for y2
real<lower=0> abs_z[n]; // abs value of latent variable
real<lower=0> sigma1; // sd for y1's normal distribution
}
transformed parameters {
real v[n];
vector[n] mu1; // location for y1
vector[n] mu2; // location for z
for (i in 1:n) {
v[i] <- sign[i] * abs_z[i];
}
mu1 <- x1 * b1;
mu2 <- x2 * b2;
}
model {
b1 ~ normal(0, 100);
b2 ~ normal(0, 100);
phi ~ normal(0, 10);
increment_log_prob(copulapdf_log(y1, v, mu1, mu2, sigma1, phi, n));
}