So I'm considering a simple linear regression model with $p = 1$ predictor $$y = \beta x + \epsilon$$ where $\epsilon \sim N(0,\sigma^2)$. I want to use a generalised beta distribution of the second kind(GB2) as a prior for $\beta$ and an inverse gamma distribution for $\sigma^2$, and then sample from the posterior distribution. Would this be doable using Stan? If so, could you please show me how? I noticed GB2 is not a built in.

Please keep in mind I am completely new to using Stan.


1 Answer 1


In Stan, you can declare that a variable has a particular probability distribution, by either using a "sampling statement", e.g.

x ~ normal(mu, sigma);

or equivalently incrementing the built-in variable target by the log-likelihood of the distribution (see the Stan reference manual for more information):

target += normal_lpdf(x | mu, sigma);

The latter form is more flexible, as you can use any custom distribution, even if it is not built-in, by incrementing the target variable with its log-likelihood.

In your case, if we use the definition of the GB2 distribution from Wikipedia, the probability density function is (assuming non-negative $a$ for simplicity):

$$ f(x; a, b, p, q) = \frac{ax^{ap-1}}{b^{ap} B(p, q) \left(1 + \left(\frac{x}{b} \right)^a \right)^{p+q}}. $$

The log-likelihood is then:

$$ L(x; a, b, p, q) = \log a + (ap - 1) \log x - ap \log b - \log B(p, q) - (p + q) \log \left( 1 + \left ( \frac{x}{b} \right)^a \right). $$

So we just need to increment the target variable by this quantity.

An implementation of the complete model could be the following.

// model.stan

data {
  int<lower=0> n;         // number of samples
  vector[n] x;            // x values
  vector[n] y;            // y values

parameters {
  real slope;             // the slope of the regression line
  real<lower=0> sigma2;   // the squared standard deviation

model {
  real a = 1;             // the parameters of the GB2 prior
  real b = 2;
  real p = 2;
  real q = 1;
  // the GB2 prior as an increment of the log-likelihood
  target += log(a) + (a * p - 1) * log(slope) - a * p * log(b) - lbeta(p, q) - (p + q) * log (1 + (slope / b)^a);
  // the inverse Gamma prior for sigma2 using a sampling statement
  sigma2 ~ inv_gamma(2, 1);
  // and finally, the regression model
  y ~ normal(slope * x, sqrt(sigma2));

Note that the lbeta function is the logarithm of the $B(.,.)$ function. You can also wrap the GB2 log-likelihood as a function to easily reuse it if needed, as is done here (although I am not sure if the expression used in the link produces the same output as the expression used here).

An example usage of the model in R with simulated data is the following.

# main.R



slope <- 1.5        # the parameters we wish to estimate
sigma2 <- 1

# simulated data
n <- 20
x <- runif(n, 0, 5)
y <- rnorm(n, slope * x, sqrt(sigma2))

mdl <- stan("model.stan", data = list(n = n, x = x, y = y))

# ...
#         mean se_mean   sd ...
# slope   1.47    0.00 0.07 ...
# sigma2  0.64    0.00 0.20 ...
# lp__   -8.01    0.03 1.05 ...
# ...

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.