Non-compensatory IRT model

I played with this model a bit with generated. I think the frac_rt1s were basically all non-identified.

I have a version of the model here where I pass the frac_rt1s in and it’s kinda sketchy, but it’ll get in the ballpark.

I think the problem is that unless we make some assumptions about how the different items are related, there’s no way to figure out how expensive each part of the test is.

Like consider the noise-free situation:

r_1 + r_2 = r_t

We measure r_t and need to solve for r_1 and r_2. If you had alternate tests and they had alternate test times r_1' and r_2', then you could measure:

r_1 + r_2 = r_{ta}
r_1 + r_2' = r_{tb}
r_1' + r_2 = r_{tc}
r_1' + r_2' = r_{td}

And then you could solve for the individual test times (four equations, four unknowns).

So if there are shared components in the different items, maybe we have a chance. Otherwise I’m not sure it’ll work. Anyway, code is here, play around with it. I could definitely be wrong.

library(tidyverse)
library(cmdstanr)

I = 1 # number of items
J = 2 # number of examinees
N = 200 # number of observed data points
ii = sample(c(1:I), N, replace = TRUE) # item id
jj = sample(c(1:J), N, replace = TRUE) # person id

skill_test1 = rnorm(J)
skill_test2 = rnorm(J)

alpha1 = rnorm(I, 2.0, 1.0)
alpha2 = rnorm(I, 2.0, 1.0)

sigma1 = abs(rnorm(I))
sigma2 = abs(rnorm(I))

phi = abs(rnorm(2))

rt1 = rlnorm(N, alpha1[ii] - phi[1] * skill_test1[jj], sigma1[ii])
rt2 = rlnorm(N, alpha2[ii] - phi[2] * skill_test2[jj], sigma2[ii])

rtt = rt1 + rt2

mod = cmdstan_model("non_comp_irt.stan")

fit = mod$sample(list(I = I,
                      J = J,
                      N = N,
                      ii = ii,
                      jj = jj,
                      rt = rtt,
                      frac_rt1 = rt1 / rtt),
                 chains = 2,
                 parallel_chains = 2,
                 iter_warmup = 500, iter_sampling = 500)
data {
  int<lower = 1> I;                      // number of items
  int<lower = 1> J;                      // number of examinees
  int<lower = 1> N;                      // number of observed data points
  int<lower = 1> ii[N];                  // item id
  int<lower = 1> jj[N];                  // person id
  vector[N] rt;                          // response times
  vector<lower=0, upper=1>[N] frac_rt1;
}

parameters {
  vector[I] alpha1;
  vector[I] alpha2;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
  vector<lower = 0.0>[2] phi;
  vector[J] skill_test1;
  vector[J] skill_test2;
}

model {
  vector[N] rt1 = rt .* frac_rt1;
  vector[N] rt2 = rt  - rt1;
  vector[N] mu1 = alpha1[ii] - phi[1] * skill_test1[jj];
  vector[N] mu2 = alpha2[ii] - phi[2] * skill_test2[jj];
  
  skill_test1 ~ normal(0, 1.0);
  skill_test2 ~ normal(0, 1.0);
 
  sigma1 ~ normal(0, 1);
  sigma2 ~ normal(0, 1);
  
  alpha1 ~ normal(2, 1);
  alpha2 ~ normal(2, 1);
  
  phi ~ normal(0, 1);

  rt1 ~ lognormal(mu1, sigma1[ii]);
  rt2 ~ lognormal(mu2, sigma2[ii]);
}
1 Like

The best way to check for these kinds of identifiability is to an estimate a model where everything is observed and then replace observed data with parameters. That’ll pinpoint where it becomes unidentified. Otherwise can be pretty hard to tell!

2 Likes

Dear Ben, Dear Robert:

Thank you very much for your help!

I took your advice @saudiwin and specified the model with only observed data, and replaced the data with parameters successively (I hope I understood you correctly, Robert). In any case, sampling became problematic once frac_rt1 was switched from data to parameter. So this essentially confirms @bbbales2, who came to the same conclusion.

