Hello Dr. Rick,
I used do loop for
- I simulated data sets of Weibull (a, b) each includes 15 observations for 5000 times.
- MLEs for each sample would be calculated using nonlinear optimization Newton Raphson subroutine (call nlpnra (rc,xres,"f_weib2",x,optn,con); which should be resulted in 5000 MLEs groups. Each MLEs group includes two estimates for a and b.
The problem is time consumption with nested do loop.
You recommend to use vectorization instead of do loop. I tried a simple program that includes only one do loop in different ways to change do loop to matrix work in different ways such that:
- I simulated data sets of Weibull (a, b) each includes 15 observations for 5000 times by constructing a matrix (5000,15) such that each row represents a sample of 15 observations of Weibull (a,b).
- I used a module to calculate a matrix of which likelihood functions “f” that would be used to get a vector of MLEs for 5000 samples “xopt”.
The program did not run and I got errors.
The question is how I can fix this program to get MLEs for 100 samples using vectorization instead of do loop.
The program is
proc iml;
*options nonotes;
a=12.1543;
b=4;
start f_weib2(x) global(g_x,thet); /* x[1]=b and x[2]=a and g_x is global variable, which can be x1 or x2*/
f = p*log(x[2]) - p*x[2]*log(x[1]) + (x[2]-1)*sum1 - sum2;
return(f);
finish;
count=0;
N = 15; /* number of obs in each sample */
NumSamples = 5000; /* number of samples */
call randseed(12345);/* set seed for random number stream */
x1 = j(NumSamples,N ); /* each row is sample of size N */
call randgen(x1, 'weibull', a, b); /* fill it*/
g_x = x1;
*print g_x;
n = 2; thet = 0.;
sum1=j(5000,1);
sum2=j(5000,1);
p = ncol(g_x);
temp = g_x - thet;
logtemp=log(temp);
x = j(1,n,.5); /* initial guess for solution (a, b) */
tempdev=(temp / x[1])##x[2];
sum1=logtemp[,+];/*a vector of sum over columns of the matrix logtemp to be used in LHF "f"*/
sum2=tempdev[,+];/*a vector of sum over columns of the matrix tempdev to be used in LHF "f"*/
*The following statements set up the initial conditions, options, and constraints, and they call the optimization
routine;
optn = {1 2}; /* 1: find max of a function and 2: print information*/
con = { 1.e-6 1.e-6 ,
. . }; /*lower bounds of parms >0 and upper bounds of parms < infinity (.)*/
call nlpnra(rc,xres,"f_weib2",x,optn,con);/*nonlinear optimization Newton Raphson subroutine*/
/*--- Save result in xopt, fopt ---*/
xopt = xres`; fopt = f_weib2(xopt);
quit;
run;
The log of the program is
NOTE: AUTOEXEC processing completed.
1 proc iml;
NOTE: IML Ready
2 *options nonotes;
3 a=12.1543;
4 b=4;
5 start f_weib2(x) global(g_x,thet);
5 ! /* x[1]=b and x[2]=a and g_x is global variable, which can
5 ! be x1 or x2*/
6 f = p*log(x[2]) - p*x[2]*log(x[1]) + (x[2]-1)*sum1 - sum2;
7 return(f);
8 finish;
NOTE: Module F_WEIB2 defined.
9 count=0;
10 N = 15;
10 ! /* number of obs in each sample */
11 NumSamples = 5000;
11 ! /* number of samples */
12 call randseed(12345);
12 ! /* set seed for random number stream */
13 x1 = j(NumSamples,N );
13 ! /* each row is sample of size N */
14 call randgen(x1, 'weibull', a, b);
14 ! /* fill it*/
15 g_x = x1;
16 *print g_x;
17 n = 2;
17 ! thet = 0.;
18 sum1=j(5000,1);
19 sum2=j(5000,1);
20 p = ncol(g_x);
21 temp = g_x - thet;
22 logtemp=log(temp);
23 x = j(1,n,.5);
23 ! /* initial guess for solution (a, b) */
24 tempdev=(temp / x[1])##x[2];
25 sum1=logtemp[,+];
25 ! /*a vector of sum over columns of the matrix logtemp to be used in LHF "f"*/
26 sum2=tempdev[,+];
26 ! /*a vector of sum over columns of the matrix tempdev to be used in LHF "f"*/
27 *The following statements set up the initial conditions, options, and constraints, and they
27 ! call the optimization
28 routine;
29 optn = {1 2};
29 ! /* 1: find max of a function and 2: print information*/
30 con = { 1.e-6 1.e-6 ,
31 . . };
31 ! /*lower bounds of parms >0 and upper bounds of parms < infinity
31 ! (.)*/
32 call nlpnra(rc,xres,"f_weib2",x,optn,con);
ERROR: (execution) Matrix has not been set to a value.
operation : * at line 6 column 6
operands : P, _TEM1002
P 0 row 0 col (type ?, size 0)
_TEM1002 1 row 1 col (numeric)
-0.693147
statement : ASSIGN at line 6 column 1
traceback : module F_WEIB2 at line 6 column 1
ERROR: (execution) Matrix has not been set to a value.
operation : NLPNRA at line 32 column 1
operands : *LIT1025, X, OPTN, CON
*LIT1025 1 row 1 col (character, size 7)
f_weib2
X 1 row 2 cols (numeric)
0.5 0.5
OPTN 1 row 2 cols (numeric)
1 2
CON 2 rows 2 cols (numeric)
1E-6 1E-6
. .
statement : CALL at line 32 column 1
32 ! /*nonlinear optimization Newton Raphson
32 ! subroutine*/
33 /*--- Save result in xopt, fopt ---*/
34 xopt = xres`;
ERROR: (execution) Matrix has not been set to a value.
operation : ` at line 34 column 12
operands : XRES
XRES 0 row 0 col (type ?, size 0)
statement : ASSIGN at line 34 column 1
34 ! fopt = f_weib2(xopt);
ERROR: (execution) Matrix has not been set to a value.
operation : F_WEIB2 at line 34 column 29
operands : XOPT
XOPT 0 row 0 col (type ?, size 0)
statement : ASSIGN at line 34 column 15
35 quit;
NOTE: Exiting IML.
NOTE: 5 workspace compresses.
NOTE: The SAS System stopped processing this step because of errors.
NOTE: PROCEDURE IML used (Total process time):
real time 0.10 seconds
cpu time 0.10 seconds
36 run;
