Self-study: Finding the maximum likelihood estimates of the parameters of a density function - UPDATED

UPDATED

I am trying to find maximum likelihood estimation of a probability distribution function given below

\begin{equation} g(x)=\frac{1}{\Gamma \left( \alpha \right)\gamma^{2\alpha}2^{\alpha-1}}x^{2\alpha-1}\exp\left\{{-\frac{x^2}{2\gamma^{2}}}\right\}I_{{\rm I\!R}^{+}}(x) \end{equation}

where $\alpha >0$ is the shape parameter, $\sigma >0$ is scale parameter.

The likelihood function is then given by

\begin{equation} L(\alpha,\gamma/x)=\prod\limits_{i=1}^{n}\frac{1}{\Gamma \left( \alpha \right)\gamma^{2\alpha}2^{\alpha-1}}x_i^{2\alpha-1}\exp\left\{{-\frac{x_i^2}{2\gamma^{2}}}\right\} \end{equation}

Thus, the complete likelihood function is then \begin{equation} L(\alpha,\gamma/x)=\frac{1}{[\Gamma \left( \alpha \right)]^{n}\gamma^{2\alpha n}{2^{n\alpha-n}}} \exp\left\{{-\frac{1}{2\gamma^{2}}\sum\limits_{i=1}^{n}x_{i}^{2}}\right\}\left(\prod\limits_{i=1}^{n}x_{i}\right)^{2\alpha-1} \end{equation}

Now, the log-likelihood function denoted by $\ell$ is

\begin{equation} \begin{aligned} \ell &=\log[L(\alpha,\gamma/x)]\\ &=-n\log(\Gamma \left( \alpha \right))-2\alpha n \log(\gamma)-n(\alpha-1)\log(2)-\frac{1}{2\gamma^{2}}\sum\limits_{i=1}^{n}x_{i}^{2}+(2\alpha-1)\sum\limits_{i=1}^{n}\log(x_{i}) \end{aligned} \end{equation}

The entries of the score function are given by

\begin{equation} \begin{aligned} \frac{\partial \ell}{\partial \alpha}=-n\psi(\alpha)-2n\log(\gamma)-n\log(2)+2\sum\limits_{i=1}^{n}\log(x_{i}) \end{aligned} \end{equation} where $\psi(\alpha)$ is the digamma function and

\begin{equation} \begin{aligned} \frac{\partial \ell}{\partial \gamma}=-\frac{2\alpha n}{\gamma}+\frac{\sum\limits_{i=1}^{n}x_{i}^{2}}{\gamma^{3}} \end{aligned} \end{equation}

Setting these two equations to zero and solving them simultaneously results in maximum likelihood estimates (MLE) of parameters, $\hat{\alpha}$ and $\hat{\gamma}$. However, the equations obtained by setting the above partial derivatives to zero are not in closed form and the values of parameters $\alpha$ and $\gamma$ must be found using iterative methods.

Fisher information matrix is defined as $I_{ij}=-E\left\{\frac{\partial^{2} \ell}{\partial \tau_i \partial \tau_j} \log[L(x_i, \vec{\tau})]\ \right\}$ where $\tau_1=\alpha$ and $\tau_2=\gamma$. Thus, information matrix for gamma-rayleigh distribution is given by,

\begin{equation} I=n \left[ \begin{array}{cc} \psi_{1}(\alpha) & 2/\gamma\\ 2/\gamma & 4\alpha/\gamma^2 \end{array} \right] \end{equation}

I am trying to use Fisher Scoring to find MLEs of the parameters. Here is my MATLAB code. I first generate 1000 random observations from gamma-distribution and run this code. My starting values and the rest are given in the code.

clear all;
clc;

%Simulate 1000 sample from Gamma Distribution
n=1000;
alpha=3;
lambda=0.05;
x=gamrnd(alpha,1/lambda,1,n);

figure(1)
histfit(x,8,'gam');

sumlogx=sum(log(x)); sumxsquare=sum(x.^2);

%Initial Values
alpha=mean(x)^2/var(x);
gam=mean(x)/var(x);
theta=[alpha; gam];
S=Inf;

while sum(abs(S) > 10^(-5)) > 0
S=[-n*psi(theta(1))-2*n*log(theta(2))-n*log(2)+2*sumlogx;...
(-2*theta(1)*n/theta(2))+(sumxsquare/(theta(2)^3))];
FIM=n*[psi(1, theta(1)), 2/theta(2);...
2/theta(2), 4*theta(1)/(theta(2)^2)];
theta=theta + FIM\S;
end

alpha_hat=theta(1)
gam_hat=theta(2)

fprintf('alpha_hat=%g, gamma_hat=%g \n', theta(1),theta(2))

But for some reasons I cannot figure out, I am getting "Error using psi X must be nonnegative." error. My $\alpha$ values are being negative in the iteration somehow and I do not know how to fix it!

I am also running Newton-Raphson whose MATLAB code is given below

clear all;
clc;

