Chapter Contents |
Previous |
Next |

Language Reference |

**CALL QUAD(***r*, "fun", points <,**EPS=**eps><,**PEAK=**peak>

*<,***SCALE=**scale><,**MSG=**msg><,**CYCLES=**cycles>**);**

The QUAD subroutine returns the following value:

*r*- is a numeric vector containing the results of the integration.
The size of
*r*is equal to the number of subintervals defined by the argument*points*. Should the numerical integration fail on a particular subinterval, the corresponding element of*r*is set to missing.

The inputs to the QUAD are as follows:

*"fun"*- specifies the name of an IML module
used to evaluate the integrand.

*points*- specifies a sorted vector that provides the
limits of integration over connected subintervals.
The simplest form of the vector provides the
limits of the integration on one interval.
The first element of
*points*should contain the left limit. The second element should be the right limit. A missing value of`.M`in the left limit is interpreted as , and a missing value of`.P`is interpreted as . For more advanced usage of the QUAD call,*points*can contain more than two elements. The elements of the vector must be sorted in an ascending order. Each two consecutive elements in*points*defines a subinterval, and the subroutine reports the integration over each specified subinterval. The use of subintervals is important because the presence of internal points of discontinuity in the integrand will hinder the algorithm.

*eps*- is an optional scalar specifying the desired relative accuracy.
It has a default value of 1E-7.
You can specify
*eps*with the keyword EPS.

*peak*- is an optional scalar that is the approximate
location of a maximum of the integrand.
By default, it has a location of 0 for infinite
intervals, a location that is one unit away from
the finite boundary for semi-infinite intervals,
and a centered location for bounded intervals.
You can specify
*peak*with the keyword PEAK.

*scale*- is an optional scalar that is the approximate
estimate of any scale in the integrand along
the independent variable (see the examples).
It has a default value of 1.
You can specify
*scale*with the keyword SCALE.

*msg*- is an optional character scalar that restricts the
number of messages produced by the QUAD subroutine.
If
*msg*= "NO" then it does not produce any warning messages. You can specify*msg*with the keyword MSG.

*cycles*- is an optional integer indicating the number of
refinements allowed to achieve the required accuracy.
It has a default value of 8.
You can specify
*cycles*with the keyword CYCLES.

If the dimensions of any optional argument are 0 ×0, the QUAD call uses its default value. The QUAD subroutine quad is a numerical integrator based on adaptive Romberg-type integration techniques. Refer to Rice (1973), Sikorsky (1982), Sikorsky and Stenger (1984), and Stenger (1973a, 1973b, 1978). Many adaptive numerical integration methods (Ralston and Rabinowitz 1978) start at one end of the interval and proceed towards the other end, working on subintervals while locally maintaining a certain prescribed precision. This is not the case with the QUAD call. The QUAD call is an adaptive global-type integrator that produces a quick, rough estimate of the integration result and then refines the estimate until achieving the prescribed accuracy. This gives the subroutine an advantage over Gauss-Hermite and Gauss-Laguerre quadratures (Ralston and Rabinowitz 1978, Squire 1987), particularly for infinite and semi-infinite intervals, because those methods perform only a single evaluation.

Consider the integration

The following statements evaluate this integral:

proc iml; /* Define the integrand */ start fun(t); v = exp(-t); return(v); finish; /* Call QUAD */ a = { 0 .P }; call quad(z,"fun",a); print z[format=E21.14];

The integration is carried out over the interval , as specified by the variable A. Note that the missing value in the second element of A is interpreted as . The values of

The following code performs the integration over two subintervals, as specified by the variable A:

proc iml; /* Define the integrand */ start fun(t); v = exp(-t); return(v); finish; /* Call QUAD */ a = { 0 3 .P }; call quad(z,"fun",a); print z[format=E21.14];

Note that the elements of A are in ascending order. The integration is carried out over (0,3) and ,and the corresponding results are shown in the output. The values of

b = sum(z); print b[format=E21.14];

The purpose of the

Consider the integration

proc iml; /* Define the integrand */ start fun(t); v = exp(-10000*t); return(v); finish; /* Call QUAD */ a = { 0 .P }; /* Either syntax can be used */ /* call quad(z,"fun",a,1E-10,0.0001); or */ call quad(z,"fun",a) eps=1E-10 peak=0.0001 ; print z[format=E21.14];

Only one interval exists. The integration is carried out over . The default values of

