MLE with weighted in Gamma Distribution ( error : call NLPFDD)

[babaJB]

I am trying to find the parameters in Gamma distribution using maximum likelihood estimation adjusted to allow for the weights w_i. IML successfully gives a resulting parameter but but when I check the results of these parameters with the parameters that exist on paper the results showed a highly significant difference. So i am trying to check the derivative formula, IML gives an ERROR message saying Matrix has not been set to a value. below is the code and output error in call NLPFDD. Thanks for the help.

proc iml;

klaim = { 1268 2225 2323 3135 3639 3815 3935 4187 4792 4960 5212 5253 5395 5953 6094 6210 6251 6319 6401 6615 6801 6905 7048 7062 7418 7589 7647 7886 8573 8661 9244 9515 9553 10043 10160 11826 11905 12968 13265 13609 16659 16756 18455 19875 22278 34091};

w ={ 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 0.478 4.480 0.478 0.478 0.478 0.478 0.478 0.478 4.480 0.478 0.478 4.480 0.478 4.480 4.480};

print klaim w;

p = ncol(klaim);

/*----------------------------------------------------------------*/

start f_gamma(x) global(klaim,w);

/* x[1]=alpha dan x[2]=theta*/

p = ncol(klaim);

sum1=0.;sum2=0.;sum3=0.;

do i =1 to p;

        x_i = klaim#w; /* assign wi*xi */

        ai=klaim; /* assign for sum of klaim*/

        bi=w;

        xi_log = log(ai)*bi; /* untuk log(x_i) */

        sum1 = sum1+x_i;

        sum2 = sum2+xi_log;

        sum3 = sum3+bi;

end;

f = sum3*(-p*x[1]*log(x[2])-p*lgamma(x[1]))+ (x[1]-1)*sum2 - sum1/x[2];

return(f);

finish f_gamma;

start g_gamma(x) global(klaim,w);

p = ncol(klaim);

g = j(1,2,0.); /* matrix 1x2 with 0 entries  */

sum1=0.; sum2=0.; sum3=0.;

do i =1 to p;

        x_i = klaim#w; /* assign x_i wi*xi */

        ai=klaim; /* assign ai for sum of klaim*/

        bi=w;

        xi_log = log(ai)*bi;

        sum1 = sum1+x_i/(x[2]**2);

        sum2  = sum2+xi_log;

        sum3  = sum3+bi;

end;

g[1] =  sum3*(-p * log(x[2])-p*digamma(x[1]))+sum2;

g[2] =  sum3*( -p * x[1] /x[2] )+ sum1;

return(g);

finish g_gamma;

n=2;

x0=j(1,n,.5);

optn = {1 2};

con = { 1.e-6 1.e-6,

1.e6 1.e6};

call nlpnra(rc,xres,"f_gamma",x0,optn,con,,,,"g_gamma");

xopt=xres;

fopt=f_gamma(xopt);

w = {2.4 5000}; /* evaluate derivative at this (alpha, theta) */

grad = g_gamma(x);

call nlpfdd (func, FDDGrad, hessian, "f_gamma", w); 

diff = max(abs(grad - FDDGrad));

print grad, FDDGrad, diff;

############ Error Output

2107  w = {2.4 5000};

2107!                 /* evaluate derivative at this (mu, sigma) */

2108  grad = g_gamma(x);

ERROR: (execution) Matrix has not been set to a value.

operation : G_GAMMA at line 2108 column 15

operands  : x

x      0 row       0 col     (type ?, size 0)

statement : ASSIGN at line 2108 column 1

2109  call nlpfdd(func, FDDGrad, hessian, "f_gamma", w);

ERROR: (execution) Invalid subscript or subscript out of range.

operation : [ at line 2068 column 25

operands  : w, i

w      1 row       2 cols    (numeric)

       2.4      5000

i      1 row       1 col     (numeric)

         3

statement : ASSIGN at line 2068 column 9

traceback : module F_GAMMA at line 2068 column 9

ERROR: NLPFDD call: Objective function cannot be computed at X0.

NOTE: NLPFDD call: You must provide a feasible point X0.

ERROR: (execution) Invalid subscript or subscript out of range.

operation : NLPFDD at line 2109 column 1

operands  : *LIT1036, w

*LIT1036      1 row       1 col     (character, size 7)

f_gamma

w      1 row       2 cols    (numeric)

       2.4      5000

statement : CALL at line 2109 column 1

2110  diff = max(abs(grad - FDDGrad));

ERROR: (execution) Matrix has not been set to a value.

operation : - at line 2110 column 21

operands  : grad, FDDGrad

grad      0 row       0 col     (type ?, size 0)

FDDGrad      0 row       0 col     (type ?, size 0)

statement : ASSIGN at line 2110 column 1

2111  print grad, FDDGrad, diff;

ERROR: Matrix grad has not been set to a value.

statement : PRINT at line 2111 column 1

[Rick_SAS]

In the statement

grad = g_gamma(x);

you have not defined the matrix x.

Perhaps you mean to use

grad = g_gamma(xopt);

[babaJB]

I found some error in analytic derivative, so I  fixed and checked using Wolfram Alpha.  Then i just need to assign grad = g_gamma(p) , p = is the value that i want to check the derivatives. And it works, take a look :

                                              grad

                                      1.0960612 0.0004149

                                            FDDGrad

                                       1.096059 0.0004149

                                                  diff

                                           2.2054E-6

Thank you for your article, I do need to be more careful in reading and understanding the code

[Rick_SAS]

Great! Congratulations.

For other readers, I assume the article babaJB is refers to is

http://blogs.sas.com/content/iml/2011/10/14/hints-for-derivatives/

[babaJB]

I have a question that might be out of topic

I want to generate number according gamma distribution with mean=mu and variance= var. I tried the function: RANGAM(seed,a) but I am looking The GAMMA distribution that has two parameters.

Can I do like this  ? ( beta is the scale and alpha is shape parameter)

x=beta*rangam(seed,alpha);

gama= ;

output;

end;

[Rick_SAS]

Yes.  If X~Gamma(shape, scale), then c*X~Gamma(shape, c*scale). The SAS RAND (or RANDGAM) function sets scale=1, so you can multiply the result by the scale parameter if you want a non-unit scale.

In your formula, x will have a sample mean that is close to alpha*beta and a sample variance that is close to alpha*beta##2.

[babaJB]

What if i want to implement generating random number from gamma distribution(alpha, beta) in SAS IML, let say :

x = j(1,n1);

call beta*randgen(x,'Gamma',alpha);

This code is not working . How should I do?

[Rick_SAS]

RANDGEN is not a function: it is a subroutine that fills an entire matrix or vector with one call.

Use

x = j(1,n1);

call randgen(x,'Gamma',alpha);

x = beta*x;