Maximum Likelihood Estimator of the parameter $\alpha$ of Gamma( $\alpha$ ,1) distribution

Sarthak Bagaria

Solution to Q1: Gamma Sample Generation

We use the standard routine rgamma in R language, which uses the following method for generation of Gamma sample:

For $\alpha<1$

A rejection technique [1] is used based on the majoring function

g(x)=x^{\alpha-1}/\Gamma(\alpha)\quad\text{if }0<x\leq 1;\qquad g(x)=e^{-x}/% \Gamma(\alpha)\quad\text{if }1\leq x.

Since $e^{-x}\leq 1$ if $0<x$ and $x^{\alpha-1}\leq 1$ if $\alpha\leq 1$ and $1\leq x$ the inequality $e^{-x}x^{\alpha-1}/\Gamma(\alpha)=f(x)\leq g(x)$ is valid for all $x>0$ . The function

h(x)=x^{\alpha-1}ea/(e+a)\quad\text{if }0<x\leq 1;\qquad h(x)=e^{-x}ea/(e+a)% \quad\text{if }1\leq x

is a probability density that is proportional to g(x). Sampling from h(x) is no problem since both parts have easily invertible integrals: with a probability of $e/(e+a)$ an x below 1 (first part of h(x)) is sampled, otherwise the second part of h(x) is used. The rejection test is based on $f(x)/g(x)$ which is $e^{-x}$ or $x^{\alpha-1}$ .

For $\alpha\geq 1$

Ahrens and Dieter’s modified rejection method [2] is used.

Applying the transformation $x=(\sqrt{\alpha-1/2}+t/2)^{2}$ , the resulting transformed function g(t) is close to the standard normal density f(t). The mode of g(t) is at $t=0$ , but g(0) is a little larger than $f(0)=1/\sqrt{2\pi}$ . Also, g(t) intersects the standard normal density f(t) only once at some $t=\tau(a)<0$ . Consequently, $g(t)\geq f(t)$ for all $t\geq 0$ . This calls for the following modification of Von Neumann’s acceptance-rejection technique:
Generate a standard normal deviate T [3]. If $T\geq 0$ , accept $x=(\sqrt{\alpha-1/2}+T/2)^{2}$ as a gamma( $\alpha$ ) sample. For $T\leq\tau(\alpha)$ , where f(t) majorizes g(t), the ratio $r(T)=g(T)/f(T)$ can be compared with a (0, 1)-uniform deviate U for an ordinary rejection test. (For simplicity this test is also applied when $\tau(a)<T<0$ . In this case $r(t)>1$ and T is always accepted.) Obviously rejection occurs with probability $P(H)=\int_{-\infty}^{\tau}(f(t)-g(t))dt=\int_{\tau}^{\infty}(g(t)-f(t)).$ Hence, whenever a negative T is rejected, it must be replaced with a new $T\geq\tau(\alpha)$ , and this has to be a sample from the difference distribution whose probability density function is proportional to $g(t)-f(t)$ in [ $\tau,\infty$ ). Sampling from this difference is done by means of a double-exponential hat.

Solution to Q2: Maximum Likelihood Estimation of $\alpha$

We use the standard routine fitdistr in R language, with BFGS method. For Gamma distribution reasonable starting values for the iterative method can be computed by the routine itself.

BFGS Method [4]

In Newton’s method, we find the new iterate $x_{k+1}$ as a function of $x_{k}$ as follows. For any point x define $p=x-x_{k}$ , the second order Taylor expansion around $x_{k}$ is given by

m_{k}(p)=f_{k}+p^{T}\nabla f_{k}+\frac{1}{2}p^{T}B_{k}p

This defines a quadratic model of the function near the point $x_{k}$ . Its gradient with respect to x is $m_{k}(p)=\nabla f_{k}+B_{k}p$ , and it is minimized at $p_{k}=-B_{k}^{-1}\nabla f_{k}$ .

Working with the inverse Hessian $H_{k}$ in place of $B_{k}$ , the secant equation becomes $H_{k}y_{k-1}=s_{k-1}$ . The optimization is then : minimize $||H-H_{k-1}||_{W}$ subject to $H=H^{T}$ , $Hy_{k-1}=s_{k-1}$ , which has the unique solution

H_{k}=(I-p_{k}s_{k-1}y^{T}_{k-1})H_{k-1}(I-p_{k-1}y_{k-1}s^{T}_{k-1})+s_{k-1}p% _{k-1}s^{T}_{k-1}

where $s_{k-1}=x_{k}-x_{k-1}$ , $y_{k-1}=\nabla f_{k}-\nabla f_{k-1}$ , W is any matrix satisfying $Wy_{k-1}=s_{k-1}$ , and $||H-H_{k-1}||_{W}=||W^{\frac{1}{2}}(H-H_{k-1})W^{\frac{1}{2}}||$ .

