# Brms custom skew_generalized_t family

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;

if (is_inf(q) && is_inf(p))

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;

if (is_inf(q) && is_inf(p))

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 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 `\$stderr` for 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!!

1 Like

Here are a subset of the data I am trying to model.
SLD_num.csv (38.7 KB)

If anyone has any advice as to how to define this model family I will be very appreciative

Hi @sea83, I will take a look next week when Im back from vacation.

Looking at it right here, if you addan end parenthesis at the end of that line 125, it may work. I added the skew_t (not generalized) recently and pushed without testing.

2 Likes

Thanks for being willing to take a look at this! Adding the closing parenthesis definitely made a difference, but now it returns different errors depending on if I use cmdstanr or not.

Using cmdstanr:

``````Semantic error in 'C:\Users\salavi\AppData\Local\Temp\RtmpAVEN7r/model_e2f48b6e84298852def69cb159062099.stan', line 111, column 56 to column 57:
-------------------------------------------------
109:    int N = num_elements(x);
110:    real out = 0;
^
112:
113:    if (is_inf(q) && is_inf(p))
-------------------------------------------------

Identifier 'p' not in scope.

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 `\$stderr` for more):
-------------------------------------------------
109:    int N = num_elements(x);
110:    real out = 0;
^
112:
113:    if (is_inf(q) && is_inf(p))
-------------------------------------------------

Identifier 'p' not in scope.
---
Type .Last.error to see the more details.
``````

Without cmdstanr:

``````Compiling Stan program...
Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
0

Semantic error in 'string', line 41, column 5 to column 22:

Identifier 'mean_centered_sgt' is already in use.
``````

Apologies if the problem is an obvious one that I should understand but don’t!

See if this works

``````functions {
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;

if (is_inf(q) && is_inf(p))

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;

if (is_inf(q))

vector[N] r = mean_centered_sgt(x, sigma_adj, lambda, 1, 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 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)]);
}
}
``````
1 Like

Attempting to follow the brms vignette

providing the stan functions

``````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;

if (is_inf(q) && is_inf(p))

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;

if (is_inf(q))

vector[N] r = mean_centered_sgt(x, sigma_adj, lambda, 1, 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 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)]);
}

'
``````

defining the custom family

``````variance_adjusted_sgt <- custom_family( name =  "variance_adjusted_sgt",
dpars = c("mu", "sigma","l", "p", "q"),
links = c("identity", "log", "identity", "identity", "identity"),
type = "int")

stanvars <- stanvar(scode = stan_funs, block = "functions")
``````

The model

``````sleep_distance\$SPT=as.integer(sleep_distance\$SPT)
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 = 2000,
init = 0,
family = variance_adjusted_sgt, stanvars = stanvars)
``````

returns

``````Compiling Stan program...
Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
0

Semantic error in 'string', line 42, column 7 to column 24:

Identifier 'mean_centered_sgt' is already in use.
``````

Using cmdstanr returns

``````Semantic error in 'C:\Users\salavi\AppData\Local\Temp\RtmpAVEN7r/model_950566e9f23f6cec0af5e44fba6fbcf7.stan', line 245, column 16 to column 72:
-------------------------------------------------
243:      }
244:      for (n in 1:N) {
245:        target += variance_adjusted_sgt_lpmf(Y[n] | mu[n], sigma, l, p, q);
^
246:      }
247:    }
-------------------------------------------------

A returning function was expected but an undeclared identifier 'variance_adjusted_sgt_lpmf' was supplied.

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 `\$stderr` for more):
-------------------------------------------------
243:      }
244:      for (n in 1:N) {
245:        target += variance_adjusted_sgt_lpmf(Y[n] | mu[n], sigma, l, p, q);
^
246:      }
247:    }
-------------------------------------------------

A returning function was expected but an undeclared identifier 'variance_adjusted_sgt_lpmf' was supplied.
---
Type .Last.error to see the more details
``````

What version of cmdstan are you using?

R version 4.2.1
RStudio 2022.07.1 Build 554
brms 2.17.5
cmdstanr version 0.5.2
CmdStan version: 2.30.1

