Hello,
I have written an IRT model which involves a user-defined probability function. To compute the latter, I need to integrate the exponential of a spline. I accomplished this via “integrate_ode_rk45” and it compiled and worked.
However, my first solution involved a lot of loops and I tried to refactor the code to use vectorized expressions instead of loops. Now the new model (see below) does compile when calling “stanc”. However, when I try to fit the model, I get an error message involving the following error message (translated from German):
"Execution of command ‘D:/R-4.0.2/bin/x64/R CMD SHLIB file3588724c109.cpp 2> file3588724c109.cpp.err.txt’ resulted in Status 1".
Via uncommenting and changing some statements in the underlying stan code, I was able to narrow the error down to the line, wherein I call the integrate_ode_rk45 function. However, I have no idea why this line causes the error.
Thank you for any advices/sugesstions on how to resolve this problem!
functions
{
vector build_b_spline(real[] t, real[] ext_knots, int ind, int order);
vector build_b_spline(real[] t, real[] ext_knots, int ind, int order)
{
vector[size(t)] b_spline;
vector[size(t)] w1 = rep_vector(0, size(t));
vector[size(t)] w2 = rep_vector(0, size(t));
if (order==1)
for (i in 1:size(t))
b_spline[i] = (ext_knots[ind] <= t[i]) && (t[i] < ext_knots[ind+1]);
else {
if (ext_knots[ind] != ext_knots[ind+order-1])
w1 = (to_vector(t) - rep_vector(ext_knots[ind], size(t))) /
(ext_knots[ind+order-1] - ext_knots[ind]);
if (ext_knots[ind+1] != ext_knots[ind+order])
w2 = 1 - (to_vector(t) - rep_vector(ext_knots[ind+1], size(t))) /
(ext_knots[ind+order] - ext_knots[ind+1]);
b_spline = w1 .* build_b_spline(t, ext_knots, ind, order-1) +
w2 .* build_b_spline(t, ext_knots, ind+1, order-1);
}
return b_spline;
}
real[] exp_spline_ode(real time, real[] state, real[] theta, real[] x_r, int[] x_i)
{
real dydt[1];
real first[1];
first[1] = time;
dydt[1] = exp(theta[1]*build_b_spline(first, x_r, x_i[1], x_i[2])[1]);
return dydt;
}
real espl_lpdf(vector y, real[] theta, real coeff, data real[] ext_knots)
{
real loc[rows(y)];
real deriv_spline[rows(y)];
real ruck[rows(y)];
int x_i[2];
real times[rows(y)];
real theta_h[1];
real init_st[rows(y)] = rep_array(0.0, rows(y));
//real init_st[1] = {0.0};
x_i[1] = 1;
x_i[2] = 2;
times = to_array_1d(y);
theta_h[1] = coeff;
deriv_spline = to_array_1d(build_b_spline(to_array_1d(y), ext_knots, 1, 2));
//This call is causing the trouble!
loc = (integrate_ode_rk45(exp_spline_ode, init_st, ext_knots[1], times, theta_h, ext_knots, x_i))[,1];
return sum(std_normal_lpdf(to_vector(loc)-to_vector(theta))+coeff*to_vector(deriv_spline));
}
}
data {
int<lower=0> K; // number of items
int<lower=1> N; // number of test takers
matrix[N, K] Y; // responses of test takers to the items
int num_knots;
real ext_knots[num_knots];
}
parameters {
vector[N] z;
vector[K] intercepts;
vector[K] coeffs;
real<lower=0> sigma_theta;
}
transformed parameters {
vector[N] theta = sigma_theta*z;
}
model {
sigma_theta ~ cauchy(0,1);
intercepts ~ normal(0,1);
coeffs ~ normal(0,1);
z ~ normal(0, 1);
for(j in 1:K)
{
Y[,j] ~ espl(to_array_1d(theta-intercepts[j]),coeffs[j],ext_knots);
}
}