Hi,
I am trying to include samples from my prior in my brms model output and am running into an error. When I set sample_prior = “yes”, I receive an error that prevents compilation. Everything runs smoothly when I choose sample_prior = “no” or “only”.
The error is:
Info: left-hand side variable (name=circSD) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
Info: left-hand side variable (name=theta) occurs on right-hand side of assignment, causing inefficient deep copy to avoid aliasing.
variable definition base type mismatch, variable declared as base type: vector variable definition has base: real error in 'model158167cb2407_file158135721b25' at line 155, column 54
-------------------------------------------------
153: real prior_b_circSD_1 = normal_rng(4,0.5);
154: real prior_b_circSD_2 = normal_rng(0,0.5);
155: vector[K_theta] prior_b_theta_1 = normal_rng(0,1.5);
^
156: vector[K_theta] prior_b_theta_2 = normal_rng(0,1);
-------------------------------------------------
Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname, :
failed to parse Stan model 'file158135721b25' due to the above error.
The model is a two-level uniform + von mises mixture I created using custom_family. The model code is:
#definition
vm_uniform_mix <- custom_family(
"vm_uniform_mix",
dpars = c("mu", "circSD", "theta", "a", "b"),
links = c("identity", "log", "logit", "identity", "identity"),
lb = c(NA, 0, 0, NA, NA),
ub = c(NA, NA, 1, NA, NA),
type = "real" # response type
)
# fill in functions needed to density evaluation and rng
stan_funs <- "
real deg2rad( real d){
return(d * (pi()/180));
}
real sd2k( real sd_input) {
//sd_input needs to be radians
//original matlab code
// R <- exp(-S.^2/2);
// K = 1./(R.^3 - 4 * R.^2 + 3 * R);
// K(R < 0.85) = -0.4 + 1.39 * R(R < 0.85) + 0.43./(1 - R(R < 0.85));
// K(R < 0.53) = 2 * R(R < 0.53) + R(R < 0.53).^3 + (5 * R(R < 0.53).^5)/6;
real R;
real K;
R = exp(-(sd_input^2)/2);
K = 1/((R^3) - (4 * R^2) + (3*R));
if(R < 0.85){
K = -0.4 + (1.39 * R) + (0.43/(1-R));
}
if(R < 0.53){
K = (2 * R) + R^3 + (5 * R^5)/6;
}
return(K);
}
real vm_uniform_mix_lpdf(real y, real mu, real circSD, real theta, real a, real b) {
real kappa = sd2k(deg2rad(circSD));
//* for kappa > 100 the normal approximation is used
//* for reasons of numerial stability
if (kappa > 100){
return log_mix( theta ,
normal_lpdf(y | mu, sqrt(1/kappa)),
uniform_lpdf(y | a, b));
} else {
return log_mix( theta ,
von_mises_lpdf(y | mu, kappa),
uniform_lpdf(y | a, b));
}
}
real vm_uniform_mix_rng(real mu, real circSD, real theta, real a, real b) {
int flip = bernoulli_rng(theta);
real phi;
if (flip == 1){
phi = von_mises_rng(mu, sd2k(deg2rad(circSD)));
}else{
phi = uniform_rng(a, b);
}
return(phi);
}
"
stanvars <- stanvar(scode = stan_funs, block = "functions")
# fitting
bform <- bf(error ~ 0, #formula for first parameter mu
circSD ~ 0 + intercept + stimulation + (1 + stimulation || subj_index),
theta ~ 0 + intercept + stimulation + (1 + stimulation || subj_index),
a = -pi,
b = pi,
family = vm_uniform_mix)
# priors
bprior <- prior(normal(4, 0.5), class = "b", coef = "intercept", dpar = "circSD") +
prior(normal(0, 0.5), class = "b", coef = "stimulation", dpar = "circSD") +
prior(normal(0, 0.5), class = "sd", coef = "Intercept", dpar = "circSD", group = "subj_index") +
prior(normal(0, 0.5), class = "sd", coef = "stimulation", dpar = "circSD", group = "subj_index") +
prior(normal(0, 1.5), class = "b", coef = "intercept", dpar = "theta") +
prior(normal(0, 1), class = "b", coef = "stimulation", dpar = "theta") +
prior(normal(0, 1.5), class = "sd", coef = "Intercept", dpar = "theta", group = "subj_index") +
prior(normal(0, 1), class = "sd", coef = "stimulation", dpar = "theta", group = "subj_index")
model_fit <- brm(bform, obs_only, family = vm_uniform_mix, prior = bprior, stanvars = stanvars,
sample_prior = "yes",
warmup = 1000, iter = 2000, cores = 4, chains = 4,
control = list(adapt_delta = 0.99), inits = 0,
file = "model_fit")
Any help figuring this out would be greatly appreciated! Let me know if I can provide any more useful info.
Best,
Chris
- Operating System: macOS Mojave 10.14.6
- brms Version: 2.9.0