Thank you for this insightful comment, Ben! Yes, it is possible to argue that the components share the time intensity parameter \alpha: rt_{1i} = \alpha_i - \phi_1 + \omega_{1j} and rt_{2i} = \alpha_i - \phi_2 + \omega_{2j}, with i = 1,2,...,I the number of items (tests) and j = 1,2,...,J the number of test takers. So, you would have only one \alpha vector in the model specification:

parameters {
  vector[I] alpha;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
  vector<lower = 0.0>[2] phi;
  vector[J] skill_test1;
  vector[J] skill_test2;
}

model {
  vector[N] rt1 = rt .* frac_rt1;
  vector[N] rt2 = rt  - rt1;
  vector[N] mu1 = alpha[ii] - phi[1] * skill_test1[jj];
  vector[N] mu2 = alpha[ii] - phi[2] * skill_test2[jj];
  
  skill_test1 ~ normal(0, 1.0);
  skill_test2 ~ normal(0, 1.0);
 
  sigma1 ~ normal(0, 1);
  sigma2 ~ normal(0, 1);
  
  alpha1 ~ normal(2, 1);
  
  phi ~ normal(0, 1);

  rt1 ~ lognormal(mu1, sigma1[ii]);
  rt2 ~ lognormal(mu2, sigma2[ii]);
}

So, if we assume the \alpha parameters to be identical across the components, do I need an ODE system in the model specification to make use of this assumption?

Best wishes
Christoph

Is that what you mean?

r_{1ij} \sim logN(\alpha_i - \phi_1 \omega_{1j}, \sigma_1) \\ r_{2ij} \sim logN(\alpha_i - \phi_2 \omega_{2j}, \sigma_2)

This is sorta assuming that the time to take the first and second parts of the test are the same. Try it out with the simulated data and see what happens! That seems like a pretty strong assumption so maybe it’ll make it identifiable.

Yes, that’s what I meant. I set \alpha_{1i} = \alpha_{2i} and ran some simulations, but it did not really improve anything (still E-BMFI < 0.05, low effective sample sizes and convergence issues).

In my naivety I thought about adding an interaction term \chi_{j} = \phi_{3}\omega_{1j}\omega_{2j} to both latent response time components, so that they have something in common (if I am not mistaken, this makes the model partially compensatory, but I am not sure):

rt_{1} \sim LN(\alpha_{1i} - \phi_1\omega_{1j} + \chi_{j}, \sigma_1)
rt_{2} \sim LN(\alpha_{2i} - \phi_2\omega_{2j} + \chi_{j}, \sigma_2)

I ran some more simulations, but unfortunately this doesn’t change anything either (I have still very low effective sample sizes and convergence issues). The funny thing is that it seems to work with I = 1 and J < 5, but fails as soon as I increase the number of items and the number of test takers.

So I am at a loss again. I am wondering whether it is possible to make the frac_rt1 parameters dependent on some other parameters of the model, or if that would only shift the identification problem to another set of parameters?

Current data generation:

I = 1                                  # number of items
J = 2                                  # number of examinees
N = I*J                                # number of observed data points
ii = rep(1:I, times = J)               # item id
jj = rep(1:J, each = I)                # person id

omega1 = rnorm(J)
omega2 = rnorm(J)

alpha1 = 1.55  
alpha2 = 1.15 
sigma1 = 0.45  
sigma2 = 0.25

phi = c(0.35,0.20,0.15)

rt1 = rlnorm(N, alpha1[ii] - phi[1] * omega1[jj] + (phi[3] * (omega1[jj] * omega2[jj])), sigma1[ii])
rt2 = rlnorm(N, alpha2[ii] - phi[2] * omega2[jj] + (phi[3] * (omega1[jj] * omega2[jj])), sigma2[ii])
rtt = rt1 + rt2

rt_data = list(I = I, J = J, N = N, ii = ii, jj = jj, rt = rtt, frac_rt1 = rt1 / rtt)

And the current model:

data {
  int<lower = 1> I;                      // number of items
  int<lower = 1> J;                      // number of examinees
  int<lower = 1> N;                      // number of observed data points
  int<lower = 1> ii[N];                  // item id
  int<lower = 1> jj[N];                  // person id
  vector[N] rt;                          // response times
}

parameters {
  vector[I] alpha1;
  vector[I] alpha2;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
  positive_ordered[3] phi;
  vector[J] omega1;
  vector[J] omega2;
  vector[N] frac_rt1_raw;
}

