calling quad and nlpqn in proc iml

[Emara]

Hi all,

I want to write a SAS 9 optimization syntax in which the objective

function is a numerical integration subject to certain constraints

where the integrand and the constraints don't have closed forms , but

they are a result of  loops. The function I want to optimize is

the integration of FUN(s), subject to the constraint that FUN(s)=10

when s=0.5.The decesion variable is "M" and the integration variable is "s". Here is the code I wrote but it gave me errors . The lst file is shown after the syntax

Here is my syntax:

proc iml;

start FUN(s);

  TM=0;

  do i=1 to 10;

  X=0;

  rl=0;

  do until(X>M);

  X=normal(-1)+s;

  rl=rl+1;

  end;

  TM=TM//rl;

  end;

  bb=(sum(TM))/10;

  return(bb);

  finish;

start FUN2(M);
  a={0   2};
  call quad(k,"FUN",a)eps=1E-4 peak=0.01;
  return(bb1);
  finish;
 

start FUN3(M);
bb2=10;
  TM2=0;
  do J=1 to 10;
  X2=0;
  rl2=0;
  do until(X2>M);
  X2=normal(-1)+0.5;
  rl2=rl2+1;
  end;
  TM2=TM2//rl2;
  end;
  bb2=(sum(TM2))/10;
  return(bb2);
  finish;
 

M=3;
optn= j(1,10,.); optn[2]= 2;
CALL NLPQN(rc,Mres,"FUN2",M,optn) nlc="FUN3";
 

quit;
 

The log file is:
 

1    proc iml;
NOTE: IML Ready
2    start FUN(s);
3      TM=0;
4      do i=1 to 10;
5      X=0;
6      rl=0;
7      do until(X>M);
8      X=normal(-1)+s;
9      rl=rl+1;
10     end;
11     TM=TM//rl;
12     end;
13     bb=(sum(TM))/10;
14     return(bb);
15     finish;
NOTE: Module FUN defined.
16
17
18   start FUN2(M);
19     a={0   2};
20     call quad(k,"FUN",a)eps=1E-4 peak=0.01;
21     return(bb1);
22     finish;
NOTE: Module FUN2 defined.
23
24   start FUN3(M);
25   bb2=10;
26     TM2=0;
27     do J=1 to 10;
28     X2=0;
29     rl2=0;
30     do until(X2>M);
31     X2=normal(-1);
32     rl2=rl2+1;
33     end;
34     TM2=TM2//rl2;
35     end;
36     bb2=(sum(TM2))/10;
37     return(bb2);
38     finish;
NOTE: Module FUN3 defined.
39
40   M=3;
41   optn= j(1,10,.);
41 !                  optn[2]= 2;
42   CALL NLPQN(rc,Mres,"FUN2",M,optn) nlc="FUN3";
ERROR: (execution) Matrix has not been set to a value.
 

operation : > at line 7 column 13
operands  : X, M
 

X      1 row       1 col     (numeric)
 

0.1364511
 

M      0 row       0 col     (type ?, size 0)
 

statement : DO at line 7 column 3
traceback : module FUN at line 7 column 3
             module FUN2 at line 18 column 1
 

ERROR: Error evaluating DO UNTIL expression; DO loop exited.
Convergence could not be attained over the subinterval
             (0 , 2 )
 

The function might have one of the following:
     1) A slowly convergent or a divergent integral.
     2) Strong oscillations.
     3) A small scale in the integrand: in this case
        you can change the independent variable
        in the integrand, or,
        provide the optional vector describing roughly
        the mean and the standard deviation of the
        integrand
     4) Points of discontinuities in the interior
        in this case, you can provide a vector containing
        the points of discontinuity and try again.
 

ERROR: Execution error as noted previously. (rc=100)
 

operation : QUAD at line 20 column 3
operands  : *LIT1010, A, *LIT1011, *LIT1012
 

*LIT1010      1 row       1 col     (character, size 3)
 

FUN
 

A      1 row       2 cols    (numeric)
 

         0         2
 

*LIT1011      1 row       1 col     (numeric)
 

    0.0001
 

*LIT1012      1 row       1 col     (numeric)
 

      0.01
 

statement : CALL at line 20 column 3
traceback : module FUN2 at line 20 column 3
 

ERROR: Execution error as noted previously. (rc=100)
 

operation : NLPQN at line 42 column 1
operands  : *LIT1028, M, OPTN, , , , , , *LIT1029
 

*LIT1028      1 row       1 col     (character, size 4)
 

FUN2
 

M      1 row       1 col     (numeric)
 

         3
 

OPTN      1 row      10 cols    (numeric)
 

         .
2         .         .         .         .         .         .         .
        .
 

*LIT1029      1 row       1 col     (character, size 4)
 

FUN3
 

statement : CALL at line 42 column 1
43
44   quit;
NOTE: Exiting IML.
NOTE: The SAS System stopped processing this step because of errors.
NOTE: PROCEDURE IML used (Total process time):
      real time           0.18 seconds
      cpu time            0.17 seconds
 

Thanks,
A.Emara

[Rick_SAS]

The short answer is that M is not defined in the FUN function. From your code, I'm guessing you want to use GLOBAL(M) when defining the FUN function.

I do not understand what you are trying to accomplish, but I think there might be other problems in your code. It concerns me that the integrand (FUN) is random. It looks like you are computing the integral of a random function. That will itself be random, so I don't see how it can have a well-defined optimum.   You are using a Quasi-newton method, but this function has no derivative!

Are you trying to do a Monte Carlo simulation of something?  Are you using the QUAD function in oder to compute the CDF of some distribution?  If so, there are simpler ways to do these tasks.

