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;
  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 $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;
   111:    real sigma_adj = variance_adjusted_sgt(sigma, lambda, p, q);
                                                                 ^
   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;
   111:    real sigma_adj = variance_adjusted_sgt(sigma, lambda, p, q);
                                                                 ^
   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;
    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, 1, q);
    
    if (is_inf(q)) 
      return uniform_lpdf(x | mu - sigma_adj, mu + sigma_adj);
    
    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 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)]);
  }
}
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;
    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, 1, q);
    
    if (is_inf(q)) 
      return uniform_lpdf(x | mu - sigma_adj, mu + sigma_adj);
    
    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 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)]);
  }


'

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;
    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);
    
    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 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)]);
  }
}
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;
    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);
    
    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 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)]);
  }
'

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