single and double integration of a multivariable function in SAS

[hrmannan]

Dear all,

I have been trying to singly and doubly integrate multivariable functions using the quad routine in SAS IML but getting errors. the integration domain is finite. The result should be a function of rho rather than a constant. Here is the log file for double integration. i am concerned about the error 'the matrix has not been set to a value'.  i dont know how to assign the optional command for mean and std dev of the vector as suggested. both x and y are bivariate N(0,1) variates, r is the weighted polychoric correlation and is unknown. i can increase the integration domain length for convergence if needed. please help to debugg the errors.

0 proc iml;/*double integration for weighted polychoric correlation*/
NOTE: IML Ready
101 start inner(y) global(a);
102 return(a#exp(-.5*y#y+rho#x#y/(1-rho**2))) ;
103 finish;
NOTE: Module INNER defined.
104 start outer(x) global(a);
105 a =exp(-.5*x#x/(1-rho**2)) ;
106 yinterval = 1.5 || 3;
107 /** evaluate inner integral for the function in x **/
108 w=exp(-.5*x#x/(1-rho**2));
109 call quad(w, "inner", yinterval);
110 return (w);
111 finish;
NOTE: Module OUTER defined.
112 print w;
ERROR: Matrix w has not been set to a value.

statement : PRINT at line 112 column 1
113 xinterval= {1 2.6};
114 call quad(v, "outer", xinterval);
ERROR: (execution) Matrix has not been set to a value.

operation : ** at line 105 column 1
operands : rho, *LIT1006

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

*LIT1006 1 row 1 col (numeric)

2

statement : ASSIGN at line 105 column 1
traceback : module OUTER at line 105 column 1

Convergence could not be attained over the subinterval
(1 , 2.6 )

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 114 column 1
operands : *LIT1014, xinterval

*LIT1014 1 row 1 col (character, size 5)

outer

xinterval 1 row 2 cols (numeric)

1 2.6

statement : CALL at line 114 column 1
114! /** outer integral **/
115 print v;
ERROR: Matrix v has not been set to a value. statement : PRINT at line 115 column 1

here's the log for the single integration. the answer should be a function of rho and z. y and z are bivariate standard normal variates.

6 proc iml;
NOTE: IML Ready
117 start inner(y) global(a);
118 return(exp(-0.5*(y-rho#z)##2/(1-rho##2)));
119 finish;
NOTE: Module INNER defined.
120
121 /** evaluate integral for the parameter value, a=3 **/
122 a = 3;
123 yinterval = 1.5 || a;
124 call quad(w, "inner", yinterval);
ERROR: (execution) Matrix has not been set to a value.

operation : # at line 118 column 26
operands : rho, z

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

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

statement : RETURN at line 118 column 4
traceback : module INNER at line 118 column 4

Convergence could not be attained over the subinterval
(1.5 , 3 )

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 124 column 1
operands : *LIT1007, yinterval

*LIT1007 1 row 1 col (character, size 5)

inner

yinterval 1 row 2 cols (numeric)

1.5 3

statement : CALL at line 124 column 1
125 print w;
ERROR: Matrix w has not been set to a value.

statement : PRINT at line 125 column 1

[Rick_SAS]

> ERROR: (execution) Matrix has not been set to a value.

> rho 0 row 0 col (type ?, size 0)

In both integrands, you have not defined the RHO variable. This information is in the log, so please read the error messages. Add RHO to the GLOBAL statement and give RHO a value somewhere in the program. 

[sbxkoenk · Accepted solution]

@hrmannan 

Also ... next time you put extracts from the LOG-window in your Communities post ...

Click the "Insert Code" icon ( </> ) in the header bar and copy / paste your log in the pop-up window!

That makes your post much more readable!

Cheers,

Koen

[hrmannan]

Dear Dr Rick
Rho is unknown in my case because I want to estimate it following the double integration. Can quad still be used to find the function which I need to maximise?

[Rick_SAS]

Yes, I know that RHO is not known. But optimization works by evaluating the objective function at values of RHO, then iterating towards the value that optimizes the objective. To do that, you must be able to evaluate the function for any value of RHO. 

For an example of optimizing a function that requires evaluating an integral, see Optimizing a function that evaluates an integral - The DO Loop (sas.com)

[hrmannan]

Dear Dr Rick

as you suggested i gave a starting value for the rho matrix and included rho in the global statement. still got errors stating that starting values for w and v were not given. So i inserted starting values for w and v. still got errors. I am perplexed. Note that w and v should be functions not constants.

here is the log:

249 proc iml;/*double integration for weighted polychoric correlation*/
NOTE: IML Ready
250 start inner(y) global(rho);
251 return(exp(-.5*x#x/(1-rho**2))#exp(-.5*y#y+rho#x#y/(1-rho**2))) ;
252 finish;
NOTE: Module INNER defined.
253 start outer(x) global(rho);
254 rho=j(1,2,.5);
255 w=j(1,2,.);
256 v=j(1,2,.);
257 *w =exp(-.5*x#x/(1-rho**2)) ;
258 yinterval = 1.5 || 3;
259 /** evaluate inner integral for the function in x **/
260 *w=exp(-.5*x#x/(1-rho**2));
261 call quad(w, "inner", yinterval);
262 return (w);
263 finish;
NOTE: Module OUTER defined.
264 print w;
ERROR: Matrix w has not been set to a value.

statement : PRINT at line 264 column 1
265 xinterval= {1 2.6};
266 call quad(v, "outer", xinterval);
ERROR: (execution) Matrix has not been set to a value.

operation : * at line 251 column 18
operands : _TEM1001, x

_TEM1001 1 row 1 col (numeric)

-0.5

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

statement : RETURN at line 251 column 4
traceback : module INNER at line 251 column 4
module OUTER at line 253 column 2

Convergence could not be attained over the subinterval
(1.5 , 3 )

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 261 column 4
operands : *LIT1018, yinterval

*LIT1018 1 row 1 col (character, size 5)

inner

yinterval 1 row 2 cols (numeric)

1.5 3

statement : CALL at line 261 column 4
traceback : module OUTER at line 261 column 4

Convergence could not be attained over the subinterval
(1 , 2.6 )

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 266 column 1
operands : *LIT1020, xinterval

*LIT1020 1 row 1 col (character, size 5)

outer

xinterval 1 row 2 cols (numeric)

1 2.6

statement : CALL at line 266 column 1
266! /** outer integral **/
267 print v;
ERROR: Matrix v has not been set to a value.

statement : PRINT at line 267 column 1

[sbxkoenk]

You marked my post as the solution (a while ago), even though it actually has nothing to do with the original question on your part. That makes it all the stranger that you still don't follow the recommendation I put in there. 😉

Namely:

While editing your post / reply ... Click the "Insert Code" icon ( </> ) in the header bar and copy / paste your log in the pop-up window!

That makes your LOG-messages much more readable for us!

It's not that problematic of course, but try to think about it next time.

BR,

Koen

[Rick_SAS]

I think you need to study some of the online examples of using the QUAD function. RHO should be defined outside of all functions and x in the outer function must be provided as a global variable to the inner function. In your program, you have defined RHO to be a (1 x 2) vector. It should be a scalar. You have not provided the x value to the inner integral.

I could show you how to do this, but I think your program has a lot of other problems that need to be addressed. I think you need to carefully think about what you are trying to do. Here are a few things that seem strange to me:

1) An MLE estimate uses data to fit the parameters of the model. You have no data.

2) You are integrating the bivariate normal PDF over the region [1, 2.5] x [1.5, 3]. I don't think this integral has a maximum value as a function of rho. I think the integral approaches infinity as rho -> 1.

3) You are integrating the density (PDF), not the log-PDF. Yet, you keep talking about MLE. In computational statistics, we use the log-likelihood because the LOG-PDF is better behaved, numerically.

Please use the SAS Support Communities to discuss programming problems. I see that you posted the same questions to my blog, but the blogging software is not very good for discussing programming issues or for having back-and-forth discussions.

[hrmannan]

dear Dr Rick

my intention is to find the result of the multivariable double integral which should be a function of rho so that i can maximize it in order to get a MLE of rho. i have followed your advise and have given a starting value for rho before the inner and outer integrals and have given a starting value for x.  i only get the result for outer integral v which is a function of x but not for the multivariable function w of x and y. here is the sas log:

249 proc iml;/*double integration for weighted polychoric correlation*/
NOTE: IML Ready
250 rho={.5};
251 x={1,1,.6};
252 start inner(y) global(rho);
253 return(a#exp(-.5#y#y+rho#x#y/(1-rho##2))) ;
254 finish;
NOTE: Module INNER defined.
255 rho={.5};
256 w={.};
257
258 start outer(x) global(rho);
259 *w =exp(-.5#x#x/(1-rho##2)) ;
260 a=exp(-.5#x#x/(1-rho##2));
261 yinterval = .5 || .9;
262 /** evaluate inner integral for the function in x **/
263 *w=exp(-.5*x#x/(1-rho**2));
264 call quad(w, "inner", yinterval);
265 return (w);
266 finish;
NOTE: Module OUTER defined.
267 print w;
268 v={.};
269 xinterval= {.5 .95};
270 v=exp(-.5#x#x/(1-rho##2));
271 call quad(v, "outer", xinterval);
ERROR: (execution) Matrix has not been set to a value.

operation : # at line 253 column 28
operands : rho, x

rho 1 row 1 col (numeric)

0.5

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

statement : RETURN at line 253 column 4
traceback : module INNER at line 253 column 4
module OUTER at line 258 column 2

Convergence could not be attained over the subinterval
(0.5 , 0.9 )

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 264 column 4
operands : *LIT1013, yinterval

*LIT1013 1 row 1 col (character, size 5)

inner

yinterval 1 row 2 cols (numeric)

0.5 0.9

statement : CALL at line 264 column 4
traceback : module OUTER at line 264 column 4

Convergence could not be attained over the subinterval
(0.5 , 0.95 )

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 271 column 1
operands : *LIT1019, xinterval

*LIT1019 1 row 1 col (character, size 5)

outer

xinterval 1 row 2 cols (numeric)

0.5 0.95

statement : CALL at line 271 column 1
271! /** outer integral **/
272 print v;

[hrmannan]

in continuation note that i just sent the log file showing the errors. please note that i am still in the first part of my objective which is to find the outcome of the double integral which should be a function of rho. i havent tried to maximise it yet because i havent been able to find the function which i need to maximise.

Haider

[Rick_SAS]

Currently, this problem is marked as "Solved". If your problem is not solved, you might want to change the solution status so that more experts will read the thread.

[hrmannan]

There are two posts on the topic. The problem isn’t solved so please respond with a solution
Haider Mannan

[Rick_SAS]

This program is based on your second example. the program shows the correct way to evaluate a double integral for the function that you have defined and over the range of integration.

/* double integration for weighted polychoric correlation */
proc iml;
start inner(y) global(g_x, rho);
   x = g_x;
   return(exp(-.5*x#x/(1-rho**2))#exp(-.5*y#y+rho#x#y/(1-rho**2))) ;
finish;

/* TEST THE INNER MODULE */
rho=0.5;                 /* global value of RHO */
g_x = 0;                 /* for testing: evaluate inner integral when x=g_x */
yinterval = {1.5  3};
call quad(z, "inner", yinterval);
print z;
free g_x;                /* undefine g_x; it'll be set in the OUTER module */

/* g_x is the global variable that holds the value of x for the outer integral.
   This is a constant when we evaluate the inner integral. */
start outer(x) global(g_x, rho);  
   g_x = x;   /* save copy of x so the inner integral can use it */
   yinterval = 1.5 || 3;
   /* evaluate w = inner integral when x is held constant */
   call quad(w, "inner", yinterval);
   return (w);
finish;

/* TEST THE OUTER MODULE */
xinterval= {1 2.6};
call quad(v, "outer", xinterval);
print v;

[hrmannan]

Dear Dr Rick

thank you very much. it produced a constant value for the integrand. however, i need the integrand as a function of rho in order to maximise it using my data weights (survey). the log-likelihood is-

sum(weights*log(integrand))

i tried with nlpnra subroutine but dont know how to access the data. kindly let me know.

Haider