transformed parameters{
  vector<lower=0,upper=1>[N] frac_rt1 = inv_logit(frac_rt1_raw);
}

model {
  vector[N] rt1 = rt .* frac_rt1;
  vector[N] rt2 = rt  - rt1;

  omega1 ~ normal(0, 1);
  omega2 ~ normal(0, 1);
  
  sigma1 ~ normal(0, 1);
  sigma2 ~ normal(0, 1);
  
  alpha1 ~ normal(2, 1);
  alpha2 ~ normal(2, 1);
  
  phi ~ normal(0, 1);
  frac_rt1_raw ~ logistic(0, 1);
  
  rt1 ~ lognormal(alpha1[ii] - phi[3] * omega1[jj] + (phi[1] * (omega1[jj] .* omega2[jj])), sigma1[ii]);
  rt2 ~ lognormal(alpha2[ii] - phi[2] * omega2[jj] + (phi[1] * (omega1[jj] .* omega2[jj])), sigma2[ii]);
}

Best wishes
Christoph

If the explicit frac_rt1 formulation is not working then I recommend you try the original convolve1_kernel model again. CmdStan 2.26 fixed an issue in integrate_1d so maybe the error about exceeding relative tolerance has disappeared.

3 Likes

Hi Niko, Hi All:

Thanks for your reply!

I tried to run the model with the convolve1_kernel function (involving integrate_1d), but I still get various error messages that prevent the sampler from even starting.

The first tells me that it is rejecting the initial value,

Chain 2  Rejecting initial value:
Chain 2  Log probability evaluates to log(0), i.e. negative infinity.
Chain 2  Stan can't start sampling from this initial value.

and the second one is the message about the relative tolerance again

Exception: Exception: integrate: error estimate of integral 5.39116e-09 exceeds the given relative tolerance times norm of integral.

The current model with the convolve1_kernel function looks like this:

functions{ 
  real kernel(real x, real xc, real[] theta, data real[] x_r, data int[] x_i) {
    real out;
    real mu1    = theta[1];     // mean of (log)normal1
    real sigma1 = theta[2];     // sd of (log)normal1
    real mu2    = theta[3];     // mean of (log)normal2
    real sigma2 = theta[4];     // sd of (log)normal2

    // convolution [(log)normal1 + (log)normal2 distribution]
    out = exp( lognormal_lpdf( x | mu1, sigma1) + lognormal_lpdf( x_r[1] - x | mu2, sigma2) );
    return out;
  }
  
  real convolve1_kernel(real a, real b, real[] theta, data real x) {
    int x_i[0];
    return(integrate_1d(kernel, a, b, theta, rep_array(x, 1), x_i, 1E-6));
  }
}

data{
  int<lower = 1> I;                      // number of items
  int<lower = 1> J;                      // number of examinees
  int<lower = 1> N;                      // number of observed data points
  int<lower = 1> ii[N];                  // item id
  int<lower = 1> jj[N];                  // person id
  real<lower=0> rt[N];                   // response times
}

parameters{
  vector[J] omega1;
  vector[J] omega2;
  vector<lower=0>[I] alpha1;
  vector<lower=0>[I] alpha2;
  positive_ordered[2] phi;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
}

transformed parameters{
  vector[N] mu1 = alpha1[ii] - phi[2] * omega1[jj];
  vector[N] mu2 = alpha2[ii] - phi[1] * omega2[jj];
}

model{
  omega1 ~ std_normal();
  omega2 ~ std_normal();
  
  alpha1 ~ normal(2, 1);
  alpha2 ~ normal(2, 1);
  
  phi ~ normal(0, 1);
  
  sigma1 ~ normal(0, 1);
  sigma2 ~ normal(0, 1);
  
  for(n in 1:N) target += log( convolve1_kernel(0.0, rt[n], {mu1[n], sigma1[ii[n]], mu2[n], sigma2[ii[n]]}, rt[n]) );
}

If I try the other convolve_kernel function (with integrate_ode_rk45), the error says

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: Exception: Exception: lognormal_lpdf: Random variable is -8.88178e-16, but must be nonnegative!. 

In this case, the functions-block looks like this:

