Let me first introduce some context.

The probability density function and cumulative distribution function of a Beta random variable with parameters $\alpha>0$, $\beta>0$ and support defined on $x\in[0,1]$ are given by: \begin{align} f(x)&=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{\text{B}(\alpha,\beta)} \end{align} and \begin{align} F(x)&=\frac{\text{B}(x;\alpha,\beta)}{\text{B}(\alpha,\beta)} \end{align} respectively, where $\text{B}(\alpha,\beta)$ is the Beta function defined as: \begin{align} \text{B}(\alpha,\beta)&=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}=\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}\,dt=\int_{0}^{\infty}\frac{y^{\alpha-1}}{(1+y)^{\alpha+\beta}}\,dy \end{align} and $\Gamma(x)$ is the Gamma function defined as: \begin{align} \Gamma(x)&=\int_{0}^{\infty}t^{s-1}e^{-t}\,dt \end{align} and $\text{B}(x;\alpha,\beta)$ is the incomplete Beta function defined as: \begin{align} \text{B}(x;\alpha,\beta)&=\int_{0}^{x}t^{\alpha-1}(1-t)^{\beta-1}\,dt=\int_{0}^{x/(1+x)}\frac{y^{\alpha-1}}{(1+y)^{\alpha+\beta}}\,dy \end{align}

The expectation of a Beta random variable that is right-truncated at $x=x_{t}$ is given by: \begin{align} \mathbb{E}[X\,|\,x\le x_{t}]&=\frac{1}{F(x_{t})}\int_{0}^{x_{t}}x\cdot f(x)\,dx\\ &=\frac{1}{F(x_{t})}\int_{0}^{x_{t}}x\cdot\frac{x^{\alpha-1}(1-x)^{\beta-1}}{\text{B}(\alpha,\beta)}\,dx\notag\\ &=\frac{1}{F(x_{t})}\int_{0}^{x_{t}}\frac{x^{\alpha}(1-x)^{\beta-1}}{\text{B}(\alpha,\beta)}\,dx\notag\\ &=\frac{1}{F(x_{t})}\frac{\text{B}(x_{t};\alpha+1,\beta)}{\text{B}(\alpha,\beta)}\notag \end{align}

The second moment of the right-truncated Beta distribution is given by: \begin{align} \mathbb{E}[X^{2}\,|\,x\le x_{t}]&=\frac{1}{F(x_{t})}\int_{0}^{x_{t}}x^{2}\cdot f(x)\,dx\\ &=\frac{1}{F(x_{t})}\int_{0}^{x_{t}}x^{2}\cdot\frac{x^{\alpha-1}(1-x)^{\beta-1}}{\text{B}(\alpha,\beta)}\,dx\notag\\ &=\frac{1}{F(x_{t})}\int_{0}^{x_{t}}\frac{x^{\alpha+1}(1-x)^{\beta-1}}{\text{B}(\alpha,\beta)}\,dx\notag\\ &=\frac{1}{F(x_{t})}\frac{\text{B}(x_{t};\alpha+2,\beta)}{\text{B}(\alpha,\beta)}\notag \end{align}

Finally, the variance of the right-truncated Beta distribution is given by: \begin{align} \text{Var}(X\,|\,x\le x_{t})&=\mathbb{E}[X^{2}\,|\,x\le x_{t}]-\Big(\mathbb{E}[X\,|\,x\le x_{t}]\Big)^{2}\\ &=\frac{1}{F(x_{t})}\frac{\text{B}(x_{t};\alpha+2,\beta)}{\text{B}(\alpha,\beta)}-\bigg(\frac{1}{F(x_{t})}\frac{\text{B}(x_{t};\alpha+1,\beta)}{\text{B}(\alpha,\beta)}\bigg)^{2} \end{align}

