I was able to create the model successful but I can’t get pp_check to work on the stanfit object. “could not find gamma_rng”.
To get the succesful model, however, I had to move the custom gamma2 function in the function block to the top of the body to remove parsing errors saying gamma2 does not exist. This same error is marked in my stancode file for partial_log_lik and shows up when trying to call expose_stan_functions().
Variable "partial_log_lik" does not exist.
error in 'model45ad79d125e8_User_defined_functions' at line 77, column 40
-------------------------------------------------
75: // likelihood including all constants
76: if (!prior_only) {
77: target += reduce_sum(partial_log_lik, seq, grainsize, Y, Xc, b, Intercept, Xc_v, b_v, Intercept_v);
^
78: }
-------------------------------------------------
Error in stanc(model_code = mc, model_name = "User-defined functions", :
failed to parse Stan model 'User-defined functions' due to the above error.
// generated with brms 2.13.11
functions {
real gamma2_lpdf(real y, real mu, real v) {
return gamma_lpdf(y | mu * mu / v, mu / v);
}
real gamma2_rng(real mu, real v) {
return gamma_rng(mu * mu / v, mu / v);
}
int[] sequence(int start, int end) {
int seq[end - start + 1];
for (n in 1:num_elements(seq)) {
seq[n] = n + start - 1;
}
return seq;
}
// compute partial sums of the log-likelihood
real partial_log_lik(int[] seq, int start, int end, vector Y, matrix Xc, vector b, real Intercept, matrix Xc_v, vector b_v, real Intercept_v) {
real ptarget = 0;
int N = end - start + 1;
// initialize linear predictor term
vector[N] mu = Intercept + Xc[start:end] * b;
// initialize linear predictor term
vector[N] v = Intercept_v + Xc_v[start:end] * b_v;
for (n in 1:N) {
// apply the inverse link function
v[n] = exp(v[n]);
}
for (n in 1:N) {
// apply the inverse link function
mu[n] = exp(mu[n]);
}
for (n in 1:N) {
int nn = n + start - 1;
ptarget += gamma2_lpdf(Y[nn] | mu[n], v[n]);
}
return ptarget;
}
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int<lower=1> K_v; // number of population-level effects
matrix[N, K_v] X_v; // population-level design matrix
int grainsize; // grainsize for threading
int prior_only; // should the likelihood be ignored?
}
transformed data {
int Kc = K - 1;
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
int Kc_v = K_v - 1;
matrix[N, Kc_v] Xc_v; // centered version of X_v without an intercept
vector[Kc_v] means_X_v; // column means of X_v before centering
int seq[N] = sequence(1, N);
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
for (i in 2:K_v) {
means_X_v[i - 1] = mean(X_v[, i]);
Xc_v[, i - 1] = X_v[, i] - means_X_v[i - 1];
}
}
parameters {
vector[Kc] b; // population-level effects
real Intercept; // temporary intercept for centered predictors
vector[Kc_v] b_v; // population-level effects
real Intercept_v; // temporary intercept for centered predictors
}
transformed parameters {
}
model {
// likelihood including all constants
if (!prior_only) {
target += reduce_sum(partial_log_lik, seq, grainsize, Y, Xc, b, Intercept, Xc_v, b_v, Intercept_v);
}
// priors including all constants
target += normal_lpdf(b | 0, 3);
target += normal_lpdf(Intercept | 3, 1);
target += normal_lpdf(b_v | 0, 3);
target += normal_lpdf(Intercept_v | 3, 1);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
// actual population-level intercept
real b_v_Intercept = Intercept_v - dot_product(means_X_v, b_v);
// additionally draw samples from priors
real prior_b = normal_rng(0,3);
real prior_Intercept = normal_rng(3,1);
real prior_b_v = normal_rng(0,3);
real prior_Intercept_v = normal_rng(3,1);
}
![image](https://global.discourse-cdn.com/standard14/uploads/mc_stan/original/2X/1/1ae17fc0d5f28824bc51866b0a5525fe68593fae.png)
is there a workaround? a Cpp workaround?