[Math] Generate random numbers with a modified PERT distribution

MATLABprobability distributionsrandom

I want to generate random numbers based on the modified PERT distribution. The modified PERT distribution is a special case of the beta distribution and is defined as:

$$f_X(x) = \frac{1}{B(\alpha_1,\alpha_2)}\frac{(x – a)^{\alpha_1 – 1}(b – x)^{\alpha_2 – 1}}{(b – a)^{\alpha_1 + \alpha_2 – 1}}$$

Where $B(\alpha_1,\alpha_2)$ is the beta function and

$$\alpha_1 = 1 + \gamma \left(\frac{m – a}{b – a}\right), \alpha_2 = 1 + \gamma \left( \frac{b – m}{b – a}\right)$$

Herein $a \leq m \leq b$. $[a,b]$ is the support of the probability densitify function (where it is none zero) and $m$ is the modal value and $\gamma$ is a shape parameter which influences the variance.

You can read more about the modified PERT distribution here: vosesoftware or mathematica documentation.
From what I understand from the last link the cumulative distribution function should be

$$F_X(x) = I_{\alpha_3}(\alpha_1,\alpha_2) = \frac{B(\alpha_3,\alpha_1,\alpha_2)}{B(\alpha_1,\alpha_2)}$$

With $I_{\alpha_3}(\alpha_1,\alpha_2)$ being the regularized incomplete beta function, $B(\alpha_3,\alpha_1,\alpha_2)$ the incomplete beta function. Herein

$$\alpha_3 = \frac{x – a}{b – a}$$

Now I want to use the inverse transform sampling to generate random numbers. This means that I generate a pseudorandom number $u$ which has a standard uniform distribution, rand in Matlab will then suffice. Then I need to compute a value $x$ for which $F_X(x) = u \Leftrightarrow F_X(u)^{-1} = x$. So I do need to find $F_X(x)^{-1}$.

Now my question is can anyone compute $F_X(x)^{-1}$ for me? Further I would appreciate a Matlab implementation, this is however not necessary since I can also try it for myself of course 🙂

— EDIT —
Apparently Matlab actually has a implementation of the inverse regularized incomplete beta function, it howevers calls it inverse incomplete beta function

Best Answer

If you have $0 \le u \le 1$ you get: $$u=F_X(x)=I_{\alpha_3}(\alpha_1,\alpha_2)$$ $$⟹\alpha_3=\alpha_3(u)=\mathrm{betaincinv}(u,\alpha_1,\alpha_2)$$ $$⟹x=F_X^{−1}(u)=(b−a)\alpha_3(u)+a$$

Related Question