functions{ 
  
  real[] kernel(real x, real[] y, real[] theta, data real[] x_r, data int[] x_i) {
    real dydx[1];
    real mu1    = theta[1];     // mean of lognormal1
    real sigma1 = theta[2];     // sd of lognormal1
    real mu2    = theta[3];     // mean of lognormal2
    real sigma2 = theta[4];     // sd of lognormal2
    
    // convolve both lognormal distributions
    dydx[1] = exp( lognormal_lpdf( x | mu1, sigma1) + lognormal_lpdf( x_r[1] - x | mu2, sigma2) );
    return dydx;
  }
  
  real convolve_kernel(real a, real b, real[] theta, data real x, real rel_tol, real abs_tol) {
    int x_i[0];
    return(integrate_ode_rk45(kernel, rep_array(0.,1), a, rep_array(b,1), theta, rep_array(x,1), x_i)[1,1]);
  }
}

Not that it would change much, but I tried to upgrade this block to the new ODE interface (with the help of the case study by Ben and Sebastian), but I do not understand how to convert the convolve_kernel function. If I am not mistaken, the kernel function above in the new interface would look like this:

  vector kernel(real t, vector y, real mu1, real sigma1, real mu2, real sigma2, real x) {
    vector[1] dydx;
    // convolve both lognormal distributions
    dydx[1] = exp( lognormal_lpdf( t | mu1, sigma1) + lognormal_lpdf( x - t | mu2, sigma2) );
    return dydx;
  }

but how would the ode_rk45 solver interface look like in this case? When I try to set it up analogous to the old interface, it always tells me that kernel is not found… as mentioned before, I am aware this would change nothing in the outcome of my unsuccessful attempts to get the convolution model running, but I’d like to learn.

If you have any questions, or need additional information, please let me know!

What you can try doing is setting one value for frac_rt1 to a fixed value (say 1 or 0) and see if that allows the model to sample the parameters relative to that constant. For simulated data, just use the “true” value and then see if the other free parameters will take on the true values as well.

If it is unidentified, this doesn’t really change the model as adding the constant means you’re adding information to give it a single mode.

Exception: Exception: integrate: error estimate of integral 5.39116e-09 exceeds the given relative tolerance times norm of integral .

That’s a bug in Stan that got fixed in 2.26: Release of CmdStan 2.26.0 . Probly the easiest way to use that is cmdstanr

I did. The last results and error messages mentioned in my last post were in simulations done in cmdstan 2.26.0 via cmdstanr.

Another thing I found today: integrate_ode_rk45 works with one item and two test takers, while integrate_1d does not. Once I run the model with more items and more test takers, both do not work anymore. I am beginning to think the error is somewhere else, but I have no idea where…

Edit: That’s the output when I load the cmdstanr-package

> library(cmdstanr)
This is cmdstanr version 0.3.0.9000
- Online documentation and vignettes at mc-stan.org/cmdstanr
- CmdStan path set to: C:/Users/ckoenig/Documents/.cmdstanr/cmdstan-2.26.0
- Use set_cmdstan_path() to change the path

Sorry! Behind on some other stuff, if I don’t respond sometime tomorrow ping me again.

No worries!

In the meantime, I did some further digging, and I am wondering if the kernel in the integrate_1d version of the model is set up correctly. I added some print statements to see what the parameters of the two lognormal distributions, the value of the integral, and the response time are. Here’s what I found for the integrate_1d version of the model (it’s an excerpt of the print-output):

                  mu1 ,  sigma1 ,      mu2 ,   sigma2 ,        value,   rt[n]
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 1.67579e-06 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 4.80961e-13 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 4.2959e-23 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 6.03236e-11 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 9.92889e-11 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 1.09195e-11 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0.00719056 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 2.77921e-12 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 5.14692e-122 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 8.40563e-10 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 1.56489e-09 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0.00136742 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 8.7948e-13 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 1.09369e-07 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 8.51028e-13 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 6.20755e-57 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 1.22079e-11 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 8.6849e-240 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 2.77081e-10 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 8.83661e-10 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 0 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 4.35777e-08 , 9.04837 
    Chain 1 -0.308038 , 4.40524 , 0.623202 , 0.216137 , 6.26381e-05 , 9.04837 

