ダブルスだけの結果から個人のレーティングとペアの相性を推定する(stan)

llibrary(rstan)

scode="
data{
  int N;
  int M;
  int x[N,2];
  int y[N,2];
  int one[N];
}

parameters {
  real rate[M];
  real rate_two[M,M];
}

transformed parameters {
    real p[N];
    real tmp1;
    real tmp2;
    real b=0.00575364;
    for(i in 1:N){
      tmp1=rate[x[i,1]]+rate[x[i,2]]+rate_two[x[i,1],x[i,2]];
      tmp2=rate[y[i,1]]+rate[y[i,2]]+rate_two[y[i,1],y[i,2]];
      p[i]=inv_logit(b*(tmp1-tmp2));
      print(p[i])
    }
}

model{

  rate~normal(1500,800); #1500がデフォルト

  for(i in 1:M){
    for(j in 1:M){
      if(i<j){
        rate_two[i,j]~normal(0,500);
      }else{
        rate_two[i,j]~normal(0,500);
      }
    }
  }

  for(i in 1:N){
    one[i]~bernoulli(p[i]);
  }

}

"

#b=Solve[1/(1 + Exp[-x*100]) == 0.64, x]

d=c(1,2,3,4,
    1,4,2,3,
    2,3,4,5,
    2,3,1,5,
    1,3,2,5,
    2,3,1,4,
    1,5,2,4
    )

d2=c(1,2,3,4,
    1,4,2,3,
    2,3,4,5,
    2,3,1,5,
    1,3,2,5,
    1,2,4,5,
    1,3,4,5
)

#12が強いが14だと弱くなるペア
d3=c(1,2,3,4,
     1,2,3,4,
     1,3,2,4,
     1,3,2,5,
     1,2,3,5,
     1,2,3,4,
     2,3,4,5,
     2,3,1,4
     )
#あえてシングル
dsingle=c(1,1,2,2,
          1,1,3,3,
          2,2,3,3,
          2,2,4,4,
          3,3,4,4,
          1,1,2,2,
          1,1,3,3,
          2,2,3,3,
          2,2,4,4,
          3,3,4,4
          )

df=matrix(d3,ncol=4,byrow = T)

x=df[,1:2]
y=df[,3:4]
M=max(df)

lst=list(N=nrow(x),x=x,y=y,M=max(df),one=rep(1,nrow(x)))
fit=stan(model_code = scode,data=lst,
         #init = function(){list(rate_two=rtinit,rate=rep(1500,M))},
         chain=1,iter=2000,warmup =500)

la=extract(fit)

rate=apply(la$rate,2,mean)
M=max(df)
rate_two=matrix(0,M,M)
for(i in 1:M){
  for(j in 1:M){
    rate_two[i,j]=mean(la$rate_two[,i,j])
  }
}

#rate_two
         [,1]        [,2]       [,3]       [,4]       [,5]
[1,] 12.860810 182.5917557 161.011254 -118.22797   4.415395
[2,] 14.088220   1.7288561 139.543518  -72.28733 -87.325399
[3,] -3.235520  -9.7427856  -3.568055  -98.99620 -99.269477
[4,] 14.150042  -0.5098703  18.887168   16.07496 -18.809176
[5,] -3.794969  12.4902216   1.515641   11.94288  -7.983007

#rate
[1] 2059.3576 1887.4912 1863.8897  693.6636  973.8996