New Contributor
Posts: 2

# Problem IML and loop for estimation

Hello,

I'm a student and for a school project, i need to compute many loglikelihood of a logit model based on Spector and Mazzeo's article, and this in order to produce a 3D plot.

So for this i created a sequence of parameters and i try to compute log likelihood for each couple of parameters and used a loop for this.

But programe refuses to execute at 28th iteration and i don't understand why.

Here is my code :

``````data spector;
input grade const tuce gpa ;
cards;
0 1 20 2.66
0 1 22 2.89
0 1 24 3.28
0 1 12 2.92
0 1 17 2.86
0 1 17 2.76
0 1 21 2.87
0 1 25 3.03
0 1 20 2.63
0 1 23 3.32
0 1 23 3.57
0 1 26 3.53
0 1 19 2.74
0 1 25 2.75
0 1 19 2.83
0 1 23 3.12
0 1 22 2.06
0 1 14 2.89
0 1 26 3.51
0 1 24 2.67
0 1 21 3.10
1 1 21 4.00
1 1 29 3.92
1 1 25 3.26
1 1 25 3.16
1 1 28 3.62
1 1 24 3.54
1 1 27 2.83
1 1 17 3.39
1 1 21 3.65
1 1 23 4.00
1 1 19 2.39
;
run;

PROC IML ;
use spector;
read all var{const tuce gpa} into X;
close spector;
print X Y ;

N = nrow(X);

*step 1  : beta sequence ;
X2=do(-1.5,1.5,0.5)`;
X3=do(1,4,0.5)`;
nX2=nrow(X2);
nX3=nrow(X3);
print nX2, nX3;

*step 2 : beta occurence;
rX2=repeat(X2,nX3,1);
rX3=repeat(X3,nX2,1);
rcst=repeat(-10.656002,nrow(rX2),1);

print rX2, rX3;

*step 3 sort;
call sort(rX3);
print rX3;

*step 4 concat;
theta=rcst||rX2||rX3;
ntheta = nrow(theta);
print theta ntheta;

*step 5 compute loglikelihood for each sequence of beta;
*Compute log likelihood of the sample;

do j=1 to 49;
XB = X*(theta`)[,j];
expXB = exp(XB);
P = expXB/(1+expXB);
logli1 = (Y#log(P))[+] ;
b=1-P;
a = 1-Y;
logli2 = (log(b##a))[+];
logli = logli1 + logli2 ;
logln = logln//logli ;
end;

print logln;``````

I know that " b##a " is not a good way but it's the only one i've found.

Someone can help me please ?

Thank you in advance for helping.

SAS Super FREQ
Posts: 508

## Re: Problem IML and loop for estimation

[ Edited ]

log(a * b) = log(a) + log(b).

log(a ## b) = b * log(a)

Avoid the multiplication or exponentiation followed by log to keep the numbers smaller and more manageable.

Log(a ## b) might overflow whereas b * log(a) will be better behaved.

Thanks, Rick.  I corrected my typos.

SAS Super FREQ
Posts: 4,274

## Re: Problem IML and loop for estimation

[ Edited ]

Warren's advice is good, although he meant to write that

log(a # b) = log(a) + log(b)

log(a ## b) = b # log(a)

Here are some articles for you to read:

1. "Two simple ways to construct a log-likelihood function in SAS"
2. "Two ways to compute maximum likelihood estimates in SAS"
3. To evaluate the log-likelihood at on a uniform grid of points, use the EXPANDGRID function as shown in "How to find an initial guess for an optimization"
4. To graph a surface plot, see "Create a surface plot in SAS"
New Contributor
Posts: 2

## Re: Problem IML and loop for estimation

[ Edited ]

Yes i read this somewhere in the forum.
But in writting my log as    " b*log(a) " there is an error of compute.

``````332      logli2 = (a*log(b))[+];
ERROR: (execution) Matrices do not conform to the operation.

operation : * at line 332 column 16
operands  : a, _TEM1001
a     32 rows      1 col     (numeric)
_TEM1001     32 rows      1 col     (numeric)

statement : ASSIGN at line 332 column 5``````

There are no error when i write  " a#log(b) "   and out of the loop

``````	XB = X*(theta`)[,1];
expXB = exp(XB);
P = expXB/(1+expXB);
logli1 = (Y#log(P))[+] ;
b=1-P;
a = 1-Y;
logli2 = (a#log(b))[+];   /* HERE */
logli = logli1 + logli2 ;
logln = logln//logli ;``````

But when i want to implemant this in a loop there is a problem

``````345      do j=1 to 49;
346      XB = X*(theta`)[,j];
346!
347      expXB = exp(XB);
347!
348      P = expXB/(1+expXB);
348!
349      logli1 = (Y#log(P))[+] ;
350      b=1-P;
351      a = 1-Y;
352      logli2 = (a#log(b))[+];
353
354      logli = logli1 + logli2 ;
355      logln = logln//logli ;
356      end;
ERROR: (execution) Invalid argument to function.

count     : number of occurrences is 2
operation : LOG at line 352 column 20
operands  : b
b     32 rows      1 col     (numeric)

statement : ASSIGN at line 352 column 5``````

Moreover, my loop is ok for a less number of iteration (13)

``````429      do j=1 to 13;             /*13 iterations    OK*/
430      XB = X*(theta`)[,j];
430!
431      expXB = exp(XB);
431!
432      P = expXB/(1+expXB);
432!
433      logli1 = (Y#log(P))[+] ;
434      b=1-P;
435      a = 1-Y;
436      logli2 = (a#log(b))[+];
437
438      logli = logli1 + logli2 ;
439      logln = logln//logli ;
440      end;``````

And with 14 iterations : problems

``````441      do j=1 to 14;                 /* 14 iterations PROBLEMS */
442      XB = X*(theta`)[,j];
442!
443      expXB = exp(XB);
443!
444      P = expXB/(1+expXB);
444!
445      logli1 = (Y#log(P))[+] ;
446      b=1-P;
447      a = 1-Y;
448      logli2 = (a#log(b))[+];
449
450      logli = logli1 + logli2 ;
451      logln = logln//logli ;
452      end;
ERROR: (execution) Invalid argument to function.

count     : number of occurrences is 2
operation : LOG at line 448 column 20
operands  : b
b     32 rows      1 col     (numeric)

statement : ASSIGN at line 448 column 5``````

So i don't understand where is the problem, and this is problematic because after i would like to compute more of 160 000 loglikelihood for my 3D plot in order to obtain a smoother graph.

SAS Super FREQ
Posts: 4,274

## Re: Problem IML and loop for estimation

The error seems to indicate that you are trying to evaluate the LOG function for an invalid parameter. Note that if P=0 then LOG(P) is undefined and if P=1 then LOG(1-P) is undefined. In your case, I think you are running into the P=1 problem for certain values of j.

I don't completely understand what you are doing, but it looks like you are evaluating the logistic log-likelihood function for linear combinations of the variables for a grid of (X2, X3) coffficients. For some values of j, the value of XB is so large that logistic(XB) is numerically equal to 1.  Try using a more restricted range of coefficients, such as

X2=do(-1,1,0.5)`;

Discussion stats
• 4 replies
• 253 views
• 0 likes
• 3 in conversation