Each step of the BFGS method has the form

x_{k+1}=x_{k}-\alpha_{k}H_{k}\nabla f_{k},\qquad k=0,1,2,\ldots

where $\alpha_{k}$ is computed from a line search procedure to satisfy the Wolfe conditions:

f(x_{k}+\alpha_{k}p_{k})\leq f(x_{k})+c_{1}\alpha_{k}\nabla f_{k}^{T}p_{k},

\nabla f(x_{k}+\alpha_{k}p_{k})^{T}p_{k}\geq c_{2}\nabla f_{k}^{T}p_{k}

with $0<c_{1}<c_{2}<1$ .

Solution to Q3 & Q4: Results

The following graphs were plotted for values of MLE of $\alpha$ against the values of $\alpha_{0}$ , the parameter used for generation of sample. $\alpha_{0}$ ranges from .1 to 10 in steps of .1, and for each value of $\alpha_{0}$ , estimates for 20 samples are calculated. n is the number of Gamma variates in a single sample.

Example outputs from different samples with $\alpha$ for $\alpha_{0}=5$ :

$n=10:$
$5.3144648(0.6950958);\quad 4.5355942(0.6368612);\quad 4.6917132(0.6489457);$
$5.743960(0.725238);\qquad 4.0315324(0.5962093);\quad 5.2970827(0.6938487);$

$n=50:$
$4.6669183(0.2893656);\qquad 4.6980724(0.2904354);\qquad 4.9866184(0.3001651);$
$4.7338926(0.2916606);\qquad 4.7629318(0.2926501);\qquad 4.927377(0.298193);$

$n=100:$
$5.4493586(0.2228456);\qquad 5.1090057(0.2151015);\qquad 5.5133350(0.2242719);$
$4.9219525(0.2107262);\qquad 5.2803015(0.2190328);\qquad 5.3109010(0.2197278);$

$n=500:$
$4.99358132(0.09499359);\quad 5.15714169(0.09669357);\quad 5.15428979(0.0966641% 8);$
$5.11426871(0.09625083);\quad 5.01874406(0.09525706);\quad 4.95975300(0.0946382% 4);$

The quantities in brackets are the estimated standard errors due to numerical optimization, which, we see, decrease as the sample size n is increased. Also, the variance of the estimates is less for large sample sizes.

And though we have different algorithms for sample generation for the cases $\alpha_{0}<1$ and $\alpha_{0}\geq 1$ , the sampling and estimation methods do not differentiate between integer and non-integer $\alpha_{0}$ values, and the precision achieved is similar in two cases.

For example, when $n=50$ and $\alpha=5.1$ :
$5.3183637(0.3109812);\qquad 4.9230182(0.2980474);\qquad 4.8010873(0.2939454);$
$5.1262415(0.3047635);\qquad 5.3411929(0.3117119);\qquad 5.0274508(0.3015171);$

Examples of generated Samples:

Code for plots (in R language)


library(MASS)
n=50                       #number of variates in each sample
alpha=.1*(1:100)           #range of alpha for sample generation
s=20Ψ                      #number of samples for each value of alpha
assign("x", c())
assign("y", c())
assign("z", c())
for (a in alpha)
  for (i in 1:s) {
    x <- c(x,(a))
    y <- c(y,fitdistr(rgamma(n,(a)), "gamma", rate=1, method="BFGS"))
  }
for (i in 1:(s*length(alpha))) z[i] = y[[5*i-4]]
plot(
 Ψx,z, type="l", main=expression(paste("MLE of ",alpha)),
 Ψsub=substitute(paste("n = ", n), list(n=n)),
 Ψxlab=expression(alpha[0]), ylab=expression(alpha)
Ψ)

Code for single estimation (in R language)


n=50Ψ                      #number of variates in each sample
a=2                        #value of alpha for generation of sample
fitdistr(rgamma(n,a), "gamma", rate=1, method="BFGS")

References

[1] Ahrens, J. H. and Dieter, U. (1974). Computer methods for sampling from gamma, beta, Poisson and binomial distributions. Computing, 12, 223–246
[2] Ahrens, J. H. and Dieter, U. (1982). Generating gamma variates by a modified rejection technique. Communications of the ACM, 25, 47–54
[3] Ahrens, J. H. and Dieter, U. (1972). Computer methods for sampling from the exponential and normal distributions. Communications of the ACM, 15, 873-882
[4] Jorge Nocedal and Stephen J. Wright (2000). Numerical Optimization

Maximum Likelihood Estimator of the parameter α of Gamma(α,1) distribution