Hutch_sas Tracker

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Fri, 25 Jul 2014 15:43:54 GMT

The optimization method should only be generating feasible points in the search. It looks like you are setting conh before you define mincvf(); your statement:

conh={0.4 . .,. . .,. . .};

should come AFTER the"start mincvf(CVF1);" statement.

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Thu, 24 Jul 2014 20:22:38 GMT

You could specify a lower bound or linear constraint, see in the IML doc under Nonlinear Optimization Examples->Details->Parameter Constraints

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Wed, 23 Jul 2014 19:55:16 GMT

the warning appears to be saying that the value of (sum(le)-le) is 0. This could be due to loss

of precision if le >> le, i != x. I recommend you add some conditional print statements, like

if (sum(le) - le) = 0 then do;

print temp1, x, i, ue, le;

end;

in front of the statement: t5=(sum(ue)-ue)/(sum(le)-le);

to see what is going on.

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Wed, 23 Jul 2014 14:50:55 GMT

start outer(p,m) global(newp);

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Wed, 23 Jul 2014 13:28:19 GMT

Your subroutine "outer" needs to declare newp in a global clause. Otherwise, when it sets newp it is working with a local variable only, and the "global" newp is not initialized. That is causing the first error in your log.

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Tue, 22 Jul 2014 14:36:34 GMT

Oops, I see that you will also have to manipulate m within the loop, but you get the idea?

Re: Is it possible to redefine subroutines in a subroutine?

Hutch_sas — Tue, 22 Jul 2014 14:30:30 GMT

In IML all subroutine definitions are global, even if you define then inside another subroutine. In your example, could you not accomplish the same goal by using a global clause in mhatx2 and Kmod? That is, use

start mhatx2(h) global(p,m,pi,e);

start Kmod(x,h) global(pi,e);

and manipulate p in your "outer" routine. maybe use something like:

start outer(p_input) global(p,e,pi,m); m=nrow(p_input);

[your existing code]

p = p_input(w);

hyi=mincvf(CVF1);

...

Re: Product of functions/Variable scope

Hutch_sas — Tue, 15 Jul 2014 14:04:06 GMT

You need to modify the line:

call quad(val, "Func", yinterval); /* evaluate integral on [0, 1] */

val will not be replaced by the quad call. do

call quad(result, "Func", yinterval); /* evaluate integral on [0, 1] */

val = result;

Re: How to aggregate a tree efficiently?

Hutch_sas — Tue, 15 Jul 2014 13:39:45 GMT

What operations do you need to do on the tree? Do you need to be able to add/delete nodes or prune branches? How is the tree to be constructed?

Re: compute a formula

Hutch_sas — Fri, 11 Jul 2014 16:52:32 GMT

For a 4th degree poynomial there will be generally 4 complex solutions. Looks like the polyroot() function is what you need, check the IML documentation. I suggest you first compute z = polyroot( {1.2 1.2 1.2 1.2 -10) to get the imaginary roots, then since z = 1 + x, x = z-1, just subtract 1 from the first column of z to get x, i. e.

proc iml;

a = {1.2 1.2 1.2 1.2 -10};

z = polyroot(a);

x = z - {1 0,1 0,1 0,1 0};

print x;

Re: Computational Accuracy and Efficiency for the NLPNMS

Hutch_sas — Tue, 01 Jul 2014 14:33:15 GMT

The other default values seem to be in line with the 9.1 doc, so there could be an error in the documentation regarding that parameter. It might be helpful if you shared your program, to verify your setup for the optimization.

Re: Computational Accuracy and Efficiency for the NLPNMS

Hutch_sas — Mon, 30 Jun 2014 13:42:26 GMT

There are many different ways tolook at computational efficiency. The Nelder-Mead simplex algorithm, because it does not use gradients directly, may require more objective function evaluations than some other techniques such as the conguate gradient method. However, the computation of gradients and Hessians can sometimes be very cpu-intensive themselves, so the least cpu-intensive method generally depends on the problem constraints and objective function.

There are also several ways to estimate the accuracy of the results. Look at the SAS/IML doc under Nonlinear Optimization Examples->Details->Termination Criteria. IML provides several options for terminating the optimization process for NLPNMS, based on the relative or absolute change in value of either the objective function or the parameters between iterations of the algorithm, and that should give you some estimate of the accuracy of the results. Obviously, the tighter the termination criteria, the more iterations will be required. Since your input data or model iteslf is probably not exact in the first place, requiring extreme accuracy may not be justified and can lead to excessive run times or numerical problems. Keep in mind also that all the nlp methods will converge to a local optimum, so unless your objective function is complex it is possible that your answer might not always be the global one.

Re: using nlpcg in sas iml

Hutch_sas — Thu, 12 Jun 2014 12:35:57 GMT

Could you post your current program, and your run-time log?

Re: using nlpcg in sas iml

Hutch_sas — Tue, 10 Jun 2014 20:15:44 GMT

The problem is that when IML calls your objective function during the optimization, it is supplying a row vector for x. Try transposing x in the second line of your objective module, i.e.

sum= m [,2:3]*x` + m [,1];

and see if that doesn't work as it is supposed to.

Re: maximization in proc iml with global variable

Hutch_sas — Wed, 28 May 2014 12:09:40 GMT

Could you give us or point us to a complete mathematical description of your problem?

Re: maximization in proc iml with global variable

Hutch_sas — Tue, 27 May 2014 16:01:11 GMT

You did not define r. The objective function uses r, so if it is not defined it will give you an error.

Re: subroutine

Hutch_sas — Fri, 02 May 2014 12:21:29 GMT

proc iml;

start AddToDiag( result, A, B);

/* add B to diagonal of a square matrix A */

result = A + diag( j(nrow(A), 1, B) );

finish;

/* test it out */

m = {1 1 1,

2 2 2,

3 3 3 };

lambda = 3;

run AddToDiag( MPlusLambda, m, lambda);

print M lambda MPlusLambda;

quit;

Re: matrices

Hutch_sas — Thu, 01 May 2014 19:26:13 GMT

A quick way to add a constant to the diagonal elements of a matrix:

/* add B to diagonal of a square matrix A */

C = A + diag( j(nrow(A), 1, B) );

Re: A Problem with NLPNMS when using Modules with random generations

Hutch_sas — Mon, 21 Apr 2014 13:38:57 GMT

Oops, I missed the fact that you were still at 9.0. RANDGEN and RANDSEED were introduced in 9.1. I don't believe the old random number generation functions allow the arbitrary re-initialization of the random number stream. In order to get the same sampling, you will have to create a table of random variates at the beginning of your program, and re-use that stream at each call. If you need more samples than can be stored in memory, you might have to store them in a data set and read them in as called for.

Re: A Problem with NLPNMS when using Modules with random generations

Hutch_sas — Sun, 20 Apr 2014 16:30:39 GMT

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;