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>1E8 & 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.6510.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
