Problem IML and loop for estimation

[willi_m]

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;
read all var{grade} into Y ; 
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.

[WarrenKuhfeld]

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.

[Rick_SAS]

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"

[willi_m]

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.

[Rick_SAS · Accepted solution]

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)`;

[willi_m]

Thanks for your response, the problems was this. The X*B  hit 1 so the log equal 0 and compute was impossible.
Thank you for your help too., we won the first prize at school with this.