As can be seen, it 0 all the time, and the integral (not shown) is very small. Moreover, it does not match the integral of the R-function, or the exposed stan-function.

As a comparison, the output for the integrate_ode_rk45 version:

   mu1 ,   sigma1 ,     mu2 ,   sigma2 , value        , rt[n]
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [4.57871e-114] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.20465e-94] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [2.91694e-55] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.1389e-51] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [8.40904e-48] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [8.40904e-48] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.91545e-29] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [2.2267e-24] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.08777e-11] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.8199e-10] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.65875e-09] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.65875e-09] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.8977e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00057898] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.125052] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.153581] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.166528] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.166528] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [4.21655e-08] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.25028e-07] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [9.98243e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.87166e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.8977e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.8977e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.11154e-08] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [8.1785e-08] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [4.45558e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [8.0317e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.60671e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.60671e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [5.45134e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [9.49629e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000956631] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00134517] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00200791] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00200791] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00386584] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00520896] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0179617] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0215265] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0266011] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0266011] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00285658] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00337873] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00723843] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00818435] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00949476] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00949476] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0127254] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0146093] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0269781] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0297124] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.033358] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.033358] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0418711] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0465173] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.072892] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0779737] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0844121] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0844121] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0398428] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0433225] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0627921] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0665531] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0713496] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0713496] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0807985] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0855973] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.109534] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.113656] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.118699] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.118699] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.128151] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.132606] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.151361] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.153982] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.156924] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.156924] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.161704] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.163485] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.166306] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.165784] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.164728] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.164728] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.161289] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.162982] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.166522] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.166314] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.165717] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.165717] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.163739] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.162339] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.151784] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.149378] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.146186] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.146186] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.138538] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.134452] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.112487] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.108463] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.103439] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.103439] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0941432] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0895825] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0682218] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0647387] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0605371] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0605371] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0534174] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0500765] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0355719] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0333679] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0307647] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0307647] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0250924] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0225947] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0130143] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0117456] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0103136] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0103136] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00802044] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.0070551] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00362924] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00321222] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00275335] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00275335] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00202028] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.00172641] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000769043] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000663547] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000550939] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000550939] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000329682] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0.000253892] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [6.59106e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [5.14903e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.77056e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.77056e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.99253e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.44138e-05] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [2.71958e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [2.0046e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.36407e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.36407e-06] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [5.21191e-07] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.19072e-07] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [2.47419e-08] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.54021e-08] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [8.44038e-09] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [8.44038e-09] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.04504e-09] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.49609e-10] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [7.51767e-13] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [2.16789e-13] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [4.20522e-14] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [4.20522e-14] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [4.20522e-14] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [1.17796e-16] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.30553e-18] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.8154e-35] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [3.24338e-45] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0] , 8.24347
2.1114 , 0.289429 , 0.61694 , 0.197171 , [0] , 8.24347

The integral (not shown) does match the R-function, and the exposed stan-function. I am still not sure how to properly interpret this, and how this might be related to the error messages. I managed to run the integrate_ode_rk45 version of the model for I=2 items and J=10 test takers for 200 iterations (but only once I installed the Rstan-developer version), although I am still getting the

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: Exception: Exception: lognormal_lpdf: Random variable is -8.88178e-16, but must be nonnegative!. 

error message sporadically during warmup and during sampling.

If we do this like a missing data model, where we’re going to estimate r_1 (the first recation time), then we have something like:

r_1 \sim LN(\mu_1, \sigma_1) \\ r_t - r_1 \sim LN(\mu_2, \sigma_2)

r_2 = r_t - r_1 here. So if we want to integrate out r_1 that means doing:

p(r_2 | \theta) = \int_0^{r_t} p(r_1 | \theta) p(r_2 | \theta) d r_1

There are limits on this integral because I think there are limits on the valid values of r_1 (can’t be lower than zero, can’t be above r_t). If you get numerical problems integrating from zero, I don’t see why you couldn’t just make the assumption that r_1 > \epsilon where \epsilon is some small value that nobody is gonna finish test 1 faster than anyway.

So we should be able to compute that with something that looks very similar to what you have for the convolve kernel, the only other thing I noticed thinking about this is we better accumulate the log of the integral, so we want:

