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 = **1**to 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.

## 5 Answers

- What version of GAUSS are you running? (i.e. GAUSS 10 for Mac or GAUSS 14 for Windows 64-bit, etc)
- Where does
`y`come from in this line?bb = invpd(x’x)*x’y[.,iter];

- The variable
`nparam`is not defined in this code snippet. Where is it defined? - 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.

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.

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.

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;

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.