I include an example, below, but want to stress that essentially this is just one example of many scripts that successfully compiled and estimated under cmdstan 2.28 that no longer do under 2.29.1.
functions{
vector build_b_spline(vector t, vector ext_knots, int ind, int order);
matrix build_mux(int N, int start, int end, int K_sp, matrix X, matrix[] G,
matrix beta_x, matrix[] beta_w);
row_vector build_muxi(int K_sp, int num_basis, row_vector x_i, vector[] g_i,
matrix beta_x, matrix[] beta_w);
vector build_b_spline(vector t, vector ext_knots, int ind, int order) {
// INPUTS:
// t: the points at which the b_spline is calculated
// ext_knots: the set of extended knots
// ind: the index of the b_spline
// order: the order of the b-spline
vector[num_elements(t)] b_spline;
vector[num_elements(t)] w1 = rep_vector(0, num_elements(t));
vector[num_elements(t)] w2 = rep_vector(0, num_elements(t));
if (order==1)
for (i in 1:num_elements(t)) // B-splines of order 1 are piece-wise constant
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], num_elements(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], num_elements(t))) /
(ext_knots[ind+order] - ext_knots[ind+1]);
// Calculating the B-spline recursively as linear interpolation of two lower-order splines
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;
}
matrix build_mux(int N, int start, int end, int K_sp, matrix X, matrix[] G, matrix beta_x, matrix[] beta_w){
matrix[N,2] mu_x;
for( arm in 1:2 )
{
mu_x[1:N,arm] = X[start:end,] * to_vector(beta_x[,arm]); /* N x l for each arm */
// spline term
for( k in 1:K_sp )
{
mu_x[1:N,arm] += to_vector(beta_w[arm][,k]' * G[k][,start:end]); /* N x 1 */
} /* end loop k over K predictors */
}/* end loop arm over convenience and reference sample arms */
return mu_x;
}
row_vector build_muxi(int K_sp, int num_basis,
row_vector x_i, vector[] g_i, matrix beta_x, matrix[] beta_w){
row_vector[2] mu_xi;
for( arm in 1:2 )
{
mu_xi[arm] = dot_product(x_i,beta_x[,arm]); /* scalar */
// spline term
for( k in 1:K_sp )
{
mu_xi[arm] += dot_product(beta_w[arm][1:num_basis,k], g_i[k][1:num_basis]); /* scalar */
} /* end loop k over K predictors */
}/* end loop arm over convenience and reference sample arms */
return mu_xi;
} /* end function build_mu */
real partial_sum(int[] s,
int start, int end, real[] logit_pw, int K_sp, int n_c, int n,
int num_basis, matrix X, matrix[] G, matrix beta_x, matrix[] beta_w,
real phi_w) {
int N = end - start + 1;
matrix[N,2] mu_x;
matrix[N,2] p;
vector[N] p_tilde; // this pseudoprobability must be in [0,1]
real fred = 0;
// memo: slicing on all of mu_x[li,arm], p[li,arm], p_tilde[li] for li in 1:(end-start+1)
// where p_tilde is the mean vector for binary data vector, s, and mu_x[,2]
// is the mean vector for data vector logit_pw.
// Also slicing data vectors s and logit_pw in their respective
// log-likelihood contributions.
// for( li in 1:N)
// {
// int i = li + start - 1;
// vector[num_basis] g_i[K_sp] = G[1:K_sp,1:num_basis,i];
// mu_x[li,] = build_muxi(K_sp, num_basis, X[i,], g_i, beta_x, beta_w);
// }
mu_x = build_mux(N, start, end, K_sp, X, G, beta_x, beta_w);
p = inv_logit( mu_x );
// for( arm in 1:2 )
// {
// for( li in 1:N ) /* 'li' is slicing index */
// {
// p[li,arm] = inv_logit( mu_x[li,arm] );
// }
//
// } /* end loop i over all units */
// 1. bernoulli likelihood contribution
p_tilde = p[1:N,1] ./ ( p[1:N,1] + p[1:N,2] );
fred += bernoulli_lpmf( s[1:N] | p_tilde );
// for( li in 1:N ) // N = end - start + 1, index for slicing
// {
// p_tilde[li] = p[li,1] ./ ( p[li,1] + p[li,2] );
// fred += bernoulli_lpmf(s[li] | p_tilde[li]);
// } /* end loop i over all units */
// 2. normal likelihood contribution
// In non-threaded model, likelihood statement for n - n_c, logit_pw
// logit_pw ~ normal( mu_x[(n_c+1):n,2], phi_w );
// slicing on n length vector mu_x[,2]
// subsetting portion of mu_x[,2] linked to logit_pw
if( start > n_c ) ## all units in this chunk increment likelihood contribution for logit_pw
{
fred += normal_lpdf( logit_pw[start-n_c:end-n_c] | mu_x[1:N,2], phi_w );
// for( li in 1:N )
// {
// int i = li + start - 1;
// fred += normal_lpdf( logit_pw[i-n_c] | mu_x[li,2], phi_w );
// } /* end loop li over n_r units */
}else{ /* start <= n_c */
if( end > n_c ) /* some units in chunk < n_c and some > n_c; only those > n_c increment likelihood */
{
fred += normal_lpdf( logit_pw[1:end-n_c] | mu_x[n_c-start+2:N,2], phi_w );
// for( li in 1:N )
// {
// int i = li + start - 1;
// if( i - n_c > 0 )
// fred += normal_lpdf( logit_pw[i-n_c] | mu_x[li,2], phi_w );
// }
} /* end if statement on whether n_c \in (start,end ) */
} /* end if-else statement on whether add logit_pw likelihood contributions */
return fred;
}/* end function partial_sum() */
} /* end function block */
data{
int<lower=1> n_c; // observed convenience (non-probability) sample size
int<lower=1> n_r; // observed reference (probability) sample size
int<lower=1> N; // estimate of population size underlying reference and convenience samples
int<lower=1> n; // total sample size, n = n_c + n_r
int<lower=1> num_cores;
int<lower=1> multiplier;
int<lower=1> K; // number of fixed effects
int<lower=0> K_sp; // number of predictors to model under a spline basis
int<lower=1> num_knots;
int<lower=1> spline_degree;
matrix[num_knots,K_sp] knots;
real<lower=0> weights[n_r]; // sampling weights for n_r observed reference sample units
matrix[n_c, K] X_c; // *All* predictors - continuous and categorical - for the convenience units
matrix[n_r, K] X_r; // *All* predictors - continuous and categorical - for the reference units
matrix[n_c, K_sp] Xsp_c; // *Continuous* predictors under a spline basis for convenience units
matrix[n_r, K_sp] Xsp_r; // *Continuous* predictors under a spline basis for convenience units
int<lower=1> n_df;
} /* end data block */
transformed data{
// create indicator variable of membership in convenience or reference samples
// indicator of observation membership in the convenience sample
int grainsize = ( n / num_cores ) / multiplier;
real logit_pw[n_r] = logit(inv(weights));
int<lower=0, upper = 1> s[n] = to_array_1d( append_array(rep_array(1,n_c),rep_array(0,n_r)) );
matrix[n,K] X = append_row( X_c,X_r );
matrix[n,K_sp] X_sp = append_row( Xsp_c,Xsp_r );
/* formulate spline basis matrix, B */
int num_basis = num_knots + spline_degree - 1; // total number of B-splines
matrix[spline_degree + num_knots,K_sp] ext_knots_temp;
matrix[2*spline_degree + num_knots,K_sp] ext_knots;
matrix[num_basis,n] G[K_sp]; /* basis for model on p_c */
for(k in 1:K_sp)
{
ext_knots_temp[,k] = append_row(rep_vector(knots[1,k], spline_degree), knots[,k]);
// set of extended knots
ext_knots[,k] = append_row(ext_knots_temp[,k], rep_vector(knots[num_knots,k], spline_degree));
for (ind in 1:num_basis)
{
G[k][ind,] = to_row_vector(build_b_spline(X_sp[,k], ext_knots[,k], ind, (spline_degree + 1)));
}
G[k][num_knots + spline_degree - 1, n] = 1;
}
} /* end transformed data block */
parameters {
matrix<lower=0>[K,2] sigma2_betax; /* standard deviations of K x 2, beta_x */
/* first column is convenience sample, "c", and second column is "r" */
matrix[K,2] betaraw_x; /* fixed effects coefficients - first colum for p_c; second column for p_r */
// spline coefficients
vector<lower=0>[2] sigma2_global; /* set this equal to 1 if having estimation difficulties */
matrix<lower=0>[2,K_sp] sigma2_w;
matrix[num_basis,K_sp] betaraw_w[2]; // vector of B-spline regression coefficients for each predictor, k
// and 2 sample arms
real<lower=0> phi2_w; /* scale parameter in model for -1og(weights) */
} /* end parameters block */
transformed parameters {
matrix[K,2] beta_x;
matrix[num_basis,K_sp] beta_w[2];
matrix<lower=0>[K,2] sigma_betax;
vector<lower=0>[2] sigma_global; /* set this equal to 1 if having estimation difficulties */
matrix<lower=0>[2,K_sp] sigma_w;
real<lower=0> phi_w;
sigma_betax = sqrt( sigma2_betax );
sigma_global = sqrt( sigma2_global );
sigma_w = sqrt( sigma2_w );
phi_w = sqrt( phi2_w );
// for scale parameters for interaction effects from those for main effects to which they link
for( arm in 1:2 )
{
beta_x[,arm] = betaraw_x[,arm] .* sigma_betax[,arm]; /* Non-central parameterization */
}/* end loop arm over convenience and reference sample arms */
// spline regression coefficients
for(arm in 1:2)
{
for( k in 1:K_sp )
{
beta_w[arm][,k] = cumulative_sum(betaraw_w[arm][,k]);
beta_w[arm][,k] *= sigma_w[arm,k] * sigma_global[arm];
} /* end loop k over K predictors */
}
} /* end transformed parameters block */
model {
to_vector(sigma2_betax) ~ gamma(1,1);
to_vector(sigma2_w) ~ gamma(1,1);
sigma_global ~ gamma(1,1);
phi2_w ~ gamma(1,1);
to_vector(betaraw_x) ~ std_normal();
for(arm in 1:2)
to_vector(betaraw_w[arm]) ~ std_normal();
/* Model likelihood for y, logit_pw */
// Sum terms 1 to n in the likelihood
target += reduce_sum(partial_sum, s, grainsize,
logit_pw, K_sp, n_c, n, num_basis, X, G,
beta_x, beta_w, phi_w);
} /* end model block */
generated quantities{
matrix[n,2] p;
matrix[n,2] mu_x;
for( arm in 1:2 )
{
mu_x[,arm] = X[,] * to_vector(beta_x[,arm]); /* n x l for each arm */
// spline term
for( k in 1:K_sp )
{
mu_x[,arm] += to_vector(beta_w[arm][,k]' * G[k][,1:n]); /* n x 1 */
} /* end loop k over K predictors */
}/* end loop arm over convenience and reference sample arms */
p = inv_logit( mu_x );
// smoothed sampling weights for convenience and reference units
vector[n] weights_smooth_c = inv(p[,1]);
vector[n] weights_smooth_r = inv(p[,2]);
// inclusion probabilities in convenience and reference units for convenience units
// use for soft thresholding
vector[n_c] pi_c = p[1:n_c,1];
vector[n_c] pi_r_c = p[1:n_c,2];
// normalized weights
weights_smooth_c *= ((n_c+0.0)/(n+0.0)) * (sum(weights_smooth_r)/sum(weights_smooth_c));
weights_smooth_r *= ((n_r+0.0)/(n+0.0));
} /* end generated quantities block */