target += log(integrate_out_r1(...));

not just the result of the integral (the log looks like it’s missing in some of the earlier convolve posts, dunno if that changed).

1 Like

Thanks again for this valuable suggestion Ben! Would the model specified that way look like this:

functions{ 
  
  real[] kernel(real x, real[] y, real[] theta, data real[] x_r, data int[] x_i) {
    real dydx[1];
    real mu1    = theta[1];     
    real sigma1 = theta[2];    
    real mu2    = theta[3];  
    real sigma2 = theta[4];
    real rt1    = theta[5];
    real rt2    = theta[6];

    dydx[1] = exp( lognormal_lpdf( rt1 | mu1, sigma1) + lognormal_lpdf( rt2 | mu2, sigma2) );
    return dydx;
  }
  
  real integrate_out_r1(real a, real b, real[] theta, data real x, real rel_tol) {
    int x_i[0];
    real out1;
    out1 = integrate_ode_rk45(kernel, rep_array(0.,1), a, rep_array(b,1), theta, rep_array(x,1), x_i)[1,1];
    return out1;
  }
}

data{
  int<lower = 1> I;                      // number of items
  int<lower = 1> J;                      // number of examinees
  int<lower = 1> N;                      // number of observed data points
  int<lower = 1> ii[N];                  // item id
  int<lower = 1> jj[N];                  // person id
  real<lower=0> rt[N];                   // response times
}

parameters{
  vector[J] omega1;
  vector[J] omega2;
  vector[I] alpha1;
  vector[I] alpha2;
  positive_ordered[2] phi;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
  vector<lower=0.5>[N] rt1;
}

model{
  vector[N] mu1 = alpha1[ii] - phi[2] * omega1[jj];
  vector[N] mu2 = alpha2[ii] - phi[1] * omega2[jj];
  vector[N] rt2 = to_vector(rt) - rt1;
  
  omega1 ~ std_normal();
  omega2 ~ std_normal();
  
  alpha1 ~ normal(1,0.5);
  alpha2 ~ normal(1,0.5);
  
  phi ~ normal(0,1);
  
  sigma1 ~ exponential(5);
  sigma2 ~ exponential(5);
  
  rt1 ~ lognormal(mu1, sigma1[ii]);
  rt2 ~ lognormal(mu2, sigma2[ii]);
  
  for(n in 1:N) target += log( integrate_out_r1(0, rt[n], {mu1[n], sigma1[ii[n]], mu2[n], sigma2[ii[n]], rt1[n], rt2[n]}, rt[n], 1E-6) );
  
}

It compiles and samples fairly quickly, but I am afraid there is an error somewhere, because the results do not make sense yet. Yes, the log was added to the target += command in the current model specifications, so we accumulate the log of the integral now.

That is completely reasonable, finishing an item under 1 second, for example, would indicate an omitted item and would also be unrealistic. I am unsure where to implement this assumption, since in the previous convolve model specifications I don’t have r_1 as parameter, only implicitly in the kernel function.

The first model here is you sample rt1 like a parameter:

data{
  int<lower = 1> I;                      // number of items
  int<lower = 1> J;                      // number of examinees
  int<lower = 1> N;                      // number of observed data points
  int<lower = 1> ii[N];                  // item id
  int<lower = 1> jj[N];                  // person id
  vector<lower=0>[N] rt;                   // response times
}

parameters {
  vector[J] omega1;
  vector[J] omega2;
  vector[I] alpha1;
  vector[I] alpha2;
  vector[2] phi;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
  vector<lower=0.5, upper = rt>[N] rt1;
}

model{
  vector[N] mu1 = alpha1[ii] - phi[1] * omega1[jj];
  vector[N] mu2 = alpha2[ii] - phi[2] * omega2[jj];
  
  omega1 ~ std_normal();
  omega2 ~ std_normal();
  
  alpha1 ~ normal(2.0, 1.0);
  alpha2 ~ normal(2.0, 1.0);
  
  phi ~ normal(0,1);
  
  sigma1 ~ normal(0, 1);
  sigma2 ~ normal(0, 1);
}

Try that first and see if you get anywhere with it. I didn’t try this, so it probably has bugs.

