Parameter estimation (MLE)

[seyibote]

Trying to estimate the parameters for the following log-likelihood, have the following code but sas keeps giving me the error "ERROR: (execution) Matrix has not been set to a value." here is my code; if anyone can help will appreciate it attached is also the function I am working with

proc iml;

aarset=

{0.1 0.2 1 1 1 1 1 2 3 6 7 11 12 18 18 18 18 18

21 32 36 40 45 46 47 50 55 60 63 63 67 67 67 67

72 75 79 82 82 83 84 84 84 85 85 85 85 85 86 86};

P = ncol(aarset); m = p;

start f_gglfr(x) global(aarset);

/* use x-a, x[2]-b, x[3]-d, x[4]-delta */

p=ncol(aarset); m=p-4;

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

      do i = 1 to p;

         temp = (x[1]+(x[2]/2)*aarset##2);

         if i <= m then sum1 = sum1 + (x[1]+x[2]*aarset);

         sum2 = sum2 + (-x[3]*log(1-exp(-temp)));

  sum3 = sum3 + (x[1]*aarset+(x[2]/2)*aarset##2);

  sum4 = sum4 + (1-exp(-temp));

      end;

      f = m*log(x[3]) - m*log(Gamma(x[4])) + log(sum1) +(x[4]-1)*log(sum2) - sum3 + (x[3]-1)*log(sum4);

      return(f);

   finish f_gglfr;

   start g_gglfr(x) global(aarset);

     /* use x-a, x[2]-b, x[3]-d, x[4]-delta */

      p=ncol(aarset); m=p-4;

      g = j(1,4,0.);

      sum1 = 0.; sum2 = 0.; sum3 = 0.; sum4 = 0.; sum5 = 0.; sum6 = 0.; sum7 = 0.; sum8 = 0.; sum9 = 0.; sum10 = 0.; sum11 = 0.;

      do i = 1 to p;

         temp = (x[1]+(x[2]/2)*aarset##2);

         if i <= m then sum1 = sum1 + (x[1]+x[2]*aarset);

         sum2 = sum2 + (aarset*exp(-temp))/((1-exp(-temp)*log(1-exp(-temp)));

  sum3 = sum3 + (aarset*exp(temp))/(1-exp(temp));

  sum4 = sum4 + aarset;

  sum5 = sum5 + aarset/(x[1]+x[2]*aarset);

  sum6 = sum6 + (aarset##2*exp(-temp))/(log(1-exp(-temp)*(1-exp(-temp));

  sum7 = sum7 + aarset##2/2;

  sum8 = sum8 + (aarset##2*exp(temp))/(1-exp(temp));

  sum9 = sum9 + (log(1-exp(-temp)))/(log(1-exp(-temp)));

  sum10 = sum10 + (log(1-exp(-temp)));

  sum11 = sum11 + log(-x[3]*log(1-exp(-temp)));

      end;

  temp = (1/2)*(2*x[1]+x[2]*aarset##2);

      g[1] = 1/sum1+(x[4]-1)*sum2+sum4+(x[3]-1)*sum3;

      g[2] = sum5+(x[4]-1)/2*sum6+sum7+(x[3]-1)/2*sum8;

      g[3] = m/x[3] - (x[4]-1)/x[3]*sum9+sum10;

      g[4] = -m*psi(x[4])/Gamma(x[4])+sum11;

      return(g);

   finish g_gglfr;

   n = 4;

   x0 = j(1,n,2);

   optn = {1 2};

   con = { 0 0 0 0,

           . . . . };

   call nlptr(rc,xres,"f_gglfr",x0,optn,con,,,,"g_gglfr");

   /*--- Save result in xopt, fopt ---*/

   xopt = xres`; fopt = f_gglfr(xopt);

quit;

[Rick_SAS]

You should always call your function before you try to optimize. The call gives an error:

x = {2 2 2 2};

f = f_gglfr(x);

The reason is that your TEMP variable has some large numbers (greater than 7000), so 1-exp(-temp) is numerically 0.

I think there might be another problem.When aarset is a vector, the function f_ggfr returns a vector.  I assume you intend for f_ggfr to be a scalar function.