Date: Mon, 23 Mar 2009 14:29:26 -0500
Reply-To: OR Stats <stats112@GMAIL.COM>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: OR Stats <stats112@GMAIL.COM>
Subject: Re: Numerical Integration
In-Reply-To: <6eca73440903110630r23a4d8c7sec6a8237f850d0dc@mail.gmail.com>
Content-Type: text/plain; charset=ISO-8859-1
SAS tech support says that CALL QUAD does perform numerical integration but
it does not handle functions with imaginary numbers.... what are our
options in SAS or elsewhere then?
On Wed, Mar 11, 2009 at 8:30 AM, OR Stats <stats112@gmail.com> wrote:
> The pdf is derive from the characteristic function, but unfortunately does
> not have closed form. So to do MLE for the distn parameters (theta1,
> theta2), the integral would be evaluated at each observed data point and the
> maximized product (or sum of log integrals) would be calculated and
> (mle_theta1,mle_theta2) would be saved. The pdf is the following
>
> theta1*c1* Integral { e^(-iwx) * sin [theta2*sqrt(2iw)] / sin
> [-pi*sqrt(2iw)] dw }
>
> where c1 is some known constant; theta2<0; i is imaginary number. The
> integrand is integrated over w and has imaginary number; therefore it is
> unclear how to approach this in SAS as well.
>
> On Tue, Mar 10, 2009 at 6:11 PM, Lewis Jordan <lewjord@uga.edu> wrote:
>
>> OR. I may be able to put a small data set together and send you my PROC
>> MODEL code. If your interested let me know. I am interested in knowing what
>> pdf your using? I may be able to think of another solution/work around.
>>
>> Lewis
>>
>> ---- Original message ----
>> >Date: Tue, 10 Mar 2009 14:52:23 -0500
>> >From: OR Stats <stats112@gmail.com>
>> >Subject: Re: Numerical Integration
>> >To: Lewis Jordan <lewjord@uga.edu>
>> >
>> > This sounds promising. The likelihood function
>> > (objective function) that I am trying to maximize in
>> > my PROC MODEL routine is an integral b.c. the
>> > density function is an integral (no closed form
>> > available). Unfortunately I do not have PROC
>> > IML. And the alternative is solve the integral in
>> > C++ and then bring it back in MATLAB environment,
>> > b.c. I do not know how SAS would handle this
>> > issue. Thank you!
>> >
>> > On Tue, Mar 10, 2009 at 11:41 AM, Lewis Jordan
>> > <lewjord@uga.edu> wrote:
>> >
>> > Here is a numerical integration routine I adapted
>> > to be called in a data step or in a proc such as
>> > PROC MODEL. We wrote code using the SAS/MACRO
>> > function to evaluate the integrals presented in
>> > the second example using the trapezoid rule, since
>> > the terms outside of the integrand are constants
>> > and parameters to be estimated. We then
>> > implemented this MACRO within the SAS/ETS Proc
>> > Model procedure, which evaluates the integral and
>> > simultaneously estimates the parameters. An
>> > iterative procedure is used, with a criterion for
>> > convergence being agreement of the estimated
>> > integral between successive iterations to within
>> > 1X^10-6.
>> >
>> > I give two examples below and compare the routine
>> > to the "QUAD" method in IML. This can be
>> > generalized to any other function.
>> >
>> > options symbolgen mlogic mprint formdlim="*";
>> >
>> > *the 2-point closed Newton-Cotes formula called
>> > the 'trapezoid rule' ;
>> >
>> > *problem 1;
>> > *f(x) = int [ x^2 , h=LL to UL ] ;
>> >
>> > options symbolgen mlogic mprint;
>> >
>> > %macro func(x,beta,fvalue,evals) ;
>> > &fvalue = &x**2; ;
>> > &evals = &evals + 1 ;
>> > %mend func ;
>> > %macro numint(x,beta,ll,ul);
>> > eps = 1e-6 ;
>> > jmax = 15 ;
>> > reldiff = eps + 1 ;
>> > a = &ll ; b = &ul ;
>> > evals = 0 ;
>> > %func(a,&beta, fa, evals) ;
>> > %func(b, &beta,fb, evals) ;
>> > newint = .5 * (b - a) * (fa + fb) ;
>> > oldint = newint ;
>> > reldiff = 1 + eps ;
>> > it = 1 ;
>> > do j = 1 to jmax while (reldiff > eps) ;
>> > tnm = it ;
>> > del = (b - a) / tnm ;
>> > sumfun = 0 ;
>> > z = a + .5 * del ;
>> > do i = 1 to it ;
>> > %func(z,&beta,fx, evals) ;
>> > sumfun = sumfun + fx ;
>> > z = z + del ;
>> > end ;
>> > newint = .5 * (newint + (b - a) * sumfun /
>> > tnm) ;
>> > reldiff = abs((newint - oldint) / oldint) ;
>> > put j it oldint newint reldiff ;
>> > oldint = newint ;
>> > it = it * 2 ;
>> > end ;
>> > if (j >= jmax & reldiff > eps) then do ;
>> > put "No convergence after" jmax "
>> > iterations." ;
>> > stop ;
>> > end ;
>> > put "Number of function evaluations: " evals ;
>> > put "Estimated integral = "
>> > newint ;
>> > %mend numint;
>> >
>> > %macro job;
>> > data numint;
>> > beta=0;
>> > ll=0;
>> > ul=5;
>> > %numint(x,beta,ll,ul);
>> > proc print;
>> > title1 "Numerical Integration In a Data Step";
>> > title2 "f(x) = x^2";
>> > var newint;
>> > run;
>> >
>> > %mend;
>> > %job;
>> >
>> > title1 "Numerical Integration In IML";
>> > title2 "f(x) = x^2";
>> > proc iml;
>> > start fun(x);
>> > v = x##2;
>> > return(v);
>> > finish;
>> >
>> > /* Call QUAD */
>> > a = { 0 5 };
>> > call quad(z,"fun",a);
>> > print z;
>> > run;quit;
>> >
>> > *problem 2;
>> >
>> > *y = int [ { (H-h)**(2*beta4)*exp(phi1*h/TH) ,
>> > h=LL to UL ] ;
>> >
>> > *The term resulting from the integration of the
>> > above equation is an incomplete gamma function.
>> > It should be noted that even though a closed
>> > form solution exists for the integral presented
>> > above,
>> > evaluation of it will result in a complex
>> > argument, i.e., an imaginary number.
>> > This arises from the negative coefficient in the
>> > second argument of the incomplete gamma function
>> > and
>> > takes one into the realm of complex analysis.
>> > Several methods are available for evaluating
>> > this
>> > complex function including a hypergeometric
>> > power series, Laurent series, or power series of
>> > the half argument,
>> > and numerical integration (Mathar 2004). We
>> > chose the 2-point closed Newton-Cotes formula
>> > called the
>> > 'trapezoid rule' for numerical integration. ;
>> >
>> > %macro func(TH,h,beta4,phi1, fvalue, evals) ;
>> > &fvalue =
>> > (&Th-&h)**(2*&beta4)*exp(&phi1*&h/&TH) ;
>> > &evals = &evals + 1 ;
>> > %mend func ;
>> > %macro numint(TH,beta4,phi1,ll,ul);
>> > eps = 1e-6 ;
>> > jmax = 15 ;
>> > reldiff = eps + 1 ;
>> > a = &ll ; b = &ul ;
>> > evals = 0 ;
>> > %func(&th,a,&beta4,&phi1, fa, evals) ;
>> > %func(&th,b, &beta4,&phi1,fb, evals) ;
>> > newint = .5 * (b - a) * (fa + fb) ;
>> > oldint = newint ;
>> > reldiff = 1 + eps ;
>> > it = 1 ;
>> > do j = 1 to jmax while (reldiff > eps) ;
>> > tnm = it ;
>> > del = (b - a) / tnm ;
>> > sumfun = 0 ;
>> > x = a + .5 * del ;
>> > do i = 1 to it ;
>> > %func(&th,x, &beta4,&phi1,fx, evals) ;
>> > sumfun = sumfun + fx ;
>> > x = x + del ;
>> > end ;
>> > newint = .5 * (newint + (b - a) * sumfun /
>> > tnm) ;
>> > reldiff = abs((newint - oldint) / oldint) ;
>> > put j it oldint newint reldiff ;
>> > oldint = newint ;
>> > it = it * 2 ;
>> > end ;
>> > if (j >= jmax & reldiff > eps) then do ;
>> > put "No convergence after" jmax "
>> > iterations." ;
>> > stop ;
>> > end ;
>> > put "Number of function evaluations: " evals ;
>> > put "Estimated integral = "
>> > newint ;
>> > %mend numint;
>> >
>> > %macro job;
>> > data numint;
>> > TH=20;
>> > h=20;
>> > phi1=-1.113;
>> > beta4=0.649943;
>> > ll=0;
>> > ul=h;
>> > %numint(TH,beta4,phi1,ll,ul);
>> > proc print;
>> > title1 "Numerical Integration In a Data Step";
>> > title2 "f(h) = (Th-h)**(2*beta4)*exp(phi1*h/TH) ";
>> > var newint;
>> > run;
>> > %mend;
>> > %job;
>> >
>> > %let phi1=-1.113;
>> > %let beta4=0.649943;
>> > %let TH=20;
>> > %let ll=0;
>> > %let UL=20;
>> > title1 "Numerical Integration In IML";
>> > title2 "f(h) = (Th-h)**(2*beta4)*exp(phi1*h/TH) ";
>> > proc iml;
>> > start fun(h);
>> > v = (&Th-h)##(2*&beta4)*exp(&phi1*h/&TH) ;
>> > return(v);
>> > finish;
>> >
>> > /* Call QUAD */
>> > a = { &LL &UL };
>> > call quad(z,"fun",a);
>> > print z;
>> > run;quit;
>> >
>> > ---- Original message ----
>> > >Date: Tue, 10 Mar 2009 05:59:11 -0700
>> > >From: Shiling Zhang <shiling99@YAHOO.COM>
>> > >Subject: Re: Numerical Integration
>> > >To: SAS-L@LISTSERV.UGA.EDU
>> > >
>> > >On Mar 10, 6:22 am, stats...@GMAIL.COM (OR Stats)
>> > wrote:
>> > >> Hello:
>> > >>
>> > >> Has anyone done numerical integration in SAS,
>> > other than routines like FFT?
>> > >> What about Gauss-Laguerre etc.? Look forward
>> > to your input, or your trials
>> > >> & tribulations.
>> > >>
>> > >> Many thanx,
>> > >> OR Stats
>> > >
>> > >You can define your own numerical integration
>> > with SAS/IML call quad
>> > >routine. Here is the web link.
>> > >
>> > >
>> http://support.sas.com/onlinedoc/913/getDoc/en/imlug.hlp/langref_sect206.htm
>>
>
>
| 
