Hi users,
I try to execute my model but I find a lot of value Inf. So then I would like to know if my implemented model is correctly written.
In fact I have 2 types datas. One includes discrete values (0,1) and other continues values.
In each line of each data, it can have either one change point or not.
Firstly, I want to detect for each line, these change points.
Secondly, the first data can be modelized by Rasch model (includes one parameter for an examinee - P1 and one parameter for an item - I1 ) and the second by lognormal (includes one a parameter for examinee - P2 and one parameter for an item - I2). However, there is a relation between parameters of these models.
From these hypotheses, I implemented the algorithms below :
Thanks for your helping.
[model <- 'data {
int J; // number of items
int I; // number of examinees
int items[I,J]; // responses
real times[I,J]; // log responses-times
real<lower = 0> eta; // shape parameter for LKJ ditribution
}
transformed data {
real log_unif;
log_unif = -log(J);
}
parameters {
vector[2] paramP[I]; // parameters of examinees
vector<lower=0>[2] paramPnew[I]; // new parameters of examinees for decreasing abilities of examinees and increasing speeds after the change point
vector[2] paramI[J]; // parameters of items
real<lower=0> lambda[I];
// relation between the parameters of examinees (P1,P2)
vector[2] mu_P;
cholesky_factor_corr[2] L_Omega_P;
vector<lower=0>[2] L_sigma_P;
// relation between the parameters of examinees (I1,I2)
vector[2] mu_I;
cholesky_factor_corr[2] L_Omega_I;
vector<lower=0>[2] L_sigma_I;
}
transformed parameters {
matrix[I,J] lpitems;
matrix[I,J] lptimes;
vector<lower=0>[I] sigma;
for(i in 1:I)
sigma[i] = inv_sqrt(lambda[i]);
//Change point model with response
lpitems = rep_matrix(log_unif, I, J);
for(i in 1:I){
for (s in 1:J){
for (j in 1:J)
lpitems[i,s] = lpitems[i,s] + bernoulli_logit_lpmf(items[i,j] | j < s ? (paramP[i,1] - paramI[j,1]) : (paramP[i,1] - paramI[j,1] - paramPnew[i,1]));
}
}
//Change point model with response-times
lptimes = rep_matrix(log_unif, I, J);
for(i in 1:I){
for (d in 1:J){
for (j in 1:J)
lptimes[i,d] = lptimes[i,d] + normal_lpdf(times[i,j] | j < d ? (paramI[j,2] - paramP[i,2]) : (paramI[j,2] - paramP[i,2] - paramPnew[i,2]), sigma[i]);
}
}
}
model {
// Group Means
mu_P ~ normal(0, 5);
mu_I ~ normal(0, 5);
// Common Precision
lambda ~ cauchy(0, 2.5);
// Which Side is Time of Change Point ? Responses come from Rasch model and Responses-times from Gaussian
target += log_sum_exp(lpitems[I]);
target += log_sum_exp(lptimes[I]);
// Assuming a slight correlation between parameters of examinees
L_Omega_P ~ lkj_corr_cholesky(eta);
L_sigma_P ~ cauchy(0, 2.5);
mu_P ~ normal(0, 1);
paramP ~ multi_normal_cholesky(mu_P, diag_pre_multiply(L_sigma_P, L_Omega_P));
// Assuming a slight correlation between parameters of items
L_Omega_I ~ lkj_corr_cholesky(eta);
L_sigma_I ~ cauchy(0, 2.5);
mu_I ~ normal(0, 1);
paramI ~ multi_normal_cholesky(mu_I, diag_pre_multiply(L_sigma_I, L_Omega_I));
}
'
]
[edit: escaped code]