SAS/IML Software and Matrix Computations

Statistical programming, matrix languages, and more
BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
AndreMenezes
Fluorite | Level 6

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.

 

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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

View solution in original post

3 REPLIES 3
Rick_SAS
SAS Super FREQ

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

AndreMenezes
Fluorite | Level 6

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;
Rick_SAS
SAS Super FREQ

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"

sas-innovate-white.png

Special offer for SAS Communities members

Save $250 on SAS Innovate and get a free advance copy of the new SAS For Dummies book! Use the code "SASforDummies" to register. Don't miss out, May 6-9, in Orlando, Florida.

 

View the full agenda.

Register now!

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 3 replies
  • 2895 views
  • 3 likes
  • 2 in conversation