MATLAB Central Newsreader - MLE - use with Discrete Weibull mean function

MLE - use with Discrete Weibull mean function

Erik Gadzinski — Sat, 07 Nov 2009 18:31:02 -0500

I am having trouble with an MLE function that was passed along to me. I am new to MATLAB and am having trouble shooting problems. I am taking 'v' values that are randomly generated that are discrete numbers that range from 0 to a couple thousand. In the code, "v' can have a record length of 10, 50, 100, 1000, and 10,000. In the Llikef function I would like to take the log of the function, but I have been getting log zero when large numbers usually over 1,022 are entered into the function. So I took the log out. The thing I need help with is that I want my q and N (eta) values to be positive. The function sometimes outputs negative q's and eta's. Does anybody know how to do this? Thank you in advance!

In the function that(1) = q and that(2) = eta.
vsim are whole numbers that range from 0 to a few thousand

%find the Mean MLE Fitted Average.
%1) take each record length of 'v' values and plug it into the MLE routine
%2) use the q and eta values generated from that function and plug them
%in the DW mean function

slength=1;
elength=recordlength;
meanindex=1;
vsim;
while meanindex<=100
    if sum(vsim(slength:elength))>0
        x=vsim(slength:elength);

        Llikef = @(that)-sum(max(that(1),.00001).^(x.^max(that(2),.00001))-max(that(1),.00001).^((x+1).^max(that(2),.00001)));
        %Llikef = @(that)-sum(log(max(that(1),.00001).^(x.^max(that(2),.00001))-max(that(1),.00001).^((x+1).^max(that(2),.00001))))

        that0 = [0.5 1];
        [that,Llikef] = fminsearch(Llikef,that0);

        if Llikef<=0
            that(1)=.5;
            that(2)=1;
        else end
            if that(1)>0
                  q2(meanindex)=that(1);
               else q2(meanindex)=.5;
               end
               if that(2)>0;
                  eta2(meanindex)=that(2);
               else eta2(meanindex)=1;
               end

        N=eta2(meanindex);
        q=q2(meanindex);
        M=1e5;
        vm=1:M;
        mmlefm(meanindex)= sum(q.^(vm.^N))+(gamma(1/N)*gammainc((M+1)^N*(-log(q)),1/N,'upper'))/(N*(-log(q))^(1/N));
        meanindex=meanindex+1;
        slength=elength+1;
        elength=elength+recordlength;
    else
        q2(meanindex)=.5;
        eta2(meanindex)=1;
        mmlefm(meanindex)=0;
        slength=elength+1;
        elength=elength+recordlength;
        meanindex=meanindex+1;
    end
end

Re: MLE - use with Discrete Weibull mean function

Erik Gadzinski — Sat, 07 Nov 2009 18:59:01 -0500

"Erik Gadzinski" <ergadz@yahoo.com> wrote in message <hd4ed6$l92$1@fred.mathworks.com>...
> I am having trouble with an MLE function that was passed along to me. I am new to MATLAB and am having trouble shooting problems. I am taking 'v' values that are randomly generated that are discrete numbers that range from 0 to a couple thousand. In the code, "v' can have a record length of 10, 50, 100, 1000, and 10,000. In the Llikef function I would like to take the log of the function, but I have been getting log zero when large numbers usually over 1,022 are entered into the function. So I took the log out. The thing I need help with is that I want my q and N (eta) values to be positive. The function sometimes outputs negative q's and eta's. Does anybody know how to do this? Thank you in advance!

In short, how do I make sure that "that" comes back postive in the two lines below: (the rest of the code is below that)

that0 = [0.5 1];
[that,Llikef] = fminsearch(Llikef,that0)
>
> In the function that(1) = q and that(2) = eta.
> vsim are whole numbers that range from 0 to a few thousand
>
> %find the Mean MLE Fitted Average.
> %1) take each record length of 'v' values and plug it into the MLE routine
> %2) use the q and eta values generated from that function and plug them
> %in the DW mean function
>
>
> slength=1;
> elength=recordlength;
> meanindex=1;
> vsim;
> while meanindex<=100
> if sum(vsim(slength:elength))>0
> x=vsim(slength:elength);
>
> Llikef = @(that)-sum(max(that(1),.00001).^(x.^max(that(2),.00001))-max(that(1),.00001).^((x+1).^max(that(2),.00001)));
> %Llikef = @(that)-sum(log(max(that(1),.00001).^(x.^max(that(2),.00001))-max(that(1),.00001).^((x+1).^max(that(2),.00001))))
>
> that0 = [0.5 1];
> [that,Llikef] = fminsearch(Llikef,that0);
>
> if Llikef<=0
> that(1)=.5;
> that(2)=1;
> else end
> if that(1)>0
> q2(meanindex)=that(1);
> else q2(meanindex)=.5;
> end
> if that(2)>0;
> eta2(meanindex)=that(2);
> else eta2(meanindex)=1;
> end
>
> N=eta2(meanindex);
> q=q2(meanindex);
> M=1e5;
> vm=1:M;
> mmlefm(meanindex)= sum(q.^(vm.^N))+(gamma(1/N)*gammainc((M+1)^N*(-log(q)),1/N,'upper'))/(N*(-log(q))^(1/N));
> meanindex=meanindex+1;
> slength=elength+1;
> elength=elength+recordlength;
> else
> q2(meanindex)=.5;
> eta2(meanindex)=1;
> mmlefm(meanindex)=0;
> slength=elength+1;
> elength=elength+recordlength;
> meanindex=meanindex+1;
> end
> end

Re: MLE - use with Discrete Weibull mean function

Peter Perkins — Wed, 11 Nov 2009 15:17:19 -0500

Erik Gadzinski wrote:
> I am having trouble with an MLE function that was passed along to me. I am new to MATLAB and am having trouble shooting problems. I am taking 'v' values that are randomly generated that are discrete numbers that range from 0 to a couple thousand. In the code, "v' can have a record length of 10, 50, 100, 1000, and 10,000. In the Llikef function I would like to take the log of the function, but I have been getting log zero when large numbers usually over 1,022 are entered into the function. So I took the log out. The thing I need help with is that I want my q and N (eta) values to be positive. The function sometimes outputs negative q's and eta's. Does anybody know how to do this? Thank you in advance!

Erik, I can't exactly tell what your code is intending to do, but you seem to be describing two problems, and perhaps confusing them (they both seem to involve logs):

> In the Llikef function I would like to take the log of the function, but I have been getting log zero

You may be referring to a problem where you want to compute the log-likelihood, and do so by first computing the liklihood,and then compute it's log, and find that the likelihood has underflowed to zero. People usually work _directly_ with the log-likeliohood for just that reason. So, for eaxmple, the normal log-liklihood is not computed as

log(prod(exp(-.5*(x-mu).^2)/sigma))

but as sum(-.5*(x-mu).^2 - log(sigma))

or somesuch.

> The function sometimes outputs negative q's and eta's.

If that's where the (log-)likelihood surface has it's maximum, that's where FMINSEARCH will go. It may also be that FMINSEARCH is just taking too big of a step at times. In any case, reparameterize your model to maximize the log-likelihood with respect to log(q) and log(eta). Thos new parameters are unrestricted, and when you exponetiate them to get q and eta, you are guaranteed positive values.

Hope this helps.