３．４． ２群のデータ ぱらぱらめくる『はじめての統計データ分析』

• t検定が、等分散性の仮定の有無を区別するように、また、対応のありなしを区別するように、事後分布推定もそれぞれの仮定の下で行う
• 仮定はmodelとして指定する
  • 等分散性の有無は、２群の分散を１つのパラメタで表すか否かで書き分ける
  • 対応があることは、分散共分散行列の係数推定を入れることで表す

data {
  int<lower=0> n1;     int<lower=0> n2;       #データ数
  real x1[n1];         real x2[n2];           #データ
  real EQU;                                   #ＳＤ、1 - 共通、0 - 異なる
  real mL;real mH;real sL;real sH;            #事前分布
}
parameters {
  vector<lower=mL,upper=mH> [2] mu;          #平均(範囲指定)
  real<lower=sL,upper=sH> sigma1;            #標準偏差(範囲指定)
  real<lower=sL,upper=sH> dummy;             #ダミーのＳＤ(範囲指定)
}
transformed parameters {
  real<lower=0>           sigma2;             #標準偏差2
  sigma2<-if_else(EQU>0.5,sigma1,dummy);
}
model {
  x1 ~ normal(mu[1],sigma1);                  #正規分布1
  x2 ~ normal(mu[2],sigma2);                  #正規分布2
}
generated quantities{
  real xaste[2];
  real log_lik;
  xaste[1]<- normal_rng(mu[1],sigma1);        #予測分布1
  xaste[2]<- normal_rng(mu[2],sigma2);        #予測分布2
  log_lik<-normal_log(x1,mu[1],sigma1)+
           normal_log(x2,mu[2],sigma2);       #対数尤度
}

data { 
  int<lower=0> n;                             #データ数   
  vector[2] x[n];                             #データ     
  real EQU;                                   #ＳＤ、1 - 共通、0 - 異なる
  real mL;real mH;real sL;real sH;            #事前分布
}
parameters {
  vector<lower=mL,upper=mH> [2] mu;          #平均(範囲指定)
  real<lower=sL,upper=sH> sigma1;            #標準偏差(範囲指定)
  real<lower=sL,upper=sH> dummy;             #ダミーのＳＤ(範囲指定)
  real<lower=-1,upper=1>  rho;                #相関
}
transformed parameters {
  real<lower=0>           sigma2;             #標準偏差2
  cov_matrix[2]           Sigma;
  sigma2<-if_else(EQU>0.5,sigma1,dummy);
  Sigma[1,1] <- pow(sigma1,2);                #共分散行列
  Sigma[2,2] <- pow(sigma2,2);
  Sigma[1,2] <- sigma1*sigma2*rho;
  Sigma[2,1] <- Sigma[1,2];
}
model {
  for(i in 1:n){x[i]~multi_normal(mu,Sigma);} #2変量正規分布
}
generated quantities{
  vector[2] xaste;
  real log_lik;
  xaste <-  multi_normal_rng(mu,Sigma);       #予測分布
  log_lik <- multi_normal_log(x, mu, Sigma);  #対数尤度
}