rStan beta

# fake data
library(bayesAB)
library(MCMCpack)

plotInvGamma(2, 2.5) + geom_vline(xintercept = 1)

set.seed(8675309)
n <- 60
p <- 500

person <- sort(rep(1:p,n))
item <- rep(1:n,p)
itemMu <- rnorm(n,0,2)
personMu <- rnorm(p,0,2)
itemDisc <- rinvgamma(n,2,2.5)

df <- data.frame(person,item)
df$discrimination <- itemDisc[df$item]
df$difficulty <- itemMu[df$item]
df$ability <- personMu[df$person]

df$scoreexp <- exp(df$discrimination * (df$ability - df$difficulty) + rnorm(nrow(df),0,1)) 
df$score <- 1/(1+df$scoreexp)

dfStan <- df[,c("person", "item", "score")]


# model

twoplcovstan_yo <- '
data {
  int<lower=1> I;               // # word
  int<lower=1> J;               // # persons
  int<lower=1> N;               // # observations
  int<lower=1, upper=I> ii[N];  // word for n
  int<lower=1, upper=J> jj[N];  // person for n
  real<lower=0, upper=1> y[N];  // brier
}
parameters {
  vector<lower = 0>[I] alpha;     // discrimination for item i
  vector[I] beta;                 // difficulty for item i
  vector[J] theta;              // ability for person j
}
model {
  vector[N] eta;
  alpha ~ normal(1,3);
  beta ~ normal(0,2);
  theta ~ normal(0,2);

  for (n in 1:N) {
    eta[n] = inv_logit(alpha[ii[n]] * (theta[jj[n]] - beta[ii[n]]));
    print(eta[n]);
  }
  
  y ~ beta_proportion(eta,6);
}

'
standata_yo <- list(I = length(unique(dfStan$item)),
                    J = length(unique(dfStan$person)),
                    N = nrow(dfStan),
                    ii = dfStan$item, 
                    jj = dfStan$person,
                    y = dfStan$score)

twopl.fit_yo <- stan(model_code = twoplcovstan_yo,
                     data = standata_yo,
                     iter = 4000,
                     chains = 4,
                     seed = 4204)

twopl.summary_yo <- summary(twopl.fit_yo)$summary