# Estimate MPT parameters as random Rasch subject ability (B) and fixed item difficulty (D)

# datamatrixNK - item level data
# N - number of subjects
# K - number of items
# D[kk,pp] - item difficulty values
# WordErrT[kk] - item neighborhood density
# WordErrL - language neighborhood density

model{

    # Priors on subjects
    for (nn in 1:N){
        for(pp in 1:6){
            B[nn,pp] ~ dnorm(0,1)
        }
    }

    for (nn in 1:N){                       
        for(kk in 1:K){
            for(pp in 1:5){
                param[nn,kk,pp] <- exp(B[nn,pp]-D[kk,pp]) / (1+exp(B[nn,pp]-D[kk,pp]))
            }
                param[nn,kk,6] <- exp(B[nn,6]) / (1+exp(B[nn,6]))
                param[nn,kk,7] <- WordErrT[kk]
                param[nn,kk,8] <- WordErrL

            pk[nn,kk,1] <- param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*param[nn,kk,3]*param[nn,kk,4]*param[nn,kk,5]
            pk[nn,kk,2] <- param[nn,kk,1]*param[nn,kk,6]*(1-param[nn,kk,2])*param[nn,kk,5]
            pk[nn,kk,3] <- (param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*param[nn,kk,3]*param[nn,kk,4]*(1-param[nn,kk,5])*param[nn,kk,7])+(param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*param[nn,kk,3]*(1-param[nn,kk,4])*(1-param[nn,kk,5])*param[nn,kk,7])+(param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*(1-param[nn,kk,3])*param[nn,kk,5])+(param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*(1-param[nn,kk,3])*(1-param[nn,kk,5])*param[nn,kk,7])
            pk[nn,kk,4] <- param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*param[nn,kk,3]*(1-param[nn,kk,4])*param[nn,kk,5]
            pk[nn,kk,5] <- (param[nn,kk,1]*param[nn,kk,6]*(1-param[nn,kk,2])*(1-param[nn,kk,5])*param[nn,kk,8])+(param[nn,kk,1]*(1-param[nn,kk,6])*param[nn,kk,5])+(param[nn,kk,1]*(1-param[nn,kk,6])*(1-param[nn,kk,5])*param[nn,kk,8])
            pk[nn,kk,6] <- (param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*param[nn,kk,3]*param[nn,kk,4]*(1-param[nn,kk,5])*(1-param[nn,kk,7]))+(param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*param[nn,kk,3]*(1-param[nn,kk,4])*(1-param[nn,kk,5])*(1-param[nn,kk,7]))+(param[nn,kk,1]*param[nn,kk,6]*param[nn,kk,2]*(1-param[nn,kk,3])*(1-param[nn,kk,5])*(1-param[nn,kk,7]))
            pk[nn,kk,7] <- (param[nn,kk,1]*param[nn,kk,6]*(1-param[nn,kk,2])*(1-param[nn,kk,5])*(1-param[nn,kk,8]))+(param[nn,kk,1]*(1-param[nn,kk,6])*(1-param[nn,kk,5])*(1-param[nn,kk,8]))
            pk[nn,kk,8] <- 1-param[nn,kk,1]

            # Observed Counts
            datamatrixNK[nn,kk,1:8] ~ dmulti(pk[nn,kk,1:8],1)
            datamatrixNK_post[nn,kk,1:8] ~ dmulti(pk[nn,kk,1:8],1)
        }
    }

}