Hi,
I’m trying to implement reduce_sum with a graded response model but can’t get it to compile.
I have read the Reduce-sum section in the Stan User’s Guide, related topics here on discourse, such as https://discourse.mc-stan.org/t/testing-reduce-sum/14376, and @martinmodrak has already provided some helpful comments here.
Below is the model without reduce_sum using ordered_logistic_lpmf. This compiles, runs without warning messages, and produces the expected results. Data is attached. stan_data.Rdata (6.3 KB)
I have cmdstan version 2.27.0 installed.
library(cmdstanr)
load("./stan_data.Rdata")
// GRADED RESPONSE MODEL - assign as gIRT
data {
int<lower=1> I; // Number of items
int<lower=1> J; // Number of respondents
int<lower=1> N; // Total responses
int<lower=1,upper=I> ii[N]; // Item ID
int<lower=1,upper=J> jj[N]; // Person ID
int<lower=0> y[N]; // Responses
int<lower=1> K; // Number of covariates for location (mu) of ability (theta)
matrix[J,K] X; // Covariate matrix (including intercept) for location (mu) of ability (theta)
}
transformed data {
int r[N]; // modified response; r = 1, 2, ... m_i + 1
int m[I]; // number of difficulty parameters per item
int pos_alpha[I]; // first position in difficulty for item
int pos_delta[I]; // first position in between-step change for item
// Allowing for differing numbers of response categories
m = rep_array(0, I);
for(n in 1:N) {
r[n] = y[n] + 1;
if(y[n] > m[ii[n]]) m[ii[n]] = y[n];
}
pos_alpha[1] = 1;
pos_delta[1] = 1;
for(i in 2:(I)) {
pos_alpha[i] = pos_alpha[i-1] + m[i-1];
pos_delta[i] = pos_delta[i-1] + m[i-1] - 1;
}
}
parameters {
vector[I] eta; // First diffulty parameter for each item
vector<lower=0>[sum(m)-I] delta; // Between-step changes in difficulty
real<lower = 0> beta[I]; // Constrain all discrimination parameters to be positive
vector[J] theta; // Ability parameter
vector[K] gamma; // Covariate coefficients
}
transformed parameters {
vector[sum(m)] alpha; // Composite item difficulties
for(i in 1:I) {
if (m[i] > 0)
alpha[pos_alpha[i]] = eta[i];
for (k in 2:m[i])
alpha[pos_alpha[i]+k-1] = alpha[pos_alpha[i]+k-2] + delta[pos_delta[i]+k-2];
}
}
model{
eta ~ normal(0,1);
delta ~ normal(0,1);
beta ~ lognormal(1,1);
gamma ~ student_t(3, 0, 1);
target += normal_lpdf(alpha | 0, 1);
theta ~ normal(X*gamma,1);
mean(theta) ~ normal(0, 0.001); // Mean of ability constrained to 0 for identification
for (n in 1:N) {
target += ordered_logistic_lpmf(r[n] | theta[jj[n]] * beta[ii[n]],
segment(alpha, pos_alpha[ii[n]], m[ii[n]]));
}
}
"
write(gIRT, "gIRT.stan") # save above model
compiled_gIRT <- cmdstan_model(stan_file = "gIRT.stan") # compile
compiled_gIRT$check_syntax()
fit_gIRT <- compiled_gIRT$sample(
data = stan_data,
chains = 4,
parallel_chains = 4,
iter_warmup = 500,
iter_sampling = 100,
refresh = 50
)
Following from the docs, I’d have expected this reduce_sum implementation to work (just the functions and model blocks):
functions {
real partial_sum_lpmf(int[] slice_r,
int start, int end,
int[] ii,
int[] jj,
vector theta,
real beta,
vector alpha,
int[] m,
int[] pos_alpha) {
return ordered_logistic_lupmf(slice_r |
theta[jj[start:end]] * beta[ii[start:end]],
segment(alpha, pos_alpha[ii[start:end]],
m[ii[start:end]]));
}
}
model{
eta ~ normal(0,1);
delta ~ normal(0,1);
beta ~ lognormal(1,1);
gamma ~ student_t(3, 0, 1);
target += normal_lpdf(alpha | 0, 1);
theta ~ normal(X*gamma,1);
mean(theta) ~ normal(0, 0.001); // Mean of ability constrained to 0 for identification
target += reduce_sum(partial_sum_lupmf, r, grainsize,
ii, jj, theta, beta, alpha, m, pos_alpha);
}
which throws the following error:
Semantic error in line 13, column 55 to column 74:
-------------------------------------------------
11: int[] pos_alpha) {
12: return ordered_logistic_lupmf(slice_r |
13: theta[jj[start:end]] * beta[ii[start:end]],
^
14: segment(alpha, pos_alpha[ii[start:end]],
15: m[ii[start:end]]));
-------------------------------------------------
Too many indexes, expression dimensions=0, indexes found=1
Adjusting partial_sum_lpmf along suggestions by @martinmodrak here:
functions {
real partial_sum_lpmf(int[] slice_r,
int start, int end,
int[] ii,
int[] jj,
vector theta,
real[] beta,
vector alpha,
int[] m,
int[] pos_alpha) {
return ordered_logistic_lupmf(slice_r |
theta[jj[start:end]] * to_row_vector(beta[ii[start:end]]),
segment(alpha, pos_alpha[ii[start:end]],
m[ii[start:end]]));
}
}
model{
eta ~ normal(0,1);
delta ~ normal(0,1);
beta ~ lognormal(1,1);
gamma ~ student_t(3, 0, 1);
target += normal_lpdf(alpha | 0, 1);
theta ~ normal(X*gamma,1);
mean(theta) ~ normal(0, 0.001); // Mean of ability constrained to 0 for identification
target += reduce_sum(partial_sum_lupmf, r, grainsize,
ii, jj, theta, beta, alpha, m, pos_alpha);
}
yields this error message:
Semantic error in line 14, column 32 to line 15, column 65:
-------------------------------------------------
12: return ordered_logistic_lupmf(slice_r |
13: theta[jj[start:end]] * to_row_vector(beta[ii[start:end]]),
14: segment(alpha, pos_alpha[ii[start:end]],
^
15: m[ii[start:end]]));
16: }
-------------------------------------------------
Ill-typed arguments supplied to function 'segment'. Available signatures:
(vector, int, int) => vector
(row_vector, int, int) => row_vector
(array[] int, int, int) => array[] int
(array[] real, int, int) => array[] real
(array[] vector, int, int) => array[] vector
(array[] row_vector, int, int) => array[] row_vector
(array[] matrix, int, int) => array[] matrix
(array[,] int, int, int) => array[,] int
(array[,] real, int, int) => array[,] real
(array[,] vector, int, int) => array[,] vector
(array[,] row_vector, int, int) => array[,] row_vector
(array[,] matrix, int, int) => array[,] matrix
(array[,,] int, int, int) => array[,,] int
(array[,,] real, int, int) => array[,,] real
(array[,,] vector, int, int) => array[,,] vector
(array[,,] row_vector, int, int) => array[,,] row_vector
(array[,,] matrix, int, int) => array[,,] matrix
Instead supplied arguments of incompatible type: vector, array[] int, array[] int
Alternatively, I also tried replacing segment:
functions {
real partial_sum_lpmf(int[] slice_r,
int start, int end,
int[] ii,
int[] jj,
vector theta,
real[] beta,
vector alpha,
int[] m,
int[] pos_alpha) {
return ordered_logistic_lupmf(slice_r |
theta[jj[start:end]] * to_row_vector(beta[ii[start:end]]),
alpha[pos_alpha[ii[start:end]]:(pos_alpha[ii[start:end]] + m[ii[start:end]])]);
}
}
model{
eta ~ normal(0,1);
delta ~ normal(0,1);
beta ~ lognormal(1,1);
gamma ~ student_t(3, 0, 1);
target += normal_lpdf(alpha | 0, 1);
theta ~ normal(X*gamma,1);
mean(theta) ~ normal(0, 0.001); // Mean of ability constrained to 0 for identification
target += reduce_sum(partial_sum_lupmf, r, grainsize,
ii, jj, theta, beta, alpha, m, pos_alpha);
}
which produces:
Semantic error in line 14, column 38 to column 62:
-------------------------------------------------
12: return ordered_logistic_lupmf(slice_r |
13: theta[jj[start:end]] * to_row_vector(beta[ii[start:end]]),
14: alpha[pos_alpha[ii[start:end]]:(pos_alpha[ii[start:end]] + m[ii[start:end]])]);
^
15: }
16: }
-------------------------------------------------
Range bound must be of type int. Instead found type array[] int
Any help on this would be much appreciated! Maybe @andrjohns and @wds15 have some suggestions?