Given the model, the data and the formula for the model,
this function will return the combined stan code generated by bmm and
brms
Arguments
- object
A
bmmformulaobject- data
An object of class data.frame, containing data of all variables used in the model. The names of the variables must match the variable names passed to the
bmmodelobject for required argurments.- model
A description of the model to be fitted. This is a call to a
bmmodelsuch asmixture3p()function. Every model function has a number of required arguments which need to be specified within the function call. Callbmm_models()to see the list of supported models, and?modelname(for example?imm) for their required arguments.- prior
One or more
brmspriorobjects created bybrms::set_prior()or related functions and combined using the c method or the + operator. See alsodefault_prior()for more help. Not necessary for the default model fitting, but you can provide prior constraints to model parameters- ...
Further arguments passed to
brms::stancode(). See the description ofbrms::stancode()for more details
Examples
scode1 <- stancode(bmf(c ~ 1, kappa ~ 1),
data = oberauer_lin_2017,
model = sdm(resp_error = "dev_rad")
)
cat(scode1)
#> // generated with brms 2.23.0 and bmm 1.4.2
#> functions {
#> /* compute the tan_half link
#> * Args:
#> * x: a scalar in (-pi, pi)
#> * Returns:
#> * a scalar in (-Inf, Inf)
#> */
#> real tan_half(real x) {
#> return tan(x / 2);
#> }
#> /* compute the tan_half link (vectorized)
#> * Args:
#> * x: a vector in (-pi, pi)
#> * Returns:
#> * a vector in (-Inf, Inf)
#> */
#> vector tan_half(vector x) {
#> return tan(x / 2);
#> }
#> /* compute the inverse of the tan_half link
#> * Args:
#> * y: a scalar in (-Inf, Inf)
#> * Returns:
#> * a scalar in (-pi, pi)
#> */
#> real inv_tan_half(real y) {
#> return 2 * atan(y);
#> }
#> /* compute the inverse of the tan_half link (vectorized)
#> * Args:
#> * y: a vector in (-Inf, Inf)
#> * Returns:
#> * a vector in (-pi, pi)
#> */
#> vector inv_tan_half(vector y) {
#> return 2 * atan(y);
#> }
#>
#>
#>
#> // log of the numerator of the sdm likelihood
#> real sdm_simple_lpdf(vector y, vector mu, vector c, vector kappa) {
#> int N = size(y);
#> vector[N] num = exp(fma(kappa,cos(y-mu)-1,c)) ;
#> real out = dot_product(num, sqrt(kappa));
#> out *= inv(sqrt2()) * inv_sqrt(pi());
#> return(out);
#> }
#>
#> // Minus the log of the mean of exp(eta(u)) over the circle, where
#> // eta(u) = exp(c) sqrt(kappa / (2 pi)) exp(kappa (cos(u) - 1)) and c is on the
#> // log scale. eta is even, so [0, pi] carries the whole integral. The spike at
#> // u = 0 narrows with kappa and its peak grows with c (#615). Its squared
#> // inverse width is about peak * kappa at the top, at most 2 kappa log(peak)
#> // where it meets the floor, and kappa for a small peak, so kappa (peak + 1)
#> // measures how sharp it is. Up to 4000, which takes in exp(c) <= 50 with
#> // kappa <= 30, a fixed rule with nodes crowded toward u = 0 resolves it;
#> // wide_cosm1 holds cos(u) - 1 at its nodes. Sharper spikes get Gauss-Legendre
#> // on panels cut at multiples of two widths: w, where eta falls to about 1,
#> // and s, the SD of exp(eta) near the peak, which is far narrower than w once
#> // the peak is large. The last w cut is at 5 w because past 3 w the spike of a
#> // small peak still carries up to 1.2e-6 of the integral above its floor,
#> // and 1.2e-13 past 5 w. eta - peak is formed as
#> // peak expm1(kappa (cos(u) - 1)), with cos(u) - 1 free of cancellation,
#> // because cos(u) - 1 and eta - peak both cancel near the peak.
#> real sdm_simple_ldenom(real c, real kappa, data vector gl_x, data vector gl_w,
#> data vector wide_cosm1, data vector wide_w) {
#> int n_nodes = num_elements(gl_x);
#> real log_peak = c + 0.5 * log(kappa / (2 * pi()));
#> real peak = exp(log_peak);
#> if (kappa * (peak + 1) <= 4000) {
#> return log(pi()) - peak - log(dot_product(wide_w, exp(peak * expm1(kappa * wide_cosm1))));
#> }
#> real w = fmin(pi(), sqrt(2 * fmax(log_peak, 1) / kappa));
#> real s = exp(-0.5 * (log_peak + log(kappa)));
#> vector[9] cuts = [0, pi(), 0.5 * w, w, 2 * w, 5 * w, s, 3 * s, 9 * s]';
#> real log_int = negative_infinity();
#> for (i in 1:9) {
#> cuts[i] = fmin(fmax(cuts[i], 0), pi());
#> }
#> cuts = sort_asc(cuts);
#> for (i in 1:8) {
#> if (cuts[i + 1] > cuts[i]) {
#> real half = (cuts[i + 1] - cuts[i]) / 2;
#> vector[n_nodes] u = (cuts[i] + cuts[i + 1]) / 2 + half * gl_x;
#> vector[n_nodes] d = peak * expm1(-2 * kappa * square(sin(u / 2)));
#> real d_max = max(d);
#> log_int = log_sum_exp(log_int, log(half) + d_max + log(dot_product(gl_w, exp(d - d_max))));
#> }
#> }
#> return log(pi()) - peak - log_int;
#> }
#>
#> real sdm_simple_run_ldenom(vector c, vector kappa, data vector gl_x,
#> data vector gl_w, data vector wide_cosm1,
#> data vector wide_w, int G_runs,
#> array[] int run_start, array[] int run_count) {
#> real out = 0;
#> for (g in 1:G_runs) {
#> int n = run_start[g];
#> out += run_count[g] * sdm_simple_ldenom(c[n], kappa[n], gl_x, gl_w, wide_cosm1, wide_w);
#> }
#> return(out);
#> }
#>
#> real sdm_simple_run_ldenom_slice(vector c, vector kappa, data vector gl_x,
#> data vector gl_w, data vector wide_cosm1,
#> data vector wide_w, int start, int end,
#> int G_runs, array[] int run_start,
#> array[] int run_count) {
#> real out = 0;
#> for (g in 1:G_runs) {
#> if (run_start[g] >= start && run_start[g] <= end) {
#> int n = run_start[g] - start + 1;
#> out += run_count[g] * sdm_simple_ldenom(c[n], kappa[n], gl_x, gl_w, wide_cosm1, wide_w);
#> }
#> }
#> return(out);
#> }
#> }
#> data {
#> int<lower=1> N; // total number of observations
#> vector[N] Y; // response variable
#> int prior_only; // should the likelihood be ignored?
#> vector[16] sdm_gl_x;
#> vector[16] sdm_gl_w;
#> vector[32] sdm_wide_cosm1;
#> vector[32] sdm_wide_w;
#> int G_sdm_runs;
#> array[G_sdm_runs] int sdm_run_start;
#> array[G_sdm_runs] int sdm_run_count;
#> }
#> transformed data {
#> // fail fast if the run metadata computed in R does not match the data Stan received
#> {
#> int sdm_run_total = 0;
#> for (g in 1:G_sdm_runs) {
#> if (sdm_run_start[g] < 1 || sdm_run_start[g] > N) {
#> reject("bmm error: sdm_run_start[", g, "] = ", sdm_run_start[g],
#> " is outside the data range [1, ", N, "]. The SDM run metadata ",
#> "does not match the model data. Please report this at ",
#> "https://github.com/popov-lab/bmm/issues");
#> }
#> if (g > 1 && sdm_run_start[g] <= sdm_run_start[g - 1]) {
#> reject("bmm error: sdm_run_start must be strictly increasing. The SDM ",
#> "run metadata does not match the model data. Please report ",
#> "this at https://github.com/popov-lab/bmm/issues");
#> }
#> sdm_run_total += sdm_run_count[g];
#> }
#> if (sdm_run_total != N) {
#> reject("bmm error: sum(sdm_run_count) = ", sdm_run_total, " but N = ", N,
#> ". The SDM run metadata does not match the model data. Please ",
#> "report this at https://github.com/popov-lab/bmm/issues");
#> }
#> }
#> }
#> parameters {
#> real Intercept_c; // temporary intercept for centered predictors
#> real Intercept_kappa; // temporary intercept for centered predictors
#> }
#> transformed parameters {
#> real Intercept; // temporary intercept for centered predictors
#> // prior contributions to the log posterior
#> real lprior = 0;
#> Intercept = 0;
#> lprior += student_t_lpdf(Intercept_c | 5, 2, 0.75);
#> lprior += student_t_lpdf(Intercept_kappa | 5, 1.75, 0.75);
#> }
#> model {
#> // likelihood including constants
#> if (!prior_only) {
#> // initialize linear predictor term
#> vector[N] mu = rep_vector(0.0, N);
#> // initialize linear predictor term
#> vector[N] c = rep_vector(0.0, N);
#> // initialize linear predictor term
#> vector[N] kappa = rep_vector(0.0, N);
#> mu += Intercept;
#> c += Intercept_c;
#> kappa += Intercept_kappa;
#> mu = inv_tan_half(mu);
#> kappa = exp(kappa);
#> target += sdm_simple_lpdf(Y | mu, c, kappa);
#> target += sdm_simple_run_ldenom(c, kappa, sdm_gl_x, sdm_gl_w, sdm_wide_cosm1,
#> sdm_wide_w, G_sdm_runs, sdm_run_start,
#> sdm_run_count);
#> target += -(log2()+log(pi()))*N;
#> }
#> // priors including constants
#> target += lprior;
#> }
#> generated quantities {
#> // actual population-level intercept
#> real b_Intercept = Intercept;
#> // actual population-level intercept
#> real b_c_Intercept = Intercept_c;
#> // actual population-level intercept
#> real b_kappa_Intercept = Intercept_kappa;
#> }