%Simulate 100 sample from Gamma Distribution
n=1000;
alpha=3;
lambda=0.05;
x=gamrnd(alpha,1/lambda,1,n);

figure(1)
histfit(x,8,'gam');

sumlogx=sum(log(x)); sumxsquare=sum(x.^2);

%tuning parameters scale=gamma; shape=alpha
itermin=10^-7;
maxiter=10^7;
sc_init=0.000001;
sh_init=0.000001;
converged=[0;0;sc_init;sh_init];

% pdf
pdf=@(x,gam,alpha) 1/(gamma(alpha)*(gam^(2*alpha))*(2^(alpha-1)))*(x^(2*alpha-1))*exp(-(x^2)/(2*(gam^2)));

%score function is the first partial derivative of the log likelihood function
score=@(gam,alpha) -n*psi(alpha)-2*n*log(gam)-n*log(2)+2*sumlogx;

%Hessian function is the negative of the 2nd
hessian=@(gam,alpha) psi(1, alpha);

sc_loop=2;
scale_hat=zeros(1,maxiter);
scale_hat(1)=sc_init;

while 1==1
sh_loop=2;
shape_hat=zeros(1,maxiter);
shape_hat(1)=sh_init;

while 1==1
%calculate chat as chat_prev+score(chat_prev)/hessian(chat_prev)
shape_hat(sh_loop)=shape_hat(sh_loop-1)+score(scale_hat(sc_loop-1),shape_hat(sh_loop-1))/hessian(scale_hat(sc_loop-1),shape_hat(sh_loop-1));
%test for a convergence
if abs(shape_hat(sh_loop)-shape_hat(sh_loop-1))<itermin
break %the process converged to a c value
elseif sh_loop>maxiter
disp(['max iteration on \alpha achieved:', num2str(maxiter)]);
return
end
sh_loop=sh_loop+1;
end

scale_hat(sc_loop)=(sum(x.^shape_hat(sh_loop-1))/n)^(1/shape_hat(sh_loop-1));
%test for a convergence
if abs(scale_hat(sc_loop)-scale_hat(sc_loop-1))<itermin
break %the process converged to a gamma value
end

converged=[converged,[sc_loop-1;sh_loop-1;scale_hat(sc_loop);shape_hat(sh_loop)]];
sc_loop=sc_loop+1;
end