[Emara]

Hi,

Thanks alot for your reply. I have some points:

1-What do you mean by Global M? Is this something  I could add to the code?

2-Is there a theoritical problem in finding an optimum for an integral with a random integrand?Can't I find a local or global optimum for such problem?

3- Can i use other optimization methods than Quasi-method that doesn't assume a derivative?

4- What I am trying to do is the following:

Minimize the integration of FUN(s,M)

subject to

FUN(0,M)=10

where

1-M is the decesion variable I want to find

2-s is the variable representing the known limits of integration

3-FUN doesn't have a closed form and is computed using Monte Carlo simulation

4- I am not trying to find the CDF of any distribution

Could I solve this problem using SAS 9?

A.Emara

[Rick_SAS]

I'll start by saying that it looks like FUN computes the average (expected) number of draws that are required before X~N(s,1) is great than M. I point this out in case anyone else want to help you make that function more efficient, which is required if you want to change 10 trials to a larger number.

1) start Fun(s) global(M);

See http://support.sas.com/documentation/cdl/en/imlug/65547/HTML/default/viewer.htm#imlug_programstateme...

2) Yes, there is a theoretical problem, but it has nothing to do with integrals.  The problem is that

Minimize X

when X is a random variable does not have meaning. For example, add the GLOBAL statment to your code and then run the following:

M = 0;
s = do(-1,1,0.1);
f = j(1,ncol(s));
do i = 1 to ncol(s);
   f = FUN(s);
end;

create Viz var{s f}; append; close Viz;

proc sgplot data=Viz;
series x=s y=f;
run;

Notice that your function does not have a well-defined minimum as a function of s.

3) No, because the problem isn't well-defined.

4) Yes, you can solve the problem you describe in SAS, but you need (A) to rethink what you are doing and how you are doing it, (B) vectorize your IML code so that it runs faster

Here are some blog posts you might find useful:

http://blogs.sas.com/content/iml/2010/09/22/efficient-sampling/

Shows how to generate random numbers without using loops in SAS/IML. This is essential.

Pre-allocate arrays to improve efficiency - The DO Loop

If you can help it, avoid concatenating values in a loop.

Simulation of Buffon’s needle in SAS - The DO Loop This is equivalent to a Monte-Carlo integration, although it might not be obvious

Simulating the Coupon Collector's Problem - The DO Loop

A time-to-event simulation. This seems very similar to your problem.

A surprising result: The expected number of uniform variates whose sum exceeds one - The DO Loop

Another time-to-event simulation that computes the expected value of some quantity. Sometimes a second example is useful.

[Emara]

Hello........

Sorry for the trouble, but I still have some problems,

1-Adding Global(M) removed some of the errors, but didn't remove this one

"41 ! optn[2]= 2;

42 CALL NLPQN(rc,Mres,"FUN2",M,optn) nlc="FUN3";

ERROR: Execution error as noted previously. (rc=100)

operation : NLPQN at line 42 column 1

operands : *LIT1028, M, OPTN, , , , , , *LIT1029

*LIT1028 1 row 1 col (character, size 4

)"
 

2- Sorry, but could you clarify more the theoritical problem in the code. Do you mean that , in general, minimizing a random function has no meaning?And since FUN(s) is a resuslt of Monte carlo simulation then it's minimization is of no meaning and thus this problem couldn't be solved properly using any optimization algorithm?i.e, it is not a matter of a code or of the programming language itself?

3-I tried to run the code you send but the procedure"sgplot" is not defined in SAS9

4- I Want to make sure of something in my code, is the constraint that FUN3=10 added correctly to my model? Is its module correct and has NLPQN recognized it as a constraint ?

Thanks,

A.Emara

[Rick_SAS]

1) Your nlc function is not (that I can see) returning constraints. I don't know what FUN3 is doing.

2) Now that you've explained what you are trying to do, I take back my earlier reservations. I was confused because I wasn't sure what you were trying to accomplish. I now see that your objective function is a simulation that computes an estimate of some quantity that depends on M.

3) Use GPLOT.

4) Run the following code. Are you expecting to get 10 for all values of M?

do M=0 to 3 by 0.1;

  Y=Fun3(M);

  print M Y;

end;

5) I was thinking about your problem over the weekend. Are you trying to find the distribution of first-passage times for Brownian motion with zero drift and a given mean?  If so, this is a limiting case of the inverse Gaussian distribution. The limit of zero drift is known as the Levy distribution. You can simulate the distribution of times directly, which would save you lots of time and effort.

See Inverse Gaussian distribution - Wikipedia, the free encyclopedia

and Lévy distribution - Wikipedia, the free encyclopedia

The integral (CDF) of these functions is also known exactly and can be computed from the normal CDF.

[Emara]

Hello..........

1- FUN3  represents the constraint of my model. FUN3 has the same formula as FUN "the objective function" but at s=0. Thus my problem is to  find the optimal "M" value that minimizes the integration of FUN subject to the constraint that FUN3=10. Thus my optimal M should give me the minimum integration value of FUN but at the same time FUN3 should be equal to 10 at the optimal M. That is I am restricting the value of FUN to be equal to 10 at s=0

**I should note that i have a mistake in my first post as I should have replaced   "X2=normal(-1)+0.5;"    by     "X2=normal(-1);"    in FUN3

4)When i run your code  , I didn't get a 10 but I got different values of Y for different M values. What i need is to get a 10 at the optimal M since FUN3 represents my constraint.

5) No, I am not trying to do so.

Thanks alot,

A.Emara