data {
  int<lower=1> D;
  int<lower=0> N_obs;
  int<lower=0> N_cens;
  int<lower=0> M;
  real<lower=0> sigma_tau;
  int<lower=1, upper=M> jj_obs[N_obs]; // judge_ID
  int<lower=1, upper=M> jj_cens[N_cens]; // judge_ID
  int<lower=0, upper=1> dec_obs[N_obs];
  int<lower=0, upper=1> dec_cens[N_cens];
  row_vector[D] X_obs[N_obs];
  row_vector[D] X_cens[N_cens];
  int<lower=0, upper=1> y_obs[N_obs];
}

parameters {
  vector[N_obs] Z_obs;
  vector[N_cens] Z_cens;
  real alpha_T[M];
  vector[D] beta_XT_raw;
  vector[D] beta_XY_raw;
  real<lower=0> beta_ZT_raw; // jungimainen oletus latentin pos. vaik.
  real<lower=0> beta_ZY_raw; // jungimainen oletus latentin pos. vaik.
  
  real<lower=0> tau_XT;
  real<lower=0> tau_XY;
  real<lower=0> tau_ZT;
  real<lower=0> tau_ZY;

}

transformed parameters {
  vector[D] beta_XT;
  vector[D] beta_XY;
  real<lower=0> beta_ZT; // jungimainen oletus latentin pos. vaik.
  real<lower=0> beta_ZY; // jungimainen oletus latentin pos. vaik.

  beta_XT = tau_XT * beta_XT_raw;
  beta_XY = tau_XY * beta_XY_raw;
  beta_ZT = tau_ZT * beta_ZT_raw;
  beta_ZY = tau_ZY * beta_ZY_raw;
}

model {
  Z_obs ~ normal(0, 1);
  Z_cens ~ normal(0, 1);
  
  tau_XT ~ normal(0, sigma_tau);
  tau_XY ~ normal(0, sigma_tau);
  tau_ZT ~ normal(0, sigma_tau);
  tau_ZY ~ normal(0, sigma_tau);
  
  beta_XT_raw ~ normal(0, 1);
  beta_XY_raw ~ normal(0, 1);
  beta_ZT_raw ~ normal(0, 1);
  beta_ZY_raw ~ normal(0, 1);
  
  for(i in 1:N_obs){
    dec_obs[i] ~ bernoulli_logit(alpha_T[jj_obs[i]] + X_obs[i] * beta_XT + beta_ZT * Z_obs[i]);
    y_obs[i] ~ bernoulli_logit(X_obs[i] * beta_XY + beta_ZY * Z_obs[i]);
  }
  
  for(i in 1:N_cens)
    dec_cens[i] ~ bernoulli_logit(alpha_T[jj_cens[i]] + X_cens[i] * beta_XT + beta_ZT * Z_cens[i]);
}

generated quantities {
  int<lower=0, upper=1> y_est[N_cens];
  
  for(i in 1:N_cens){
    y_est[i] = bernoulli_logit_rng(X_cens[i] * beta_XY + beta_ZY * Z_cens[i]);
  }
}