[R] : Spline integration & Gaussian quadrature (was: gauss.quad.prob)

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Fri 05 May 2006 - 04:52:16 EST

Hi, Paul, Martin, Doug, Bill, et al.:

          What facilities are available in R for spline integration in R? RSiteSearch("spline integration") led me to a suggestion from Martin to feed "splinefun" to "integrate".
(http://finzi.psych.upenn.edu/R/Rhelp02a/archive/44018.html). It also led me to "integrate.xy(..., use.spline=TRUE, ...)", which I think I should study carefully.

          We could, for example, make "integrate" generic functoin with a special method for certain types of spline objects -- and write a complementary method for "deriv".

          This could be quite powerful for many applications. This could be used with multi-level modeling, for example, with potentially expensive function evaluations being saved and combined with later evaluations to produce splines from which gradients and hessians of integrals would be almost trivial. Another class of applications could involve fitting multi-dimensional splines to functions of, say, one variable of integration plus other parameters; resulting spline could then be packaged in functions that would compute very quickly. This could provide cheap, accurate answers to integrals that are currently much more expensive and difficult to evaluation.

          I'm willing to work on this, but I could use some suggestions on where to start (and other assistance if anyone feels so inclined).

          I've skimmed Wahba(1990) Spline Models of Observational Data (SIAM), but I so far have not bridged the gap between that book and R code for splines. My literature review then led me to Paul Dierckx (1995) Curve and Surface Fitting with Splines (Oxford). This book has a companion package of Fortran routine called "dierckx" available from "www.netlib.org"; Dierckx has agreed we could distribute it under the GNU license as part of an R package (provided we acknowledge his contribution). The book also contains a chapter on "fitting with convexity constraints" plus a section on "infinite end point derivatives", which looked useful for some of my work.

          My tentative work plan at the moment is to start by reviewing "integrate.xy", then possibly other functions like D1ss and D2ss, also in "sfsmisc", plus Paul Gilbert's "numDeriv" and the other spline functions in R. My goal at this stage is to understand spline estimation algorithms and how coefficients are stored. At that point, I will decide what I should do next.

          I also wonder about the feasibility of creating "splines" that merge to forms like exp(b*x) as x -> (+/-Inf) or perhaps (x-x0)*exp(-b*(x-x0)^2. This would facilitate integration and differentiation with half- or doubly-infinite ranges.

          Comments? Suggestions?         

GAUSSIAN QUADRATURE?           I started this email as a reply to a question from Harold Doran on "gauss.quad.prob": Gaussian quadrature of order n approximates a definite integral from a to b of f(x)*w(x)*dx with a weighted sum i = 1:n of w[i]*f(x[i]), where the quadrature points and weights are chosen so the formula is exact if f(x) is a polynomial of order at most (2*n-1). If a=(-Inf), b=Inf, and w(x) = a normal density, this is called Gauss-Hermite quadrature; if the mean and standard deviation of the normal are estimated from the data, it is called "adaptive Gauss-Hermite quadrature". If a and b are finite and w(x) is a beta distribution on that range, this is called Gauss-Jacobi quadrature; for the uniform distribution, it is called Gauss-Legendre. For the gamma distribution (including exponential as a special case), the integral is from 0 to Inf, and it's called Gauss-Laguerre. Recent references that I found useful for this are the following:

          Kythe and Schaeferkotter (2005) Handbook of Computational Methods for Integration (Chapman and Hall)

          Evans and Schwartz (2000) Approximating Integrals via Monte Carlo and Deterministic Methods.

          In theory, one could use Gauss-Jacobi on any finite interval, Gauss-Laguerre on any semi-finite interval and Gauss-Hermite on any integral over the entire real line. However, in my recent attempts to do this, I failed to get reasonable rates of convergence for many integrals that I thought should be reasonably well behaved.

          However, as I mentioned above, my work in this area suggested I might more wisely invest my time in spline integration.

	  Best Wishes,
	  spencer graves

Doran, Harold wrote:

> I've written a series of functions that evaluates an integral
from -inf to a or b to +inf using equally spaced quadrature points along a normal distribution from -10 to +10 moving in increments of .01. These functions are working and give very good approximations, but I think they are computationally wasteful as I am evaluating the function at *many* points.
> Instead, I would prefer to use true Gaussan Quadrature and
evaluate the integral at fewer points. I am working with the gauss.quad.prob function and have rewritten my functions to use the nodes and weights. But, I believe these are only for definite integrals and not for evaluating from -inf to a or b to +inf.

>>From what I understand, I cannot use the tails (e.g., -4.85 or 4.85 below) to approximate inf. In fact, I've tried it and it doesn't work. Instead, I believe one needs special weights when the interest is to evaluate an integral from -inf to a etc.
> I've read through the help and searched the archives and didn't
see any relevant discussion on whether it is possible to use this function to get the weights for -inf to a. Only two hits come up in the archive and neither was relevant. I think the Laguerre quadrature in gauss.prob might be the right path, but I'm uncertain.
> Is it possible to get the nodes and weights I need from R
directly, or should I refer to a table that has derived these for the normal distribution from another source?
> thank you
> Harold

>>gauss.quad.prob(10, dist='normal') 

> $nodes
> [1] -4.8594628 -3.5818235 -2.4843258 -1.4659891 -0.4849357 0.4849357 1.4659891 2.4843258 3.5818235
> [10] 4.8594628
> $weights
> [1] 4.310653e-06 7.580709e-04 1.911158e-02 1.354837e-01 3.446423e-01 3.446423e-01 1.354837e-01 1.911158e-02
> [9] 7.580709e-04 4.310653e-06
> [[alternative HTML version deleted]]
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri May 05 16:54:23 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Fri 05 May 2006 - 18:10:01 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.