// generated with brms 2.11.1 functions { /* compute a latent Gaussian process * Args: * x: array of continuous predictor values * sdgp: marginal SD parameter * lscale: length-scale parameter * zgp: vector of independent standard normal variables * Returns: * a vector to be added to the linear predictor */ vector gp(vector[] x, real sdgp, vector lscale, vector zgp) { int Dls = rows(lscale); int N = size(x); matrix[N, N] cov; if (Dls == 1) { // one dimensional or isotropic GP cov = cov_exp_quad(x, sdgp, lscale[1]); } else { // multi-dimensional non-isotropic GP cov = cov_exp_quad(x[, 1], sdgp, lscale[1]); for (d in 2:Dls) { cov = cov .* cov_exp_quad(x[, d], 1, lscale[d]); } } for (n in 1:N) { // deal with numerical non-positive-definiteness cov[n, n] += 1e-12; } return cholesky_decompose(cov) * zgp; } /* Spectral density function of a Gaussian process * Args: * x: array of numeric values of dimension NB x D * sdgp: marginal SD parameter * lscale: vector of length-scale parameters * Returns: * numeric values of the function evaluated at 'x' */ vector spd_cov_exp_quad(vector[] x, real sdgp, vector lscale) { int NB = dims(x)[1]; int D = dims(x)[2]; int Dls = rows(lscale); vector[NB] out; if (Dls == 1) { // one dimensional or isotropic GP real constant = square(sdgp) * (sqrt(2 * pi()) * lscale[1])^D; real neg_half_lscale2 = -0.5 * square(lscale[1]); for (m in 1:NB) { out[m] = constant * exp(neg_half_lscale2 * dot_self(x[m])); } } else { // multi-dimensional non-isotropic GP real constant = square(sdgp) * sqrt(2 * pi())^D * prod(lscale); vector[Dls] neg_half_lscale2 = -0.5 * square(lscale); for (m in 1:NB) { out[m] = constant * exp(dot_product(neg_half_lscale2, square(x[m]))); } } return out; } /* compute an approximate latent Gaussian process * Args: * X: Matrix of Laplacian eigen functions at the covariate values * sdgp: marginal SD parameter * lscale: vector of length-scale parameters * zgp: vector of independent standard normal variables * slambda: square root of the Laplacian eigen values * Returns: * a vector to be added to the linear predictor */ vector gpa_coefficient(int NBgp_1, real sdgp, vector lscale, vector zgp, vector[] slambda) { vector[NBgp_1] diag_spd = sqrt(spd_cov_exp_quad(slambda, sdgp, lscale)); return diag_spd .* zgp; } } data { int N; // number of observations int Y[N]; // response variable // data related to GPs // number of sub-GPs (equal to 1 unless 'by' was used) int Kgp_1; int Dgp_1; // GP dimension // number of basis functions of an approximate GP int NBgp_1; // approximate GP basis matrices matrix[N, NBgp_1] Xgp_1; // approximate GP eigenvalues vector[Dgp_1] slambda_1[NBgp_1]; // sampling fractions int N_xi; matrix[N_xi,2] shape; int xi_id[N,2]; } transformed data { } parameters { // temporary intercept for centered predictors real Intercept; // GP standard deviation parameters vector[Kgp_1] sdgp_1; // GP length-scale parameters vector[1] lscale_1[Kgp_1]; // latent variables of the GP vector[NBgp_1] zgp_1; vector[N_xi] xi; } transformed parameters { vector[N] lxi_pair; vector[NBgp_1] coefficient = gpa_coefficient(NBgp_1, sdgp_1[1], lscale_1[1], zgp_1, slambda_1); for (n in 1:N){ lxi_pair[n] = log(xi[xi_id[n,1]]) + log(xi[xi_id[n,2]]); } } model { // initialize linear predictor term // priors including all constants target += normal_lpdf(Intercept | 0, 10); target += normal_lpdf(sdgp_1 | 0, 10); target += normal_lpdf(zgp_1 | 0, 1); target += inv_gamma_lpdf(lscale_1[1] | 2.22121,7.04478); for (i in 1:N_xi){ target += beta_lpdf(xi[i]|shape[i,1],shape[i,2]); } // likelihood including all constants target += poisson_log_glm_lpmf(Y | append_col(Xgp_1,lxi_pair), Intercept, append_row(coefficient,1)); }