Thus, given you knew $F(x_{t})$ (the amount of probability in the tail you haven't observed), you could solve for the parameters of the system $(\hat{\alpha},\hat{\beta})$ using the two equations above and the empirically observed mean and variance.


My question is: what is the best/most reliable way to solve for the above system?

My attempt so far has involved the following:

Denote \begin{align} \mu \end{align} and \begin{align} \mu_{2} \end{align} as the first and second (non-central) empirical moments, respectively. From the above, we can rearrange to obtain the following system of equations to solve: \begin{align} \mu&=\frac{1}{F(x_{t})}\frac{\text{B}(x_{t};\alpha+1,\beta)}{\text{B}(\alpha,\beta)}\notag\\ &=\frac{\text{B}(x_{t};\alpha+1,\beta)}{\text{B}(x_{t};\alpha,\beta)}\notag \end{align} and \begin{align} \mu_{2}&=\frac{1}{F(x_{t})}\frac{\text{B}(x_{t};\alpha+2,\beta)}{\text{B}(\alpha,\beta)}\notag\\ &=\frac{\text{B}(x_{t};\alpha+2,\beta)}{\text{B}(x_{t};\alpha,\beta)}\notag \end{align}

To simplify the numerical optimization problem, we can take logarithmic transforms of the system: \begin{align} \log\mu-\log\text{B}(x_{t};\alpha+1,\beta)+\log\text{B}(x_{t};\alpha,\beta)&=0\notag \end{align} and \begin{align} \log\mu_{2}-\log\text{B}(x_{t};\alpha+2,\beta)+\log\text{B}(x_{t};\alpha,\beta)&=0\notag \end{align}

We implement the numerical optimization in R using the multiRoot function provided in the rootSolve package. The following code illustrates the optimization procedure.


model=function(x,parms) {


This seems to solve for the parameters most of the time, however it can be a little unreliable for certain input values. However, I do not know if it is solving for a local optimum. Moreover, in the 2-equation system above, rather than putting $F(x_{t})$ directly in (which would be a constant e.g. 0.8), I have put the actual formula for the CDF of the Beta($\alpha$,$\beta$) in evaluated at $x_{t}$. This seemed to be more stable for some odd reason.

In response to the comment Let's generate a random sample from a Beta, cap it at some quantile and then perform the optimization you suggested.




model1=function(x,pars) {
model2=function(x,pars) {

In either case, the routine more often than not fails to solve (with any combination of arguments for the function) and sometimes when it does solve, the parameters we arrive at are probably not reasonable solutions.

  • 1
    $\begingroup$ If you know $x$ such that $F(x) = 0.8$, say, you don't need the sample mean and variance in order to derive estimates, just the sample mean. This is because that $0.8$ is also a function of the parameters of the beta distribution (writing loosely.) That's likely why putting the actual formula for the CDF in adds stability to your results. $\endgroup$
    – jbowman
    Jan 15, 2018 at 15:25
  • 1
    $\begingroup$ If you could post a few of the less reliable input parameters, that would help with evaluating alternative approaches! $\endgroup$
    – jbowman
    Jan 15, 2018 at 19:20
  • 1
    $\begingroup$ I notice that you are censoring at an empirical quantile, not the theoretical quantile. Certainly reasonable, but if that's what you intend, as opposed to censoring at qbeta(0.8, a, b) for example, then I'll change my answer to match. $\endgroup$
    – jbowman
    Jan 16, 2018 at 16:17
  • 2
    $\begingroup$ I would expect maximum likelihood to work better than method of moments in many cases. MM has the singular disadvantage that it might not have any solutions, or only ridiculous solutions. It's not clear that it even should have a unique solution. It would help to clarify your question concerning what your data are like and how you determine $x_t$ and/or $F(x_t)$. Currently those details are buried in the code and it's not clear how general that code is intended to be. $\endgroup$
    – whuber
    Jan 16, 2018 at 19:33
  • 1
    $\begingroup$ I wonder at the same thing as @whuber since maximum likelihood is of the same complexity whether or not the Beta is truncated. $\endgroup$
    – Xi'an
    Jan 17, 2018 at 15:06

1 Answer 1


Your data is drawn from a censored Beta distribution, with the censoring point unknown as well as how many observations were censored. The PDF of the distribution is:

$$p(x; a, b, c) = {x^{a-1}(1-x)^{b-1} \over \int_0^c t^{a-1}(1-t)^{b-1}\text{d}t}$$

The usual Beta functions cancel out between the numerator and the denominator. This distribution evidently has three parameters; the two parameters of the uncensored Beta distribution and the censoring point $c$. Consequently, in order to use a method-of-moments estimator, we'd need to use the first three moments.

Some simulation results reported by Dishon and Weiss (1980) indicate that for the two-parameter Beta distribution the MLE is typically more efficient than the MOM estimator even for small samples unless $a=b$, as @whuber and @xi'an expected. I'd expect that adding the third moment to the MOM requirements would worsen the relative efficiency of the MOM estimator, so will continue by developing the MLE.

Taking the log of the likelihood function results in:

$$\ln L = (a-1)\sum \ln x_i + (b-1)\sum \ln (1-x_i) - n\ln\text{B}(c,a,b)$$

where $\text{B}(c,a,b)$ is the incomplete Beta function.

For our R code, we'll instead use $\ln L = \ln(p_{\beta}(x;a,b)/P_{\beta}(c;a,b))$, as base R does not, as far as I know, have an unnormalized version of the incomplete Beta function. We'll use the L-BFGS-B multivariate function minimization technique, as it allows box constraints on the parameters; an alternative would be to transform the parameters to take on any values on the real line and transform the optimization results back.

censor <- sort(xfull)[80]

n <- length(x)

lnl <- function(parms) {
  res <- sum(log(dbeta(x, parms[1], parms[2]))) - 
  if (res == Inf | res == -Inf | is.na(res)) {
    res = -9.9e99

start <- c(a, b, min(max(x)+0.001, (1+max(x))/2))
optim(par=start, fn=lnl, 
      lower=c(1e-05, 1e-05, max(x)+1e-05), 
      upper=c(999, 999, 1-1e-05),

which gives results:

[1]  1.72653952 49.99758058  0.04634924

[1] 244.2224

function gradient 
       6        6 

Not far off the actual values of $(1.7, 50, 0.04667)$, and in only six iterations.

I ran this example 100 times (with different samples each time) and the maximum number of iterations required until convergence was 29, with no convergence failures.

  • $\begingroup$ Thanks for the answer and thanks for pointing out that given the constraints the MoM approach is inadequate without considering the third moment and that the MLE approach will be much easier and more reliable. I had thought this might be the case. Cheers. $\endgroup$
    – epp
    Jan 18, 2018 at 23:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.