After you posted your code, I tried running your objective function a few times. As one would expect, it gives slightly different answers if you run it with the same input multiple times, because the random number stream is different each time the objective function or non-linear constraint function is evaluated. This can be easily fixed. Instead of calling the NORMAL function each time for your objective and constraint functions, use the newer RANDGEN call. You can then re-set the random number seed to the same value by calling RANDSEED at the beginning of each objective/constraint module call, guaranteeing that you use the same random number stream each time. I also noticed that you could hoist a couple of operations outside your inner loop to marginally improve performance. I made these changes in your program, and it ran in about 3 hrs on my PC. I pasted it in below. Give it a try and see if you get better results. I also recommend that you try some different seeds in the RANDSEED call to see how sensitive your results are to the choice of random seed. proc iml; start ARL(hr); call randseed(1000,1); shift=0.5; * p is number of variables; p=2; * nruns is the number of simulation runs; nruns=100000; arlv=0; sh=shift//j(p-1,1,0); /* sh unchanged throughout the simulation */ x = j(p,1); /* set up x with proper dimensions, just once */ * starting the simulation; do reps=1 to nruns; rl=0; *initialy z is set as a px1 vector of zeros; z=0//j(p-1,1,0); *ti carries the sum of squares of z, i.e, each element of z is squared and then we take the sum; ti=ssq(z); do until (ti>hr[1]); *generating p univariate standard normal variables; call randgen(x, 'normal'); *adding the shift for the out of control case; x=x+sh; *computing the norm of the error; ei=sqrt(ssq(x-z)); if ei>hr[3] then w=(ei-(1-hr[2])*hr[3])/ei; else w=hr[2]; z=(1-w)*z+w*x; ti=ssq(z); rl=rl+1; end; arlv=arlv+rl; end; arlr=arlv/nruns; bb=arlr; return(bb); finish; start ARL0(hry)global(hr); hr=hry; call randseed(1000,1); shift=0; bb0=j(1,1,0.); p=2; sh=shift//j(p-1,1,0); /* sh unchanged throughout the simulation */ x = j(p,1); /* set up x with propoer dimensions, just once */ * nruns is the number of simulation runs; nruns=100000; arlv=0; * starting the simulation; do reps=1 to nruns; rl=0; *initialy z is set as a px1 vector Of zeros; z=0//j(p-1,1,0); *ti carries the sum of squares of Z, i.e, each element of z is squared and then we take the sum; ti=ssq(z); do until (ti>hr[1]); *generating p univariate standard normal variables; call randgen(x,'normal'); *adding the shift for the out of control case; x=x+sh; *computing the norm of the error; ei=sqrt(ssq(x-z)); if ei>hr[3] then w=(ei-(1-hr[2])*hr[3])/ei; else w=hr[2]; z=(1-w)*z+w*x; ti=ssq(z); rl=rl+1; end; arlv=arlv+rl; end; arlr=arlv/nruns; bb0[1]=(arlr)-200; return(bb0); finish; con = { 0 0 0, . 1 . }; *initial vector; hr={0.8856 0.17 3.85}; optn= j(1,11,.);optn[2]=5;optn[10]=1;optn[11]=0; CALL NLPNMS(rc,hrres,"ARL",hr) blc=con opt=optn nlc="ARL0"; print hrres; quit;
... View more