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

