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!