Ok, that is the most recent version. The brms stuff I don’t know, it’s putting `_lpdf` on functions that aren’t lpdf. You’ll need a brms expert to weigh in there @paul.buerkner or @Solomon.

The name of the family-related custom Stan function needs to end on lpdf for continuous responses and on lpmf for discrete responses. Does that already answer your question? I am not sure which of the original question remained still.

This is an important likelihood to add for `brms` for sure.

But still an issue to get it to work.

Firstly, the correct custom family code is of course:

``````skew_generalized_t <- custom_family("skew_generalized_t",
dpars = c("mu", "sigma","lambda", "p", "q"),
links = c("identity", "log", "identity", "identity", "identity"),
type = "real")
``````

But the stan density function provided by @spinkney above and then @sea83 for `brms`, now returns the following error that for the life of me I cannot figure out how to fix.

The error message seems very clear but how to correct the density function? As far as I can determine, a vector is specified at the lpdf, but not so according to the error message below.

Tagging @spinkney and @paul.buerkner for thoughts on the fix.

Much appreciated to all three for raising this important likelihood that is especially useful for modelling nonlinear growth functions for various marine wildlife.

``````  -------------------------------------------------
225:      }
226:      for (n in 1:N) {
227:        target += skew_generalized_t_lpdf(Y[n] | mu[n], sigma, lambda, p, q);
^
228:      }
229:    }
-------------------------------------------------

Ill-typed arguments supplied to function 'skew_generalized_t_lpdf':
(real, real, real, real, real, real)
Available signatures:
(vector, real, real, real, real, real) => real
The first argument must be vector but got real

``````

The problem is that brms assumes a non-vectorized version of the lpdf while the function above is vectorized in y (first function argument is a vector not a real). Adjusting the lpdf to work with y being a single real (i.e., a non-vectorized implementation) will fix this.

1 Like

I think it’s just

``````functions {
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(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 = ", 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 = 1;
real out = 0;

if (is_inf(q) && is_inf(p))

real r = mean_centered_sgt(x, sigma_adj, lambda, p, q) - mu;
real s = r < 0 ? -1 : 1;

if (is_inf(q) && !is_inf(p)) {
out = (abs(r) ./ (sigma_adj * (1 + lambda * s))) ^ p;
return log(p) - log(2) - log(sigma_adj) - lgamma(1.0 / p) - out;
} else {
out = 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_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 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)]);
}
}
``````
1 Like

@spinkney many thanks as your revised density function works with `brms`.

I found that special attention to initial values was necessary to support adequate sampling. Perhaps there are other parameterisations of the function that could be explored.

Anyway, again, many thanks.

@paul.buerkner Thanks and apologies as I had completely forgotten that a non-vectorized density function was appropriate.

Many thanks @paul.buerkner @spinkney and @MilaniC! Things all seem to be working well thanks to your help.

Great! Can you post the final code that runs now? I am sure a lot of other users would benefit from this as well.

1 Like

Of course, I should have thought to do so right away.

Providing the updated stan functions

``````library(brms)
library(cmdstanr)

options(mc.cores = parallel::detectCores())

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(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 = ", 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 = 1;
real out = 0;

if (is_inf(q) && is_inf(p))

real r = mean_centered_sgt(x, sigma_adj, lambda, p, q) - mu;
real s = r < 0 ? -1 : 1;

if (is_inf(q) && !is_inf(p)) {
out = (abs(r) ./ (sigma_adj * (1 + lambda * s))) ^ p;
return log(p) - log(2) - log(sigma_adj) - lgamma(1.0 / p) - out;
} else {
out = 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_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 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)]);
}
'
``````

Defining the custom brms family

``````skew_generalized_t <- custom_family("skew_generalized_t",
dpars = c("mu", "sigma","lambda", "p", "q"),
links = c("identity", "log", "identity", "identity", "identity"),
type = "real")

stanvars <- stanvar(scode = stan_funs, block = "functions")
``````

The model

``````SPT_distance_model <- brm(bf(SPT ~ distance + (distance | tag)),
data = sleep_distance,
save_pars = save_pars(all = TRUE),
iter = 2000,
family = skew_generalized_t, stanvars = stanvars,
backend = "cmdstanr")
``````
2 Likes