Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Optimation MLE

Subject: Optimation MLE

From: edward kabanyas

Date: 18 Jun, 2011 14:33:04

Message: 1 of 1

Hi All,

I want to calculate the parameter of modified gamma distribution with the maximum likelihood method (MLM). The equation of modified gamma distribution is:

f(D)=[(Lamda^(mu + 1))/gamma(mu +1)]*(D^mu)*exp(-Lambda*D) (1)

and the equation of lambda from MLE derivation is

Lambda = (mu + 1)/mean(D); (2)

and the equation of "mu" from MLE derivation is

ln(mu +1)-Psi(mu +1)=ln (mean(D)/(\prod D_i ^(1/N))), (3)

 lower bond of product in equation (3) i=1 and upper bond = N.

Assume "mu +1" in eq. (3) = Alpha, Eq. (3) can be solved by iteration by recursion as

Alpha_j+1 = Alpha_j * ln ((Alpha_j)-Psi(Alpha_j))/Y), where

Y =ln (mean(D)/(\prod D_i ^(1/N)));
Alpha_1 = (1+ sqrt(1+(4*exp(Y)/3)))/4*exp(Y)

Finally, mu = Alpha -1 and Lambda in Eq. (2)

I simulate the modified gamma distribution and try to recalculate the parameter as below:

N=1000;
 mu=0;
 Lambda=4;
 trial = 1000;
 
 N= repmat(N,1,trial);
 Lambda= repmat(Lambda,1,trial);
 mu= repmat(mu,1,trial);
 
 for j =1:1:trial;
    D = randraw('gamma', [0, 1/Lambda(j), mu(j)+1], N(j));
     b0 = mean(D);
    Capd =prod(D);
    zzln = log(b0/Capd^(1/1000));
    zz(j)=exp(zzln);
   
    muml(j) = recml(zz(j)); % estimated mu
    Ldaml(j)= ((muml(j)+1)./b0); % estimated Lambda
    
 end
 
 
 
 % function of recml
 %-----------------------
 function result = recml(zz)

 alpha1 = (1+sqrt(1+((4*zz)/3)))/(4*zz);

 iter=1;dcnorm=1.;
 while dcnorm>1E-8 & iter<1000
 alpha(1) = alpha1;
 if alpha(iter) >= 0;
 alpha(iter+1) = alpha(iter)*((log(alpha(iter))-psi(alpha(iter)))/log(zz));
 else
     alpha(iter+1) =NaN;
 end
     
     
 
 dcnorm=abs(alpha(iter+1)-alpha(iter));
 iter=iter+1;
 P1=[iter dcnorm alpha(iter)]; % Atlas results: 9.65-10.3exp(-0.6D)
 end
 
result =[(P1(1,3)-1)];
%-----------------------

Problem is muml always NaN and when I try to calculate by using the build in MLE code, the result is O.K. I think one problem is iteration in recursion and I try to solve equation (3) with fzero, but it is failed. Do you have any suggestion ? Thanks for nice help.

Cheers
Edwards

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us