From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Fri 05 May 2006 - 04:52:16 EST

*>
*

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

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.

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

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.
*