Four parameter beta distribution

Modeling
1

I am trying to reproduce a model in a paper I am reviewing using Stan. The author(s) used MCMC and specified a four-parameter beta distribution as prior for some of the parameters. When I searched for it, I couldn’t find anything close to it in the Stan manual. Then, I found this page describing this distribution.

https://en.wikipedia.org/wiki/Beta_distribution#Four_parameters

Is there a way to specify this prior on a parameter in Stan?

Thank you.

2

According to the wikipedia page:

The probability density function of the four parameter beta distribution is equal to the two parameter distribution, scaled by the range ( c - a ), (so that the total area under the density curve equals a probability of one), and with the “y” variable shifted and scaled

Denoting the PDF of the 4-parameter Beta (with lower bound l and upper bound u):

f(y;\alpha,\beta,l,u)

And the 2-parameter Beta:

f(x;\alpha,\beta)

Where:

x = \frac{y-l}{u-l}

The relationship between the 2-parameter and 4-parameter Beta distribution PDFs is given by:

f(y;\alpha,\beta,l,u) = \frac{f(x;\alpha,\beta)}{u-l}

Then because Stan needs the log-PDF, the distribution we want to express in Stan is:

\log\left(f(y;\alpha,\beta,l,c)\right) = \log\left(f(x;\alpha,\beta)\right) - \log(u-l)

We can then define this as a function in Stan:

functions {
  real beta_4p_lpdf(real y, real alpha, real beta, real lower, real upper) {
    // Scale 4-parameter Beta RV to 2-parameter Beta RV
    real x = (y - lower) / (upper - lower);
    
    // Return scaled 2-parameter beta lpdf
    return beta_lpdf(x | alpha, beta) - log(upper - lower);
  }
}

As a minimal example, you would then use like so:

functions {
  real beta_4p_lpdf(real y, real alpha, real beta, real lower, real upper) {
    // Scale 4-parameter Beta RV to 2-parameter Beta RV
    real x = (y - lower) / (upper - lower);
    
    // Return scaled 2-parameter beta lpdf
    return beta_lpdf(x | alpha, beta) - log(upper - lower);
  }
}

data {
  int N;
  vector[N] y;
  real lower;
  real upper;
}

parameters {
  real alpha;
  real beta;
}

model {
  for(n in 1:N) {
    y[n] ~ beta_4p(alpha, beta, lower, upper);
  }
}

Note that I haven’t tested any of this code, so I would recommend spot-checking some values against another implementation of the 4-parameter beta lpdf to make sure that all is correct

EDIT: Updated the specification for x, since I had accidentally used subtraction instead of division!

4 Likes
3

This is so cool! Thank you so much.

1 Like