How to fit a Beta distribution to three estimates from an "expert"?

Question

I'm modeling a process time, $X$, for a simulation study and have an "expert" estimate of the minimum, $\hat a$, the most likely (mode), $\hat m$, and the maximum, $\hat b$.

I'd prefer to avoid the possibly unrealistic use of the Triangular distribution, $X\sim \text{Tria}(a,m,b)$, with $a=\hat a$, $m = \hat m$, and $b = \hat b$.

Q: How do I use this information to fit a Generalized Beta distribution to this?

---

The Generalized Beta distribution has the form of $a+(b-a)\,\text{Beta}(\alpha_1,\alpha_2)$ with first shape parameter $\alpha_1$ and second shape parameter $\alpha_2$ where $\text{Beta}(\alpha_1,\alpha_2)$ is the Beta distribution on the $[0,1]$ interval.

---

Update:
Due to concerns this is unrelated to OR, some context may be warranted to confirm it is.

This situation arises frequently when conducting input modeling for a Simulation study and further data is unavailable or too expensive or troublesome to get. There are legitimate concerns over the Triangular distribution due to the "tails" and this approach provides a viable alternative that may better approximate reality in some applications.

Answer 1

Let $a = \hat a$ and $b = \hat b$.

Denote $r = \frac{\hat b - \hat m}{\hat m - \hat a}$.

Then choose the shape parameters to be

$\alpha_1 = \frac{4 + 3r + r^2}{1 + r^2}$

and

$\alpha_2 = \frac{1 + 3r + 4r^2}{1 + r^2}$.

This will produce a Generalized Beta distribution with mode $m = \hat m$ and a variance $\text{Var}(X) = \left(\frac{\hat b - \hat a}{6}\right)^2$. Then mean will be $\text{E}[X] = \mu = \frac{\alpha_1 b + \alpha_2 a}{\alpha_1 + \alpha_2}$. ^[1]

Notice there is no assumption about how the mean depends on $\hat a, \hat m, \hat b$ as would occur as in a common PERT assumption $\mu = \frac{a + 4m + b}{6}$.

---

If there are concerns the "expert" estimates of $\hat a$ and $\hat b$ are the smallest and largest (respectively) values they've observed, then they aren't the true bounds since $P(X = x)=0$ for continuous random variables. To address this concern, one can do a heuristic expansion.

For example, if the expansion coefficient is 20% ($\theta = 0.20$), and the initial estimate range was $\hat \rho = \hat b - \hat a$, then the new estimated bounds are

$a = \hat a - \frac{\theta \hat \rho - \hat\rho}{2}$

and

$b = \hat b + \frac{\theta \hat \rho - \hat\rho}{2}$.

---

Update: Someone rightly asked in the comments,

Why constrain the standard deviation to be one-sixth of the range? The problem is under-determined, so there has to be some such constraint if it's to be solved at all, but why that one in particular?

This is a good question. Since this Generalized Beta distribution has four parameters $(a,b,\alpha_1,\alpha_2)$ and you are starting with $(\hat a, \hat m, \hat b)$, there are many, many (even infinitely many for our mathy friends) ways to get this done. I believe the resulting variance $\text{Var}(X) = \left(\frac{\hat b - \hat a}{6}\right)^2$ is baked into the results by way of the choices for $\alpha_1$ and $\alpha_2$.

In addition to being an improvement over traditional PERT assumptions, this approach is simple, easy to understand, and works well when there is limited data.

I'll look into further reasoning for this constraint and hope to add some of the reasoning here. Until then, I refer interested parties to Kuhl et al. (2010).^[1]

Code Demonstrating the Formulas Are Correct:

% MATLAB R2018b
ahat = 4;
bhat = 24;
mhat = 12;

a = ahat;
b = bhat;
r = (bhat - mhat)./(mhat - ahat);
alpha1 = (4 + 3.*r + r.^2)./(1 + r.^2);
alpha2 = (1 + 3.*r + 4*r.^2)./(1 + r.^2);

% Simulation
NumSamples = 50000;
X = a + (b-a).*betarnd(alpha1,alpha2,NumSamples,1);
SimMean = mean(X)
SimVar = var(X)

% Analytical 
AnalyticalMean = (alpha1.*b + alpha2.*a)./(alpha1 + alpha2)
AnalyticalVar = ((b-a)./6).^2

---

_{Reference:
[1] Kuhl, M., Ivy, J., Lada, E., Steiger, N., Wagner, M., & Wilson, J. (2010). Univariate inputmodels for stochastic simulation. Journal of Simulation, Vol 4, No 2, 81–97.}