topic Re: Maximum likelihood estimation: Inverse Gamma Distribution in SAS/IML Software and Matrix Computations

Maximum likelihood estimation: Inverse Gamma Distribution

AndreMenezes — Sun, 07 Aug 2016 13:04:27 GMT

I'm tryng to find MLE from inverse gamma distribution using SAS/IML, however when I run optmization appear an error. I suposse the error is because the function 'l' underflow. I have seen the Rick's blog about MLE (http://blogs.sas.com/content/iml/2011/10/12/maximum-likelihood-estimation-in-sasiml.html) and write this code:

proc iml;
/*Quantile*/
start qinvgama(p,alpha,beta);
    qf = 1/quantile("GAMMA",1-p,alpha,beta);
    return(qf);
finish;
/*Variate*/
start rinvgama(n,alpha,beta);
    u  = j(n,1);                
    call randgen(u, "Uniform");
    rg     = qinvgama(u,alpha,beta);
    return(rg);
finish;

start MLE(par) global (x);
   alpha    =     par[1];
   beta     =     par[2];
   n        =    nrow(x);
   l        = n#(alpha#log(beta) - log(gamma(alpha))) - beta#sum(1/x) - (alpha + 1)#sum(log(x));
   return (l);
finish;

x = rinvgama(100,2,3);
sup = { 0   0,  
        .   .};
ini = {1.2 3};
opt = {1, 4};
call nlpnra(it, resmle, "MLE", ini, opt, sup);

print resmle;

quit;

ERROR: (execution) Invalid argument to function.

Re: Maximum likelihood estimation: Inverse Gamma Distribution

Rick_SAS — Sun, 07 Aug 2016 17:28:49 GMT

I think the error message is telling you that the parameters need to be strictly positive. Try

sup = { 1e-8 1e-8, 
         .    .};

Also, there is an easier way to generate from the inverse gamma:

If X ~ gamma(a, b) then 1/X ~ inverse-gamma(a, 1/b)

Some other facts about the inverse gamma distribution is available in the MCMC documentation:

- Definition (search for igamma)

- Potential confusion of parameters

Re: Maximum likelihood estimation: Inverse Gamma Distribution

AndreMenezes — Sun, 07 Aug 2016 19:38:48 GMT

Thank you very much Rick. I change the way to generate random variate from the inverse gamma distribution. Now I'm using this:

start rinvgama(n,alpha,beta);
	aux	= j(n,1);                
	call randgen(aux, "Gamma", alpha, 1/beta);
	rg	= 1/aux; 		
	return(rg);
finish;

Re: Maximum likelihood estimation: Inverse Gamma Distribution

Rick_SAS — Tue, 16 Aug 2016 00:51:22 GMT

While searching my blog for something, I realized that I blogged about how to simulate from the inverse gamma distriution in 2014:

"Simulating from the Inverse Gamma Distribution in SAS"