# Truncated Beta parameters - method of moments

Context

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.

Question

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.

m=mean(data)
v=var(data)
xt=max(data)
coef=c(m,v,xt)

model=function(x,parms) {
F1=log(parms[1])-log(pbeta(parms[3],x[1]+1,x[2]))-log(beta(x[1]+1,x[2]))+log(parms[4])+log(beta(x[1],x[2]))
F2=log(parms[2]+parms[1]\^{}2)-log(pbeta(parms[3],x[1]+2,x[2]))-log(beta(x[1]+2,x[2]))+log(parms[4])+log(beta(x[1],x[2]))
c(F1=F1,F2=F2)
}

ss=multiroot(f=model,start=(start1,start2),parms=coef,rtol=1e-15,atol=1e-15,ctol=1e-15,useFortran=FALSE,positive=TRUE,maxiter=10000)


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.

a=1.7
b=50
xfull=rbeta(100,a,b)
alpha=0.8

x=xfull[xfull<sort(xfull)[80]]

mu=mean(x)
v=var(x)
xt=max(x)
coef=c(mu,v,xt,alpha);coef

model1=function(x,pars) {
F1=pbeta(pars[3],x[1],x[2])*beta(x[1],x[2])-pars[4]
F2=pbeta(pars[3],x[1]+1,x[2])*beta(x[1]+1,x[2])-pars[4]*pars[1]
c(F1=F1,F2=F2)
}
model2=function(x,pars) {
F1=log(pbeta(pars[3],x[1],x[2]))+log(beta(x[1],x[2]))-log(pars[4])
F2=log(pbeta(pars[3],x[1]+1,x[2]))+log(beta(x[1]+1,x[2]))-log(pars[4])-log(pars[1])
c(F1=F1,F2=F2)
}
ss=multiroot(f=model1,start=c(a,b),pars=coef,positive=TRUE,maxiter=10000,rtol=1e-12,atol=1e-12,ctol=1e-12,useFortran=FALSE)
ss
ss=multiroot(f=model2,start=c(a,b),pars=coef,positive=TRUE,maxiter=10000,rtol=1e-12,atol=1e-12,ctol=1e-12,useFortran=FALSE)
ss


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.

• 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. – jbowman Jan 15 '18 at 15:25
• If you could post a few of the less reliable input parameters, that would help with evaluating alternative approaches! – jbowman Jan 15 '18 at 19:20
• 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. – jbowman Jan 16 '18 at 16:17
• 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. – whuber Jan 16 '18 at 19:33
• I wonder at the same thing as @whuber since maximum likelihood is of the same complexity whether or not the Beta is truncated. – Xi'an Jan 17 '18 at 15:06

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.

a=1.7
b=50
xfull=rbeta(100,a,b)
censor <- sort(xfull)[80]

x=xfull[xfull<censor]
n <- length(x)

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

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),
control=list(fnscale=-1),
method="L-BFGS-B")


which gives results:

$par [1] 1.72653952 49.99758058 0.04634924$value
[1] 244.2224

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

• 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. – epp Jan 18 '18 at 23:21