I am trying to define the skew_generalized_t distribution as a custom family in brms, but don’t fully understand what I am doing. I tried using @spinkney’s function from here: spinkney / helpful_stan_functions and following the brms custom family vignette
stan_funs <- '
real variance_adjusted_sgt(real sigma, real lambda, real p, real q) {
if (p * q <= 2)
reject("p * q must be > 2 found p * q = ", p * q);
if (is_inf(q))
return sigma
* inv_sqrt((pi() * (1 + 3 * lambda ^ 2) * tgamma(3.0 / p)
- 16 ^ (1.0 / p) * lambda ^ 2 * (tgamma(1.0 / 2 + 1.0 / p)) ^ 2
* tgamma(1.0 / p))
/ (pi() * tgamma(1.0 / p)));
return sigma
/ (q ^ (1.0 / p)
* sqrt((3 * lambda ^ 2 + 1) * (beta(3.0 / p, q - 2.0 / p) / beta(1.0 / p, q))
- 4 * lambda ^ 2 * (beta(2.0 / p, q - 1.0 / p) / beta(1.0 / p, q)) ^ 2));
}
vector mean_centered_sgt(vector x, real sigma, real lambda, real p, real q) {
if (p * q <= 1)
reject("p * q must be > 1 found p * q = ", p * q);
if (is_inf(q))
return x + (2 ^ (2.0 / p) * sigma * lambda * tgamma(1.0 / 2 + 1.0 / p)) / sqrt(pi());
return x + (2 * sigma * lambda * q ^ (1.0 / p) * beta(2 / p, q - 1.0 / p)) / beta(1.0 / p, q);
}
real mean_centered_sgt(real x, real sigma, real lambda, real p, real q) {
if (p * q <= 1)
reject("p * q must be > 1 found p * q = ", p * q);
if (is_inf(q))
return x + (2 ^ (2.0 / p) * sigma * lambda * tgamma(1.0 / 2 + 1.0 / p)) / sqrt(pi());
return x + (2 * sigma * lambda * q ^ (1.0 / p) * beta(2 / p, q - 1.0 / p)) / beta(1.0 / p, q);
}
real mean_centered_sgt(real x, real sigma, real lambda, real q) {
if (q <= 1)
reject("q must be > 1 found q = ", q);
if (is_inf(q))
return x + (4 * sigma * lambda * tgamma(1.0 / 2 + 1.0)) / sqrt(pi());
return x + (2 * sigma * lambda * q * beta(2, q - 1.0)) / beta(1.0, q);
}
real skew_generalized_t_lpdf(vector x, real mu, real sigma, real lambda, real p, real q) {
if (sigma <= 0)
reject("sigma must be > 0 found sigma = ", sigma);
if (lambda >= 1 || lambda <= -1)
reject("lambda must be between (-1, 1) found lambda = ", lambda);
if (p <= 0)
reject("p must be > 0 found p = ", p);
if (q <= 0)
reject("q must be > 0 found q = ", q);
int N = num_elements(x);
real out = 0;
real sigma_adj = variance_adjusted_sgt(sigma, lambda, p, q);
if (is_inf(q) && is_inf(p))
return uniform_lpdf(x | mu - sigma_adj, mu + sigma_adj);
vector[N] r = mean_centered_sgt(x, sigma_adj, lambda, p, q) - mu;
vector[N] s;
for (n in 1:N) s[n] = r[n] < 0 ? -1 : 1;
if (is_inf(q) && !is_inf(p)) {
out = sum( ( abs(r) ./ (sigma_adj * (1 + lambda * s))) ^ p );
return log(p) - log(2) - log(sigma_adj) - lgamma(1.0 / p) - out;
} else {
out = sum(log1p( abs(r) ^ p ./ (q * sigma_adj ^ p * pow(1 + lambda * s, p))));
}
return N * (log(p) - log2() - log(sigma_adj) - log(q) / p - lbeta(1.0 / p, q)) - (1.0 / p + q) * out;
}
real skew_t_lpdf(vector x, real mu, real sigma, real lambda, real q) {
if (sigma <= 0)
reject("sigma must be > 0 found sigma = ", sigma);
if (lambda >= 1 || lambda <= -1)
reject("lambda must be between (-1, 1) found lambda = ", lambda);
if (q <= 0)
reject("q must be > 0 found q = ", q);
int N = num_elements(x);
real out = 0;
real sigma_adj = variance_adjusted_sgt(sigma, lambda, p, q);
if (is_inf(q) && is_inf(p))
return uniform_lpdf(x | mu - sigma_adj, mu + sigma_adj);
vector[N] r = mean_centered_sgt(x, sigma_adj, lambda, q) - mu;
vector[N] s;
for (n in 1:N) s[n] = r[n] < 0 ? -1 : 1;
if (is_inf(q)) {
out = sum( ( abs(r) ./ (sigma_adj * (1 + lambda * s))) );
return -log(2) - log(sigma_adj) - out;
} else {
out = sum(log1p( abs(r) ./ (q * sigma_adj * (1 + lambda * s)));
}
return N * (-log2() - log(sigma_adj) - log(q) - lbeta(1.0, q)) - (1.0 + q) * out;
}
real skew_generalized_t_lcdf(real x, real mu, real sigma, real lambda, real p, real q) {
if (sigma <= 0)
reject("sigma must be > 0 found sigma = ", sigma);
if (lambda >= 1 || lambda <= -1)
reject("lambda must be between (-1, 1) found lambda = ", sigma);
if (p <= 0)
reject("p must be > 0 found p = ", p);
if (q <= 0)
reject("q must be > 0 found q = ", q);
if (is_inf(x) && x < 0)
return 0;
if (is_inf(x) && x > 0)
return 1;
real sigma_adj = variance_adjusted_sgt(sigma, lambda, p, q);
real x_cent = mean_centered_sgt(x, sigma_adj, lambda, p, q);
real r = x_cent - mu;
real lambda_new;
real r_new;
if (r > 0) {
lambda_new = -lambda;
r_new = -r;
} else {
lambda_new = lambda;
r_new = r;
}
if (!is_inf(p) && is_inf(q) && !is_inf(x))
return log_sum_exp([
log1m(lambda_new) + log2(),
log(lambda_new - 1) + log2() + beta_lcdf((r_new / (sigma * (1 + lambda_new))) ^ p | p, 1)]);
if (is_inf(p) && !is_inf(x))
return uniform_lcdf(x | mu, sigma);
if (is_inf(x) && x < 0)
return 0;
if (is_inf(x) && x > 0)
return 1;
return log_sum_exp([
log1m(lambda_new) + log2(),
log(lambda_new - 1) + log2() + beta_lcdf(1.0 / (1 + q * (sigma * (1 - lambda_new) / (-r_new)) ^ p) | 1.0 / p, q)
]);
}
'
variance_adjusted_sgt <- custom_family(
"variance_adjusted_sgt ", dpars = c("mu", "sigma","l", "p", "q"),
links = c("identity", "log", "identity", "identity", "identity"),
#lb = c(NA, NA), ub = c(NA, NA),
type = "int", vars = "vint1[n]"
)
stanvars <- stanvar(scode = stan_funs, block = "functions")
SPT_distance_model <- brm(bf(SPT ~ distance + (distance | tag)),
data = sleep_distance[complete.cases(sleep_distance[,c("distance", "SPT")]),],
save_pars = save_pars(all = TRUE),
iter = 20000,
init = 0,
family = variance_adjusted_sgt, stanvars = stanvars,
backend = "cmdstanr")
But I get the following error when I try to run the actual model
Syntax error in ‘C:\Users\salavi\AppData\Local\Temp\RtmpiYm3F3/model_6fb0009136d23a52ef244b78652b31ee.stan’, line 125, column 66 to column 67, parsing error:
123: return -log(2) - log(sigma_adj) - out;
124: } else {
125: out = sum(log1p( abs(r) ./ (q * sigma_adj * (1 + lambda * s)));
^
126: }
127:Ill-formed function application. Expect comma-separated list of expressions followed by “)” after “(”.
Error in
processx::run(command = stanc_cmd, args = c(stan_file, stanc_flags), …:
! System command ‘stanc.exe’ failed
Exit status: 1
Stderr (last 10 lines, see$stderrfor more):
123: return -log(2) - log(sigma_adj) - out;
124: } else {
125: out = sum(log1p( abs(r) ./ (q * sigma_adj * (1 + lambda * s)));
^
126: }
127:Ill-formed function application. Expect comma-separated list of expressions followed by “)” after “(”.
Type .Last.error to see the more details.
Maybe @spinkney or @paul.buerkner have some insights? Any advice as to what I am doing wrong will be much appreciated!!
