[Lme4-commits] r1790 - in pkg/Gqr: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 20 18:27:50 CEST 2012
Author: mmaechler
Date: 2012-10-20 18:27:49 +0200 (Sat, 20 Oct 2012)
New Revision: 1790
Added:
pkg/Gqr/man/integrateGQ.Rd
Modified:
pkg/Gqr/DESCRIPTION
pkg/Gqr/NAMESPACE
pkg/Gqr/R/GaussQuad.R
pkg/Gqr/man/GaussQuad.Rd
pkg/Gqr/tests/GH_test.R
Log:
add integrateGQ() and some interesting examples
Modified: pkg/Gqr/DESCRIPTION
===================================================================
--- pkg/Gqr/DESCRIPTION 2012-10-03 14:03:37 UTC (rev 1789)
+++ pkg/Gqr/DESCRIPTION 2012-10-20 16:27:49 UTC (rev 1790)
@@ -1,6 +1,6 @@
Package: Gqr
-Version: 0.1.0
-Date: 2011-11-30
+Version: 0.2-0
+Date: 2012-10-20
Title: Determine knots and weights for Gaussian quadrature rules
Author: Douglas Bates <bates at stat.wisc.edu>,
Martin Maechler <maechler at R-project.org> and
@@ -12,7 +12,8 @@
kernels can be used.
Depends: R(>= 2.14.0)
LinkingTo: Rcpp
+Suggests: sfsmisc
+SuggestsNote: only for examples, etc
LazyLoad: yes
License: GPL (>=2)
-URL: http://lme4.r-forge.r-project.org/
-Packaged: 2011-11-29 20:07:15 UTC; bates
+
Modified: pkg/Gqr/NAMESPACE
===================================================================
--- pkg/Gqr/NAMESPACE 2012-10-03 14:03:37 UTC (rev 1789)
+++ pkg/Gqr/NAMESPACE 2012-10-20 16:27:49 UTC (rev 1790)
@@ -1,2 +1,2 @@
useDynLib(Gqr, .registration=TRUE)
-export(GaussQuad)
+export(GaussQuad, integrateGQ)
Modified: pkg/Gqr/R/GaussQuad.R
===================================================================
--- pkg/Gqr/R/GaussQuad.R 2012-10-03 14:03:37 UTC (rev 1789)
+++ pkg/Gqr/R/GaussQuad.R 2012-10-20 16:27:49 UTC (rev 1790)
@@ -1,15 +1,32 @@
GaussQuad <-
- function(n, rule, a=0, b=1, alpha=0, beta=0) {
- rules <- c("Legendre", "Chebyshev", "Gegenbauer", "Jacobi",
- "Laguerre", "Hermite", "Exponential", "Rational", "Type2")
- if (is.na(rr <- pmatch(rule, rules)))
- stop("Unknown quadrature rule, ", rule)
- stopifnot(length(n <- as.integer(n)) == 1L,
- length(a <- as.numeric(a)) == 1L,
- length(b <- as.numeric(b)) == 1L,
- length(alpha <- as.numeric(alpha)) == 1L,
- length(beta <- as.numeric(beta)) == 1L,
- n > 0L)
- .Call(Gauss_Quad, n, rr, a, b, alpha, beta)
- }
+ function(n, rule = c("Legendre", "Chebyshev", "Gegenbauer", "Jacobi",
+ "Laguerre", "Hermite", "Exponential", "Rational"),
+ a=0, b=1, alpha=0, beta=0)
+{
+ ## Note: rule "Type2" is not yet supported in C code
+ rules <- eval(formals()$rule)
+ if (is.na(rr <- pmatch(rule, rules)))
+ stop("Unknown quadrature rule, ", rule)
+ stopifnot(length(n <- as.integer(n)) == 1L,
+ length(a <- as.numeric(a)) == 1L,
+ length(b <- as.numeric(b)) == 1L,
+ length(alpha <- as.numeric(alpha)) == 1L,
+ length(beta <- as.numeric(beta)) == 1L,
+ n > 0L)
+ .Call(Gauss_Quad, n, rr, a, b, alpha, beta)
+}
+integrateGQ <- function(f, lower, upper, ...,
+ n = 13, rule = "Legendre",
+ a=lower, b=upper, alpha=0, beta=0)
+{
+ f <- match.fun(f)
+ stopifnot(length(a) == 1, is.finite(a),
+ length(b) == 1, is.finite(b),
+ length(alpha) == 1, is.finite(alpha),
+ length(beta) == 1, is.finite(beta))
+ rules <- eval(formals(GaussQuad)$rule)
+ rule <- match.arg(rule, choices = rules)
+ GQ <- GaussQuad(n, rule=rule, a=a, b=b, alpha=alpha, beta=beta)
+ with(GQ, sum(weights * f(knots, ...)))
+}
Modified: pkg/Gqr/man/GaussQuad.Rd
===================================================================
--- pkg/Gqr/man/GaussQuad.Rd 2012-10-03 14:03:37 UTC (rev 1789)
+++ pkg/Gqr/man/GaussQuad.Rd 2012-10-20 16:27:49 UTC (rev 1790)
@@ -11,12 +11,14 @@
generalize the rules.
}
\usage{
-GaussQuad(n, rule, a = 0, b = 1, alpha = 0, beta = 0)
+GaussQuad(n, rule = c("Legendre", "Chebyshev", "Gegenbauer", "Jacobi",
+ "Laguerre", "Hermite", "Exponential", "Rational"),
+ a = 0, b = 1, alpha = 0, beta = 0)
}
\arguments{
- \item{n}{integer - number of quadrature points}
- \item{rule}{character - a partial match to one of the quadrature rule
- names. See the details section.}
+ \item{n}{integer - number of quadrature points.}
+ \item{rule}{character - a partial match to one of the quadrature
+ rules, defaulting to \code{"Legendre"}, see the details section.}
\item{a}{numeric scalar - an optional shift parameter or interval endpoint}
\item{b}{numeric scalar - an optional scale parameter or interval endpoint}
\item{alpha}{numeric scalar - an optional power}
@@ -34,6 +36,7 @@
\item{\dQuote{Hermite}}{\eqn{\int_{-\infty}^\infty f(x)|x-a|^\alpha \exp(-b*(x-a)^2)\,dx}{int_-inf^inf f(x)|x-a|^alpha*exp(-b*(x-a)^2) dx}}
\item{\dQuote{Exponential}}{\eqn{\int_a^b f(x)|x-(a+b)/2.0|^\alpha\,dx}{int_-inf^inf f(x)|x-(a+b)/2.0|^alpha dx}}
\item{\dQuote{Rational}}{\eqn{\int_a^\infty f(x)(x-a)^\alpha*(x+b)^\beta\,dx}{int_a^inf f(x)(x-a)^alpha*(x+b)^beta}}
+ \item{\dQuote{Type2}}{\eqn{\int_a^\infty f(x)(x-a)^\alpha*(x+b)^\beta\,dx}{int_a^inf f(x)(x-a)^alpha*(x+b)^beta}}
}
}
\value{
@@ -44,7 +47,7 @@
\references{
The original FORTRAN implementation is from
Sylvan Elhay, Jaroslav Kautsky (1987),
- \dQuote{Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
+ \dQuote{Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
Interpolatory Quadrature}, \bold{ACM Transactions on Mathematical Software},
Volume \bold{13}, Number 4, December 1987, pages 399--415.
Added: pkg/Gqr/man/integrateGQ.Rd
===================================================================
--- pkg/Gqr/man/integrateGQ.Rd (rev 0)
+++ pkg/Gqr/man/integrateGQ.Rd 2012-10-20 16:27:49 UTC (rev 1790)
@@ -0,0 +1,104 @@
+\name{integrateGQ}
+\alias{integrateGQ}
+\title{Numerical Integration via Gauss-Hermite}
+\description{
+ Numerical integration of one-dimensional functions, using
+ Gauss-Hermite, with respect to eight different kernel (\code{rule}s).
+}
+\usage{
+integrateGQ(f, lower, upper, ..., n = 13, rule = "Legendre",
+ a = lower, b = upper, alpha = 0, beta = 0)
+}
+\arguments{
+ \item{f}{an \R \code{\link{function}} taking a numeric first argument,
+ the \sQuote{integrand} of the numerical integration (aka
+ \dQuote{quadrature}).}
+ \item{lower, upper}{the limits of integration. Currently \emph{must}
+ be finite.}
+ \item{\dots}{additional arguments to be passed to \code{f}.}
+
+ \item{n}{number of Gauss-Hermite quadrature points, see
+ \code{\link{GaussQuad}}.}
+ \item{rule}{a character string specifying one of the quadrature rules,
+ see \code{\link{GaussQuad}}.}
+ \item{a,b, alpha, beta}{numeric scalars, specifying the integral to be
+ computed; the meaning of these depends much on the \code{rule}, see
+ the \code{Details} in \code{\link{GaussQuad}}.}
+}
+\value{
+ a number, simply
+ \deqn{sum_{i=1}^n w_i f(x_i),}{sum(i=1:n; w[i] * f(x[i])),}
+ where \eqn{x_i}{x[i]} are the \code{knots}
+ and \eqn{w_i}{w[i]} are the \code{weights} returned from
+ \code{\link{GaussQuad}(n, rule, ..)}.
+}
+%%\author{Martin Maechler}
+\seealso{
+ \code{\link{integrate}()}, \R's \dQuote{built-in} numerical integration.
+ \code{\link{GaussQuad}}, also for references.
+}
+\examples{
+ f1 <- function(x) x^2
+ (i1 <- integrate (f1, -1, 3))
+ (i2 <- integrateGQ(f1, -1, 3))
+ stopifnot(all.equal(i2, 28/3, tol = 1e-12), all.equal(i1$value, i2))
+ ## n does not matter for such a simple f(): always ~ 14 digits accuracy:
+ for(n in 2:20) cat(sprintf("\%3d: \%5.2f\n", n,
+ log10(abs(28/3 - integrateGQ(f1, -1, 3, n=n)))))
+
+ ## However, it *does* for Hermite, and a "hard" case:
+ gi3 <- function(x) 1/(1+ abs(x)^3)
+ n. <- c(3:130,150, 170, 200, 250, 350, 500, 1000)
+ Eg. <- sapply(n., function(n)
+ integrateGQ(gi3, n=n, rule = "Hermite", a=0, b=1/2))
+ ## (it *is* very fast)
+ plot(Eg. ~ n., type = "l", col=2,
+ subset = n. <= 130, ylim = c(1.65, 1.73))
+
+ ## "Classical" numerical integration:
+ (Icl <- integrate(function(u) gi3(u)*exp(-u^2/2), -Inf,Inf,
+ rel.tol= 1e-13))
+ ## 1.6810227521877 with absolute error < 5.2e-14
+ Ic <- Icl$value # (let's assume it's "true")
+
+ xax <- require("sfsmisc")
+ plot(abs(Eg. - Ic) ~ n., type = "o", col=2, log="xy",
+ axes = !xax,
+ main = "Gauss-Hermite approximation error [log-log scale]")
+ if(xax) { eaxis(1); eaxis(2) }
+
+ ## Compute E[ g(X) ] where X ~ N(0, 1) via Gauss-Hermite
+ E.norm <- function(f, n = 32)
+ integrateGQ(f, n=n, rule = "Hermite", a=0, b=1/2) / sqrt(2*pi)
+
+ integrands <- list(
+ g1 = function(x) (x-2)^2,
+ g2 = function(x) (x+3)^3,
+ g4 = function(x) x^4,
+ gg = function(x) 1/(1+ abs(x)^3) )
+
+ op <- par(mfrow= c(2,2), mgp=c(1.5,.6, 0), mar=.1+c(3,3,3,0))
+ for(i in seq_along(integrands)) {
+ g <- integrands[[i]]
+ nm <- names(integrands)[i]
+ cat(nm,":", deparse(expr <- body(g)), ":\n")
+ curve(g, -3, 3, main = expr); abline(h=0, lty=3)
+ cat("classical integrate():\n")
+ print(i1 <- integrate(function(u) g(u)*dnorm(u), -Inf, Inf))
+ cat("\n Gauss-Hermite (n= 5):", i2 <- E.norm(g, n= 5),"\n")
+ cat( " Gauss-Hermite (n=16):", i3 <- E.norm(g, n=16),"\n")
+ cat( " Gauss-Hermite (n=64):", i4 <- E.norm(g, n=64),"\n")
+ cat( " Gauss-Hermite (n=128):",i5 <- E.norm(g,n=128),"\n\n")
+ legend("top",
+ c(paste("E[g(X)] {integr.()} = ", formatC(i1$value)),
+ paste("E[g(X)] {Hermite_5 } =", formatC(i2)),
+ paste("E[g(X)] {Hermite_16} =", formatC(i3)),
+ paste("E[g(X)] {Hermite_64} =", formatC(i4)),
+ paste("E[g(X)] {Hermite_128}=", formatC(i5))),
+ bty="n")
+ }
+ par(op)
+
+}
+\keyword{math}
+\keyword{utilities}
Modified: pkg/Gqr/tests/GH_test.R
===================================================================
--- pkg/Gqr/tests/GH_test.R 2012-10-03 14:03:37 UTC (rev 1789)
+++ pkg/Gqr/tests/GH_test.R 2012-10-20 16:27:49 UTC (rev 1790)
@@ -1,10 +1,24 @@
library(Gqr)
-wfmt <- "http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_hermite_probabilist/hermite_probabilist_%03d_w.txt "
-xfmt <- "http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_hermite_probabilist/hermite_probabilist_%03d_x.txt "
+
+## Do not do the checks per default "on CRAN"
+(doExtras <- {
+ interactive() ||
+ grepl("bates", Sys.getenv("USER")) ||
+ nzchar(Sys.getenv("R_Gqr_check_extra")) ||
+ identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))
+})
+
+if(doExtras) {
+ urldir <- "http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_hermite_probabilist/"
+ wfmt <- paste0(urldir,"hermite_probabilist_%03d_w.txt")
+ xfmt <- paste0(urldir,"hermite_probabilist_%03d_x.txt")
+}
nn <- 2^(1:7) - 1L
for (n in nn) {
loc <- GaussQuad(n, "H", b=0.5)
- webw <- scan(url(sprintf(wfmt, n)))
- webx <- scan(url(sprintf(xfmt, n)))
- stopifnot(all.equal(webw, loc$weights), all.equal(webx, loc$knots))
+ if(doExtras) {
+ webw <- scan(url(sprintf(wfmt, n)))
+ webx <- scan(url(sprintf(xfmt, n)))
+ stopifnot(all.equal(webw, loc$weights), all.equal(webx, loc$knots))
+ }
}
More information about the Lme4-commits
mailing list