# Generate random numbers with a modified PERT distribution

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

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$$