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

TokyoR#53で発表しました。

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

Tokyo r53 from Hiroshi Shimizu

スライドは画像で張り付けちゃったので，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;
    }
  }
}

data {

int K; // number of rank

int N; // number of data points

int M; // number of items

int y[N,M]; // data

}

parameters{

ordered[K] logit_p[M];

real lambda;

real sigma;

real eta;

vector[K] q;

real alpha;

}

transformed parameters{

simplex[K] pi;

vector[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点とした場合の期待得点を計算しています。

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