I got rid of the positive ordered constraint on phi cause I wasn’t sure where it was coming from and I change some other things just to make it match the code I was using for data generation. Still problems. I think it would be useful to simplify this model more to work out how integrating stuff out will work.

And then in the second you integrate it out. Something like this. I wasn’t able to get this working so it definitely has problems somewhere:

functions{ 
  real integrand(real x, real xc, real[] theta, data real[] x_r, data int[] x_i) {
    real mu1    = theta[1];     
    real sigma1 = theta[2];    
    real mu2    = theta[3];  
    real sigma2 = theta[4];
    real rt    = x_r[x_i[1]];

    return exp( lognormal_lpdf( x | mu1, sigma1) + lognormal_lpdf( rt - x | mu2, sigma2) );
  }
  
  vector integrand2(real t, vector x, real mu1, real sigma1, real mu2, real sigma2, data real rt) {
    return [ exp( lognormal_lpdf( t | mu1, sigma1) + lognormal_lpdf( rt - t | mu2, sigma2) ) ]';
  }
}

data{
  int<lower = 1> I;                      // number of items
  int<lower = 1> J;                      // number of examinees
  int<lower = 1> N;                      // number of observed data points
  int<lower = 1> ii[N];                  // item id
  int<lower = 1> jj[N];                  // person id
  real<lower=0> rt[N];                   // response times
}

parameters {
  vector[J] omega1;
  vector[J] omega2;
  vector[I] alpha1;
  vector[I] alpha2;
  vector[2] phi;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
}

model{
  vector[N] mu1 = alpha1[ii] - phi[1] * omega1[jj];
  vector[N] mu2 = alpha2[ii] - phi[2] * omega2[jj];
  
  omega1 ~ std_normal();
  omega2 ~ std_normal();
  
  alpha1 ~ normal(2.0, 1.0);
  alpha2 ~ normal(2.0, 1.0);
  
  phi ~ normal(0,1);
  
  sigma1 ~ normal(0, 1);
  sigma2 ~ normal(0, 1);
  
  for(n in 1:N) {
    vector[1] y0 = rep_vector(0.0, 1);
    vector[1] out[1] = ode_rk45(integrand2, y0, 0.01 * rt[n], { 0.99 * rt[n] }, mu1[n], sigma1[ii[n]], mu2[n], sigma2[ii[n]], rt[n]);
    target += log(out[1, 1]);
    //target += log(integrate_1d(integrand, 0.01 * rt[n], 0.99 * rt[n], { mu1[n], sigma1[ii[n]], mu2[n], sigma2[ii[n]] }, rt, { n }));
  }
}

I left the integrate_1d and ode_rk45 versions. Not sure why integrate_1d was being super flaky. Probably indicates a bug. (Edit: I hope a bug in the model though maybe something is still messed up with integrate_1d. If you get a chance to isolate a set of problematic parameters that would be useful)

Dear All,

After quite some while I revisited this problem and came up with an alternative specification based on the Fenton-Wilkinson approximation of the sum of i.i.d lognormal random variables. The basic idea of the model remains the same: the observed response time rt_{ji} of test taker j (j = 1,2,...,J) to item i i = 1,2,...I) is the sum of two unobserved response time components rt_{1ji} = \alpha_{1i} + \phi_{1i} * \omega_{1j} + \epsilon_{1i} and rt_{2ji} = \alpha_{2i} + \phi_{2i} * \omega_{2j} + \epsilon_{2i}. Each response time component is assumed to follow a lognormal distribution with mean \alpha_{i} + \phi_{i} * \omega_{j} and scale \epsilon_{i}. This is the three-parameter lognormal model by Klein-Entink et al. (2009).

I wrote a custom function for the Fenton-Wilkinson approximation and integrated it into the model code (see below).

functions {
  row_vector fenton1(vector m, vector s) {
    
    int D = dims(s)[1];
    matrix[D,D] term1;
    matrix[D,D] term2;
    matrix[D,D] d;
    vector[D] d1;
    real u1;
    real u2;
    row_vector[3] out;
    
    d = diag_matrix(square(s));
    d1 = diagonal(d);
    
    for(i in 1:D) {
      for(j in 1:D) { 
        term1[i,j] = m[i] + m[j];
        term2[i,j] = d1[i] + d1[j]; 
      }
    }
    
    u1 = sum(exp(m + 0.5 * d1));
    u2 = sum(exp(term1 + 0.5 * term2 + d));
    out = [2*log(u1)-0.5*log(u2), log(u2)-2*log(u1), sqrt(log(u2)-2*log(u1))];
    
    return(out);
    
  }
}