If you do not specify a

Consider the integration

proc iml; /* Define the integrand */ start fun(t); v = exp(-100000*(t-3)*(t-3)); return(v); finish; /* Call QUAD */ a = { .M .P }; call quad(z,"fun",a) eps=1E-10 peak=3 scale=0.001 ; print z[format=E21.14];

Only one interval exists. The integration is carried out over .The default value of

If you use the default value of

When you perform double integration, you must separate the variables between the iterated integrals. There should be a clear distinction between the variables of the one-dimensional integration at hand and the parameters to be passed to the integrand. Posting the correct limits of integration is also an important issue. For example, consider the binormal probability, given by

You can write the equation in the form of a function with the parameters as arguments. The following statements provide an example of this technique:

start norpdf2(t) global(yv,rho,omrho2,count); /*-----------------------------------------------------*/ /* This function is the density function and requires */ /* the variable T (passed in the argument) */ /* and a list of parameters, YV, RHO, OMRHO2, COUNT */ /* (defined in the GLOBAL clause) */ /*-----------------------------------------------------*/ count = count+1; q=(t#t-2#rho#t#yv+yv#yv)/omrho2; p=exp(-q/2); return(p); finish; start marginal(v) global(yy,yv,eps); /*--------------------------------------------------*/ /* The inner integral */ /* The limits of integration from .M to YY */ /* YV is passed as a parameter to the inner integral*/ /*--------------------------------------------------*/ interval = .M || yy; if ( v < -12 ) then return(0); yv = v; call quad(pm,"NORPDF2",interval) eps=eps; return(pm); finish; start norcdf2(a,b,rrho) global(yy,rho,omrho2,eps); /*------------------------------------------------*/ /* Post some global parameters */ /* YY, RHO, OMRHO2 */ /* EPS will be set from IML */ /* RHO and B cannot be arguments in the GLOBAL */ /* list at the same time */ /*------------------------------------------------*/ rho = rrho; yy = b; omrho2 = 1-rho#rho; /*----------------------------------------------*/ /* The outer integral */ /* The limits of integration */ /*----------------------------------------------*/ interval= .M || a;

/*----------------------------------------------*/ /*Note that EPS the keyword = EPS the variable */ /*----------------------------------------------*/ call quad(p,"MARGINAL",interval) eps=eps; /*--------------------------*/ /* PER will be reset here */ /*--------------------------*/ per = 1/(8#atan(1)#sqrt(omrho2)) * p; return(per); finish; /*----------------------------------*/ /*First set up the global constants */ /*----------------------------------*/ count = 0; eps = 1E-11; /*------------------------------------*/ /* Do the work and print the results */ /*------------------------------------*/ p = norcdf2(2,1,0.1); print p[format=E21.14]; print count;

The variable COUNT contains the number of times the NORPDF2 module is called. Note that the value computed by the NORCDF2 module is very close to that returned by the PROBBNRM function, as computed by the following statements:

/*------------------------------------------------*/ /* Calculate the value with the PROBBNRM function */ /*------------------------------------------------*/ pp = probbnrm(2,1,0.1); print pp[format=E21.14];

- The iterated inner integral cannot
have a left endpoint of .
For large values of
*v*, the inner integral does not contribute to the answer but still needs to be calculated to the required relative accuracy. Therefore, either cut off the function (when ), as in the MARGINAL module in the preceding code, or have the intervals start from a reasonable cutoff value. In addition, the QUAD call stops if the integrands appear to be identically 0 (probably caused by underflow) over the interval of integration. - This method of integration (iterated,
one-dimensional integrals) is extremely conservative
and requires unnecessary function evaluations.
In this example, the QUAD call for the inner integration
lacks information about the final value that the QUAD
call for the outer integration is trying to refine.
The lack of communication between the two
QUAD routines can cause useless computations
to be performed in the inner integration.
To illustrate this idea, let the relative error be 1E-11 and let the answer delivered by the outer integral be 0.8, as in this example. Any computation of the inner execution of the QUAD call that yields 0.8E-11 or less will not contribute to the final answer of the QUAD call for the outer integral. However, the inner integral lacks this information, and for a given value of the parameter

*yv*, it attempts to compute an answer with much more precision than is necessary. The lack of communication between the two QUAD subroutines prevents the introduction of better cut-offs. Although this method can be inefficient, the final calculations are accurate.

Chapter Contents |
Previous |
Next |
Top |

Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.