Hierarchical Modeling 2 differences types of data using change point model

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]

What values are infinity? Might you be dividing by zero somewhere?

The best thing to do is trace down where they come from, because it’s either overflow or an arithmetic error somewhere.