data {
  int<lower=1> I;                      // number of items
  int<lower=1> J;                      // number of examinees
  int<lower=1> N;                      // number of observed data points
  matrix<lower=0>[J,I] rt;             // response times
}

parameters {
  vector[J] speed1;
  vector[J] speed2;
  vector<lower=0>[I] intensity1;
  vector<lower=0>[I] intensity2;
  vector<lower=0>[I] phi1_raw;
  vector<lower=0>[I] phi2;
  vector<lower=0>[I] sigma1;
  vector<lower=0>[I] sigma2;
}

transformed parameters {
  // original identification constraint for the 3PLN (Klein-Entink, 2009)
  vector[I] phi1 = exp(phi1_raw - mean(phi1_raw));
}

model {
   // Fenton-Wilkinson Approximation
   // LP1 and LP2 are only containers for the individual means
   // mus and sds are containers for the distribution of their sum
   matrix[J,I] LP1;
   matrix[J,I] LP2;
   matrix[J,I] mus;
   matrix[J,I] sds;
  
   for(j in 1:J) {
     for(i in 1:I) {
       LP1[j,i] = intensity1[i] + phi1[i] * speed1[j];
       LP2[j,i] = intensity2[i] + phi2[i] * speed2[j];
       mus[j,i] = fenton1([LP1[j,i],LP2[j,i]]',[sigma1[i], sigma2[i]]')[1];
       sds[j,i] = fenton1([LP1[j,i],LP2[j,i]]',[sigma1[i], sigma2[i]]')[3];
     }
   }

  // prior structure
  speed1 ~ normal(0,1);
  speed2 ~ normal(0,1);
  intensity1 ~ normal(0,2);
  intensity2 ~ normal(0,2);
  phi1_raw ~ normal(0,0.5);
  phi2 ~ normal(0,0.5);
  sigma1 ~ exponential(1);
  sigma2 ~ exponential(1);
  
  // Likelihood
  for(j in 1:J) {
   for(i in 1:I) {
     rt[j,i] ~ lognormal(mus[j,i], sds[j,i]);
   }
  }
}

I am fairly certain that the approximation works, as I tested it with the expose_stan_functions tool prior to integrating it into the model code (following Cobb and Rumi, 2012), as well as within the model with the true parameters as input.

Contrary to the previous model specifications, this time the sampler runs fine more often (good diagnostics) and provides reasonable parameter estimates (well, most of the time, and reasonable means that the estimates are in the range of the true parameter values; they are still biased). I am running the sampler with two parallel chains with 1800 iterations (600 being warmup), adapt_delta = 0.9 and max_treedepth = 11. The model takes quite some time to run, so debugging is a bit hard, as I can only use a small number of items and test takers without having to wait a couple hours for the results of one run. So I guess that the bias is mostly due to the small number of items and test takers.

Moreover, I am getting many normal_lpdf: Location parameter is -inf, but must be finite! and normal_lpdf: Location parameter is nan, but must be finite! warnings during warmup and sampling, so I am wondering if there is an error in the way I integrated the Fenton-Wilkinson approximation into my model code. The warnings get fewer when I increase the number of items and test takers, so I guess this is also due to the (for the complexity of the model) inadequate number of items and test takers.

So this time, my questions are as follows:
(1) Did I integrate the Fenton-Wilkinson approximation correctly into my model code (are there obvious errors)?
(2) Are there any possibilities to speed up the sampling by optimizing the code?
(3) Where do the aforementioned warnings come from, and how can I eliminate them? After looking through the discourse, the reason might be that both LP1 and LP2, and mus and sds are not assigned values? But it is unclear to me how to do this, as they are only containers.

As always, any help/advice is highly appreciated (I would like to move on to a more extensive simulation with adequate numbers of items and test takers)! I have attached code to generate some toy data to this post.

Best wishes
Chris

data_gen.R (797 Bytes)