Generate Stan code for bmm models

Given the model, the data and the formula for the model, this function will return the combined stan code generated by bmm and brms

Usage

# S3 method for class 'bmmformula'
stancode(object, data, model, prior = NULL, ...)

Arguments

object: A bmmformula object
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 bmmodel object for required argurments.
model: A description of the model to be fitted. This is a call to a bmmodel such as mixture3p() function. Every model function has a number of required arguments which need to be specified within the function call. Call supported_models() to see the list of supported models and their required arguments
prior: One or more brmsprior objects created by brms::set_prior() or related functions and combined using the c method or the + operator. See also default_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 of brms::stancode() for more details

Value

A character string containing the fully commented Stan code to fit a bmm model.

Examples

scode1 <- stancode(bmf(c ~ 1, kappa ~ 1),
  data = oberauer_lin_2017,
  model = sdm(resp_error = "dev_rad")
)
cat(scode1)
#> // generated with brms 2.22.0 and bmm 1.0.1.9000
#> 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);
#>    }
#> 
#> 
#> 
#>   // utility function trick for converting real to integer type
#>   int bin_search(real x, int min_val, int max_val) {
#>       int mid_p;
#>       int L = min_val;
#>       int R = max_val;
#>       while(L < R) {
#>         mid_p = (R-L) %/% 2;
#>         if (L + mid_p < x) {
#>           L += mid_p + 1;
#>         } else if (L + mid_p > x) {
#>           R = L + mid_p - 1;
#>         } else {
#>           return(L + mid_p);
#>         }
#>       }
#>       return(L);
#>     }
#> 
#>   // utility function for determining optimal number of chebyshev points for the denominator approximation
#>   int get_m(real c, real kappa) {
#>     real m = floor(2 * exp(0.4*c) * kappa^(fma(c,0.0145,0.7)) + 0.5)+2;
#>     int M = bin_search(m, 2, 200);
#>     return(M);
#>   }
#> 
#>   // 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);
#>   }
#> 
#>   // log of the normalization constant, approximated by chebyshev quadrature
#>   real sdm_simple_ldenom_chquad_adaptive(real c, real kappa, matrix CN) {
#>     int m = get_m(c,kappa);
#>     vector[m] cosn = CN[1:m,m];
#>     vector[m] fn = exp(fma(kappa,cosn,c)) * sqrt(kappa) * inv(sqrt2()) * inv_sqrt(pi());
#>     real out = -log_sum_exp(fn)+log(m);
#>     return(out);
#>   }
#> }
#> data {
#>   int<lower=1> N;  // total number of observations
#>   vector[N] Y;  // response variable
#>   int prior_only;  // should the likelihood be ignored?
#> }
#> transformed data {
#>     // precompute chebyshev points
#>   matrix[200,200] COSN;
#>   for (m in 1:200) {
#>     for (i in 1:m) {
#>       COSN[i,m] = cos((2*i-1)*pi()/(2*m))-1;
#>     }
#>   }
#> }
#> 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
#>   real lprior = 0;  // prior contributions to the log posterior
#>   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);
#>       // adaptive calculation of the normalization constant
#>     real z;
#>     for (n in 1:N) {
#>     	if (n == 1 || c[n] != c[n-1] || kappa[n] != kappa[n-1]) {
#>     		z = sdm_simple_ldenom_chquad_adaptive(c[n],kappa[n],COSN);
#>     	}
#>     	target += z;
#>     }
#>     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;
#> }

Usage

Arguments

Value

See also

Examples