# Using a Generalized Beta Distribution of the Second Kind as a Prior in Stan Linear Regression

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.

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

library(rstan)

set.seed(100)

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))
print(mdl)

# ...
#         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 ...
# ...