%final display
disp(repmat('-',[1,30])),disp(' Iteration Scale Shape'),disp(repmat('-',[1,30]))
disp(num2str(converged','%6.4f')),disp(repmat('-',[1,30]))
disp(['Real values: gamma=', num2str(gam),',alpha=',num2str(alpha)])

I am getting the same "Error using psi, X must be nonnegative." error.

Could you help me about it? Something is wrong with psi function and I do not know. Maybe I should use approximation but I am not sure how much of the information that I will loose!

• There is no latent variable that is obvious for this model, so it seems difficult to apply EM... Dec 19 '14 at 18:05
• This is what I think but there are many many articles finding new distribution and using either EM or Newton-Raphson to find maximum likelihood estimations. At least, they say so but they do not give any clues about it. they only mention that they pick some sort of initial values.
– ARAT
Dec 19 '14 at 18:22
• Do you mean there are a lot of distributions for which EM applies or a lot of papers applying EM to your distribution? Dec 19 '14 at 20:18
• I mean every single day a new distribution is found or generated. For the parameter estimations, the authors use either one of these algorithms, at least they say so but they do not give details such as which latent variable they opt to choose for em.
– ARAT
Dec 19 '14 at 20:37
• EM only applies to cases where there is an exponential family on a latent variable, so this is a very special case. Dec 19 '14 at 20:48

[Note: This is my answer to the Dec. 19, 2014, version of the question.]

If you operate the change of variable $y=x^2$ in your density $$f_X(x|\alpha,\beta,\sigma)=\frac{1}{\Gamma \left( \alpha \right)\beta^{\alpha}}\exp\left\{{-\frac{x^2}{2\sigma^{2}}\frac{1}{\beta}}\right\}\frac{x^{2\alpha-1}}{2^{\alpha-1}\sigma^{2\alpha}}\mathbb{I}_{{\mathbb{R}}^{+}}(x)$$ the Jacobian is given by $\dfrac{\text{d}y}{\text{d}x}= 2x = 2y^{1/2}$ and hence \begin{align*} f_Y(y|\alpha,\beta,\sigma)&=\frac{1}{\Gamma \left( \alpha \right)\beta^{\alpha}}\exp\left\{{-\frac{y}{2\sigma^{2}}\frac{1}{\beta}}\right\}\frac{y^{\frac{2\alpha-1}{2}}}{2^{\alpha-1}\sigma^{2\alpha}}\frac{1}{2 y^{1/2}}\mathbb{I}_{{\mathbb{R}}^{+}}(y)\\ &=\frac{1}{\Gamma \left( \alpha \right)\beta^{\alpha}}\exp\left\{{-\frac{y}{2\sigma^{2}}\frac{1}{\beta}}\right\}\frac{y^{{\alpha-1}}}{2^{\alpha}\sigma^{2\alpha}}\mathbb{I}_{{\mathbb{R}}^{+}}(y) \end{align*} This shows that

1. This is a standard $\mathcal{G}(\alpha,2\sigma^2\beta)$ model, i.e. you observe $$(x_1^2,\ldots,x_n^2)=(y_1,\ldots,y_n)\stackrel{\text{iid}}{\sim}\mathcal{G}(\alpha,\eta);$$
2. the model is over-parametrised since only $\eta=2\sigma^2\beta$ can be identified;
3. EM is not necessary to find the MLE of $(\alpha,\eta)$, which is not available in closed form but solution of$$\hat\eta^{-1}=\bar{y}/\hat{\alpha}n\qquad\log(\hat{\alpha})-\psi(\hat{\alpha})=\log(\bar{y})-\frac{1}{n}\sum_{i=1}^n\log(y_i)$$ where $\psi(\cdot)$ is the di-gamma function. This paper by Thomas Minka indicates fast approximations to the resolution of the above equation.
• Much obliged, Xi'an. But how this overparametrization effects the model estimation of the density function of mine? Secondly, about maximum likelihood estimation of the $\alpha$ parameter of gamma distribution, I still need some sort of iteration method to find it, is not it right? Does not it make things more complicated?
– ARAT
Dec 19 '14 at 18:15
• I edited my answer: (a) you can only estimate $\alpha$ and $\eta$, but cannot separate $\beta$ from $\sigma$; (b) the likelihood equations are now given; (c) it certainly makes the problem easier. Dec 19 '14 at 19:49
• I wonder at why you reopen this problem when considering this is a standard Gamma model, for which MLEs are provided by standard software. Oct 26 '15 at 9:18

In the case of the EM algorithm, the initial values can be set arbitrarily since the iterations are guaranteed to converge to the maximum:

We have seen that both the E and the M steps of the EM algorithm are increasing the value of a well-defined bound on the log likelihood function and that the complete EM cycle will change the model parameters in such a way as to cause the log likelihood to increase (unless it is already at a maximum, in which case the parameters remain unchanged).

There are several strategies to pick the initial values that will improve the overall performance of the algorithms:

Note that the EM algorithm takes many more iterations to reach (approximate) convergence compared with the K-means algorithm, and that each cycle requires significantly more computation. It is therefore common to run the K-means algorithm in order to find a suitable initialization for a Gaussian mixture model that is subsequently adapted using EM. The covariance matrices can conveniently be initialized to the sample covariances of the clusters found by the K-means algorithm, and the mixing coefficients can be set to the fractions of data points assigned to the respective clusters.

Bishop, C. M. (2006). Pattern recognition and machine learning (Vol. 1, p. 740). New York: springer.

You can read on this in chapter 9 of the book.

Using EM to estimate latent variables in the context of Gaussian Mixture Models with two components ($\Delta_i = 0$ and $\Delta_i = 1$), loglikelihood is given by (I'm sorry I'm in a hurry and can't type this myself)

ESLII gives some advice on how to select the initial values:

A good way to construct initial guesses for $\mu_1$ and $\mu_2$ is simply to choose two of the $y_i$ at random. Both $\sigma_1^2$ and $\sigma_2^2$ can be set equal to the overall sample variance. The mixing proportion $\hat{\pi}$ can be started at the value 0.5.

Hastie, T., Tibshirani, R.,, Friedman, J. (2008). The elements of statistical learning: data mining, inference and prediction. Springer.

I insist, you can pick any arbitrary initial values, but some will converge faster.

• Thanks mugen for the reference. I will check it out. But what about the latent variable? How to define it for this problem?
– ARAT
Dec 19 '14 at 17:50
• @MuratArat Please note that the vector of weights in the mixture model $\pi$ is just another parameter of the distribution. That means that you use the same iterative process to discover all the vectors weights, the vector of means and the variance and covariance matrix. You start with an arbitrary value and improve the estimation after each iteration. I'll edit the answer to include this comment. Dec 19 '14 at 18:02
• I will follow these advices about initial values. however, I am sory if I understand so slowly about choosing latent variable. What you wrote is a Gaussian mixture and I am well aware of the EM method for it. How about the problem I have? It's not a mixture of any kind of distributions.
– ARAT
Dec 19 '14 at 18:18
• @MuratArat oh, you're clearly right! I jumped directly into the answer without paying enough attention to the question. I assumed you had a more general question but you're clearly looking for a solution for this specific exercise. Sorry! Dec 19 '14 at 18:21
• Actually, not for this specific exercise. This was a whole trouble from the beginning I learned EM. Because as I wrote in the comments of my very first post, I saw some articles that they say they use EM algorithm but they do not give any clue about which latent variable they created (is it the right verb?) to find MLEs of the parameters.
– ARAT
Dec 19 '14 at 18:25