Resources

A looping problem

0

Dear Gaussians,

I would like to estimate a 10-separate probit models within a simple loop. Staying with the same structure, I am getting results from Proc Iml in SAS, but Gauss terminates my program and goes off. Here is the programs in Gauss and SAS with proc Iml;

iter=10;
bstart= zeros(24,iter);
bprbt = zeros(24,iter);
for iter(1,iter,1);
bb            = invpd(x’x)*x’y[.,iter];
bstart[.,iter]= bb;
{b0,fmean,grad,covp,retcode} = maxlik( “mydata”,vars,&loglik,bstart[.,iter] );
call maxprt(b0,fmean,grad,covp,retcode);
proc loglik(betas,xx);
local y,x,xb,ll;
y = xx[.,iter];
x = xx[.,11:nyx];
xb= x*betas[1:nparam];
ll= y.*lncdfn(xb) + (1-y).*ln(1-cdfn(xb));
retp(ll);
endp;
bprbt[.,iter]=b0;
endfor;

SAS proc Iml;

prbest = j(kx,kd,0);

bstart = j(kx,kd,0);

* start doing loop;

do j = 1to kd;

dd = d[,j];

bols = inv(x`*x)*x`*dd;

bstart[,j] = bols;

start LogLik(param) global (dd,x);

k = ncol(x);

n = nrow(x);

xb = x*param[1:k];

f = sum ( (dd=1) # log( probnorm(xb) ) + (dd=0) # log( probnorm(-xb) ) );

return ( f );

finish;

opt = {3,

2};

call nlpnra(rc, beta, “LogLik”, bstart[,j], opt);

prbest[,j] = beta`;

call nlpfdd(f,g,h,“LogLik”,beta);

end;

I would greatly appreciate for any valuable comments and suggestions.

with best wishes.

asked June 1, 2014

5 Answers

0
  1. What version of GAUSS are you running? (i.e. GAUSS 10 for Mac or GAUSS 14 for Windows 64-bit, etc)
  2. Where does y come from in this line?
    bb = invpd(x’x)*x’y[.,iter];
    
  3. The variable nparam is not defined in this code snippet. Where is it defined?
  4. It will be very helpful, I think, to remove the definition of the proc loglik from inside the loop. This will make the code easier to understand.
aptech
342
0

Hi,

I am running Gauss 14 for windows 64-bit. “y” is a n by 10 matrix, 0/1(matrix), and “nparam” is the number of covariates defined previously.  The program does not work when I remove the proc loglik outside the loop and if I remove “proc”, then no need to compute the program with maxlik procedure, instead I should write a loop program with gradient and hessian, an usual newton algorithm, to compute my parameters. or as I did:

mm=seqa(1,1,10);

if mm == mm[1,.];
proc loglik1(betas,xx);
local y,x,b,xb;
y = xx[.,1];
x = xx[.,11:nyx];
b = betas[1:nparam];
xb= x*b;
retp( y.*lncdfn(xb)+(1-y).*ln(1-cdfn(xb)) );
endp;
elseif mm == mm[2,.];
proc loglik2(betas,xx);
local y,x,b,xb;
y = xx[.,2];
x = xx[.,11:nyx];
b = betas[1:nparam];
xb= x*b;
retp( y.*lncdfn(xb)+(1-y).*ln(1-cdfn(xb)) );
endp;

and so on.

0

I think your understanding of GAUSS proc’s is not quite correct. In GAUSS, procedures are defined at compile time (i.e. when GAUSS turns the text in a file into a program that can be executed). The if and for statements, etc happen at run-time.

Since the proc’s are already defined at compile time, whether they are inside of a particular if clause does not make any difference. For example:

x = 0;

if x == 0;
   proc (1) = foo(x);
      retp(x.*x);
   endp;
else;
   proc (1) = foo2(x);
      retp(2.*x);
   endp;
endif;

print foo(4);
print foo2(4);

Will return exactly the same output as:

x = 0;

if x == 0;

else;

endif;

print foo(4);
print foo2(4);

proc (1) = foo2(x);
   retp(2.*x);
endp;

proc (1) = foo(x);
   retp(x.*x);
endp;

This is because, even though GAUSS is flexible enough to let you define a procedure most anywhere, the procedure will not be changed during run-time.

I think you can accomplish what you are trying to do with an array of function pointers. I will have to look over the code more. I will post something later.

aptech
342
0

Here is an example using an array of function pointers so that you can pass in a different function for each iteration from the loop. In this toy example, the procedure callProcs would be like maxlik in your example.

//Create an array of function pointers
fn_pointer = &square|&half|&cube;

//Pass in a different element of the
//function pointer array on each iteration
for i(1, 3, 1);
   print callProcs(fn_pointer[i], 4);
endfor;

proc (0) = callProcs(local_fn_pointer, x);
    local local_fn_pointer:proc;
    print local_fn_pointer(x);
endp;

proc (1) = square(x);
   retp(x.*x);
endp;

proc (1) = cube(x);
   retp(x.*x.*x);
endp;

proc (1) = half(x);
   retp(x./2);
endp;
aptech
342
0

Dear All,

I am very grateful to all of you for your kind help and guidance. I finally come up with results after your suggestions as follows:

bstart= zeros(kz,kd);
bprbt = zeros(kz,kd);
b0      = zeros(kz,kd);
covb  = zeros(kz,kd*kz);
lls      = zeros(1,kd);
i=0;
do while i<=10;
bb         = invpd(z’z)*z’d[.,i];                               /*starting values from linear prob*/
bstart[.,i]= bb;
i = i+1;
{b0,fmean,grad,covp,retcode} = maxlik( “mydata1″,varbdz,&loglik,bstart[.,i] );
call maxprt(b0,fmean,grad,covp,retcode);
bprbt[.,i]=b0;
covb[.,(i-1).*kz+1:(i).*kz] = covp;
lls[1,i]  = nd.*fmean;
endo;
proc loglik(betas,dz);
local d,z,zb,ll;
d = dz[.,i];
z = dz[.,11:kdz];
zb = z*betas[1:kz];
ll = d.*lncdfn(zb) + (1-d).*ln(1-cdfn(zb));
retp(ll);
endp;

It produces all results but with an error indication (G0058 : Index out of range [maxlik.src]). I will also try other suggestions. Thanks.