某研究会で発表したモデルのコード

TokyoR#53で発表しました。

発表スライドは以下です。最後まで見てもらえばわかりますが,実際にこの方法で合否決定はしてません。ネタスライドとして楽しんでもらえればと思います。

 

スライドは画像で張り付けちゃったので,stanコードを置いておきます。簡単な解説も下に書いてます。

data {
  int<lower=1> K; // number of rank
  int<lower=1> N; // number of data points
  int<lower=1> M; // number of items
  int y[N,M]; // data
}

parameters{
  ordered[K] logit_p[M];
  real<lower=0> lambda;
  real<lower=0> sigma;
  real<lower=0> eta;
  vector<lower=0, upper=1>[K] q;
  real<lower=0> alpha;
}

transformed parameters{
  simplex[K] pi;
  vector<lower=0>[K] q2;
  matrix[K,K] C;
  real log_lik[N];
  matrix[N,K] subrank;
  
  for(i in 1:K){
    for(j in 1:K){
      C[i,j] <- eta^2*exp(-(i-j)^2/((K-1)*lambda)^2)+if_else(i==j,sigma,0);
    }
  }
  
  q2[1] <- q[1];
  for (k in 2:K){
    q2[k] <- q[k]*(1 - q[k-1])*q2[k-1]/q[k-1];
  }
  for (k in 1:K){
    pi[k] <- q2[k]/sum(q2);
  }

  {
    real ps[K];
    real temp;
    for (n in 1:N) {
      temp <- 0;
      for (k in 1:K) {
        ps[k] <- log(pi[k]);
        for(m in 1:M){
            ps[k] <- ps[k] + bernoulli_logit_log(y[n,m],logit_p[m,k]);
        }
        temp <- temp + exp(ps[k])*pi[k];
      }
      log_lik[n] <- (log_sum_exp(ps));
      for(k in 1:K){
        subrank[n,k] <- (exp(ps[k])*pi[k])/temp;
      }
    }
  }
}

model{
  for(m in 1:M) logit_p[m] ~ multi_normal(rep_vector(0,K),C);
  q ~ beta(1, alpha);
  eta ~ cauchy(0,2.5);
  lambda ~ cauchy(0,2.5);
  sigma ~ cauchy(0,2.5);
  for(n in 1:N) increment_log_prob(log_lik[n]);
}

generated quantities{
  vector[K] p[M];
  vector[K] score;
  for(k in 1:K){
    score[k] <-0;
    for(m in 1:M){
      p[m,k] <- inv_logit(logit_p[m,k]);
      score[k] <- score[k] + p[m,k]*5;
    }
  }
}

 

データとして入れるのはランク数K,サンプルサイズN,変数の数N,そしてデータセットです。

推定するパラメータはロジット変換した正答率,ガウス過程のハイパーパラメータであるη,λ,σ,そしてqとαはディリクレ過程で使うパラメータで,これらはtransformed parametersで混合率πの計算に使います。

transformed parametersでは上から順番に,ガウス過程のカーネル関数の計算,ディリクレ過程のπの計算,最後に混合ベルヌーイ分布の尤度の計算です。今回は 2値データなのでベルヌーイ分布でしたが,任意の分布に変えることができます。もちろん,それに合わせてパラメータも変わりますが。

modelでは事前分布を指定して,あとはincrement_log_probに尤度を入れて終わりです。

generated quantitiesでは各ランクで項目ごとの正答率と,1問配点5点とした場合の期待得点を計算しています。

もし「ここ変だよ!」とかあったらご指摘ください。

This entry was posted in 心理統計学. Bookmark the permalink.