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

[Spencer Graves]

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 at 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