Commit 207dfa0a authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.7-4

parent 1f187945
Package: mgcv
Version: 1.7-3
Version: 1.7-4
Author: Simon Wood <>
Maintainer: Simon Wood <>
Title: GAMs with GCV/AIC/REML smoothness estimation and GAMMs by PQL
......@@ -12,6 +12,6 @@ Imports: graphics, stats, nlme, Matrix
Suggests: nlme (>= 3.1-64), splines, Matrix
LazyLoad: yes
License: GPL (>= 2)
Packaged: 2011-02-14 18:14:10 UTC; simon
Packaged: 2011-03-09 11:26:29 UTC; simon
Repository: CRAN
Date/Publication: 2011-02-14 19:21:19
Date/Publication: 2011-03-09 15:04:01
......@@ -2084,6 +2084,9 @@ totalPenaltySpace <- function(S,H,off,p)
## range space of the penalty, and obtain actual null space dimension
## components are roughly rescaled to avoid any dominating
Hscale <- sqrt(sum(H*H));
if (Hscale==0) H <- NULL ## H was all zeroes anyway!
if (is.null(H)) St <- matrix(0,p,p)
else { St <- H/sqrt(sum(H*H));
if (ncol(H)!=p||nrow(H)!=p) stop("H has wrong dimension")
......@@ -2,6 +2,17 @@
*** denotes really big changes
* Fix for single letter typo bug in C code called by slanczos, could
actually segfault on matrices of less than 10 by 10.
* matrix.c:Rlanczos memory error fix in convergence testing of -ve
* Catch for min.sp vector all zeroes, which could cause an ungraceful
** "ds" (Duchon splines) smooth class added. See ?Duchon.spline
......@@ -13,7 +13,7 @@ Duchon's (1977) construction generalizes the usual thin plate spline penalty as
of the squared Euclidian norm of a vector of mixed partial derivatives of the function w.r.t. its arguments. Duchon re-expresses
this penalty in the Fourier domain, and then weights the squared norm in the integral by the Euclidean norm of the fourier frequencies,
raised to the power 2s. s is a user selected constant taking integer values divided by 2. If d is the numberof arguments of the smooth,
then it is required that -d/2 < m < d/2. To obtain smooth functions we further require that m + s > n/2. If s=0 then the usual thin plate
then it is required that -d/2 < m < d/2. To obtain smooth functions we further require that m + s > d/2. If s=0 then the usual thin plate
spline is recovered.
The construction is amenable to exactly the low rank approximation method given in Wood (2003) to thin plate splines, with similar
......@@ -12,7 +12,7 @@ The first is taken to be latitude (in degrees) and the second longitude (in degr
For \code{m>0}, \code{(m+2)/2} is the penalty order, so \code{m=2} is the equivalent to the usual second
derivative penalty. \code{m=0} signals to use the 2nd order spline on the sphere, computed by
Wendelberger's (1981) method. \code{m = -1} results in a \code{\link{Duchon.spline}} being
used (with m=2 and s=1/2).
used (with m=2 and s=1/2), following an unpublished suggestion of Jean Duchon.
\code{k} (default 50) is the basis dimension.
......@@ -40,14 +40,18 @@ used (with m=2 and s=1/2).
smooth class documented under \code{\link{smooth.construct}}, this object will contain:
\item{Xu}{A matrix of the unique covariate combinations for this smooth (the basis is constructed by first stripping
out duplicate locations).}
\item{UZ}{The matrix mapping the t.p.r.s. parameters back to the parameters of a full thin plate spline.}
\item{UZ}{The matrix mapping the parameters of the reduced rank spline back to the parameters of a full spline.}
\details{ For \code{m>0}, the smooths implemented here are based on the pseudosplines on the sphere of Wahba (1981)
(there is a correction of table 1 in 1982, but the correction has a misprint in the definition of A --- the A
given in the 1981 paper is correct). For \code{m=0} then a second order spline on the sphere is used which is the analogue of
a second order thin plate spline in 2D: the computation is based on Chapter 4 of Wendelberger, 1981.
Optimal low rank approximations are obtained using exactly the approach given in Wood (2003).
Optimal low rank approximations are obtained using exactly the approach given in Wood (2003). For \code{m = -1} a smooth of the
general type discussed in Duchon (1977) is used: the sphere is embedded in a 3D Euclidean space, but
smoothing employs a penalty based on second derivatives (so that locally as the smoothing parameter tends
to zero we recover a "normal" thin plate spline on the tangent space). This is an unpublished suggestion of Jean
Note that the null space of the penalty is always the space of constant functions on the sphere, whatever the order of penalty.
#PKG_CFLAGS = -Wall -pedantic -g
# PKG_CFLAGS = -Wall -pedantic -g
## `#' out previous line for release
......@@ -2952,7 +2952,7 @@ void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm) {
if (*lm<0) { biggest=1;*lm=0;} /* get m largest magnitude eigen-values */
f_check = (*m + *lm)/2; /* how often to get eigen_decomp */
if (f_check<1) f_check ++;
kk = (int) floor(*n/10); if (kk<1) k=1;
kk = (int) floor(*n/10); if (kk<1) kk=1;
if (kk<f_check) f_check = kk;
q=(double **)calloc((size_t)(*n+1),sizeof(double *));
......@@ -3056,7 +3056,8 @@ void Rlanczos(double *A,double *U,double *D,int *n, int *m, int *lm) {
high_conv=0;i=0; while (i<=j&&err[i]<max_err&&d[i]>=0.0) { i++;high_conv++;}
for (i=0;i<low_conv;i++) /* test each of the negatives to see if it should be in the set of largest */
{ if (-d[j-i]>=d[high_conv-1]) { conv++;use_low++;}
{ if (high_conv==0) {conv++;use_low++;} else
if (-d[j-i]>=d[high_conv-1]) { conv++;use_low++;}
if (conv >= *m) /* then have enough - need to reset lm and m appropriately and break */
{ while (conv > *m) /* drop smallest magnitude terms */
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment