Commit ffd38c32 authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-3

parent 7ef10a6a
** denotes quite substantial/important changes
*** denotes really big changes
1.8-3
* Fix of two illegal read/write bugs with extended family models with no
smooths. (Thanks to Julian Faraway for reporting beta regr problem).
* bam now checks that chunk.size > number of parameters and resets the
chunk.size if not.
* Examples of use of smoothCon and PredictMat for setting up bases
for use outside mgcv (and then predicting) added to ?smoothCon.
1.8-2
* For exponential family gams, fitted by outer iteration, a warning is now
generated if the Pearson scale parameter estimate is more than twice
generated if the Pearson scale parameter estimate is more than 4 times
a robust estimate. This may indicate an unstable Pearson estimate.
* 'gam.control' now has an option 'scale.est' to allow selection of the
......
Package: mgcv
Version: 1.8-2
Version: 1.8-3
Author: Simon Wood <simon.wood@r-project.org>
Maintainer: Simon Wood <simon.wood@r-project.org>
Title: Mixed GAM Computation Vehicle with GCV/AIC/REML smoothness
......@@ -14,7 +14,7 @@ Suggests: splines, parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2014-07-21 08:56:24 UTC; sw283
Packaged: 2014-08-29 10:10:03 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2014-07-21 12:12:11
Date/Publication: 2014-08-29 22:07:35
10b8bc7c708c2f3d89cc69d8bd67f3e3 *ChangeLog
ad04030a1a5f9a3888602dee013351e7 *DESCRIPTION
c0324227202a615e286343615936dd1b *ChangeLog
149d2f3b36649e882a15d0092cf01f55 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
870396939733c58c9d523cbc41f9da59 *NAMESPACE
ca641f4a3db11cd9494762ad2286665c *R/bam.r
b5f01a8f420b0d296352d500b75e4c35 *R/bam.r
9d5bb48cc0072b9ba638567986928a73 *R/coxph.r
d6f40f3d2acd3e30f3a59cbc85da6bc3 *R/efam.r
2fb8891700c1a813dc387cdce8dd817b *R/fast-REML.r
......@@ -15,7 +15,7 @@ ee07c79f181f6cc17edeafc8a33ae51b *R/mgcv.r
ed2cdcfa5c0f772bb75eee723b7294cd *R/misc.r
c287514d201a02f060fee98aab1e93c8 *R/mvam.r
4b0df0a2310fba6b38ea8164f76d0ccd *R/plots.r
8d3e8c0e0009b8028bf08a2b1919f850 *R/smooth.r
3f10d9fbae29ae28bf89499b874f9aa4 *R/smooth.r
90f72ef91dfc2a25c262fc002f781bf7 *R/soap.r
0c5d2acdc985d2af9de6bb736d4df2ab *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
......@@ -28,7 +28,7 @@ e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
f12468625253dd3603907de233762fd6 *man/Rrank.Rd
f6cadda5999c35800fd65c08d6812f7b *man/Tweedie.Rd
80f8763baa4987579e2aa56073a9e94e *man/anova.gam.Rd
49fe84873e296267a96ed44c044d71c4 *man/bam.Rd
9c392eb571b80b00fc6fada9d502b69f *man/bam.Rd
71bde8b8caa24a36529ce7e0ac3165d8 *man/bam.update.Rd
a2beb811b1093c5e82ef32d7de1f7d32 *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
......@@ -118,7 +118,7 @@ a8d5eb4772641567fee15714f01a4727 *man/smooth.construct.so.smooth.spec.Rd
8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
3e6d88ef6a8ab21bd6f120120602dcf6 *man/smooth.construct.tp.smooth.spec.Rd
75a042d74460dd0c6e0c0a95c1f5d3f7 *man/smooth.terms.Rd
0d12daea17e0b7aef8ab89b5f801adf1 *man/smoothCon.Rd
3064042473b909ba6bbd160c476059d3 *man/smoothCon.Rd
b55a396da77559dac553613146633f97 *man/sp.vcov.Rd
83bd8e097711bf5bd0fff09822743d43 *man/spasm.construct.Rd
700699103b50f40d17d3824e35522c85 *man/step.gam.Rd
......@@ -145,11 +145,11 @@ bd5ec325242ca6e41056c2f1f7f05bfe *po/pl.po
749ac663240fd6eaa2b725602d47ef2a *po/po.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
56d501c682dc6699ff0803f25533e849 *src/coxph.c
c8dbbf1178ace1131a5d6d7ccfc2b905 *src/gdi.c
114701c0d407f946ac693da55da864da *src/gdi.c
49af97195accb65adc75620183d39a4c *src/general.h
94a7b06d798739744ef9342a6a7051ec *src/init.c
6dbd109048a3bd18b67db6230c064c21 *src/magic.c
5b3879c26540dc9eb1d4a8a9a575e975 *src/mat.c
a6527e26dde591b838d4c3c3f5a4e036 *src/mat.c
02122df2193f86e2a08f68fe5b8ff972 *src/matrix.c
6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
7a12ef0cbd8f5312cd67669a36d75ebf *src/mgcv.c
......
......@@ -322,7 +322,7 @@ bgam.fit <- function (G, mf, chunk.size, gp ,scale ,gamma,method, coef=NULL,etas
}
} else { ## use new parallel accumulation
for (i in 1:length(arg)) arg[[i]]$coef <- coef
res <- parLapply(cl,arg,qr.up)
res <- parallel::parLapply(cl,arg,qr.up)
## single thread debugging version
#res <- list()
#for (i in 1:length(arg)) {
......@@ -767,7 +767,7 @@ predict.bam <- function(object,newdata,type="link",se.fit=FALSE,terms=NULL,
if (!missing(newdata)) rm(newdata)
rm(object)
gc()
res <- parLapply(cluster,arg,pabapr) ## perform parallel prediction
res <- parallel::parLapply(cluster,arg,pabapr) ## perform parallel prediction
gc()
## and splice results back together...
if (type=="lpmatrix") {
......@@ -916,7 +916,7 @@ bam.fit <- function(G,mf,chunk.size,gp,scale,gamma,method,rho=0,
}
} else { ## use parallel accumulation
res <- parLapply(cl,arg,ar.qr.up)
res <- parallel::parLapply(cl,arg,ar.qr.up)
## Single thread de-bugging...
# res <- list()
# for (i in 1:length(arg)) {
......@@ -1207,7 +1207,12 @@ bam <- function(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,n
if (is.null(G$offset)) G$offset <- rep(0,n)
if (ncol(G$X)>nrow(mf)) stop("Model has more coefficients than data") ##+nrow(G$C)) stop("Model has more coefficients than data")
if (ncol(G$X) > chunk.size) { ## no sense having chunk.size < p
chunk.size <- 4*ncol(G$X)
warning(gettextf("chunk.size < number of coefficients. Reset to %d",chunk.size))
}
G$cl<-cl;
G$am <- am
......
......@@ -179,10 +179,10 @@ cSplineDes <- function (x, knots, ord = 4)
## copy end intervals to start, for wrapping purposes...
knots <- c(k1-(knots[nk]-knots[(nk-ord+1):(nk-1)]),knots)
ind <- x>xc ## index for x values where wrapping is needed
X1 <- splineDesign(knots,x,ord,outer.ok=TRUE)
X1 <- splines::splineDesign(knots,x,ord,outer.ok=TRUE)
x[ind] <- x[ind] - max(knots) + k1
if (sum(ind)) {
X2 <- splineDesign(knots,x[ind],ord,outer.ok=TRUE) ## wrapping part
X2 <- splines::splineDesign(knots,x[ind],ord,outer.ok=TRUE) ## wrapping part
X1[ind,] <- X1[ind,] + X2
}
X1 ## final model matrix
......@@ -1564,7 +1564,7 @@ smooth.construct.ps.smooth.spec<-function(object,data,knots)
if (length(k)!=nk+2*m[1]+2)
stop(paste("there should be ",nk+2*m[1]+2," supplied knots"))
}
object$X <- spline.des(k,x,m[1]+2,x*0)$design # get model matrix
object$X <- splines::spline.des(k,x,m[1]+2,x*0)$design # get model matrix
if (!is.null(k)) {
if (sum(colSums(object$X)==0)>0) warning("knot range is so wide that there is *no* information about some basis coefficients")
}
......@@ -1596,12 +1596,12 @@ Predict.matrix.pspline.smooth<-function(object,data)
n <- length(x)
ind <- x<=ul & x>=ll ## data in range
if (sum(ind)==n) { ## all in range
X <- spline.des(object$knots,x,m)$design
X <- splines::spline.des(object$knots,x,m)$design
} else { ## some extrapolation needed
## matrix mapping coefs to value and slope at end points...
D <- spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
D <- splines::spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
X <- matrix(0,n,ncol(D)) ## full predict matrix
X[ind,] <- spline.des(object$knots,x[ind],m)$design ## interior rows
X[ind,] <- splines::spline.des(object$knots,x[ind],m)$design ## interior rows
## Now add rows for linear extrapolation...
ind <- x < ll
if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
......
......@@ -17,8 +17,8 @@ for large datasets. \code{bam} can also compute on a cluster set up by the \link
bam(formula,family=gaussian(),data=list(),weights=NULL,subset=NULL,
na.action=na.omit, offset=NULL,method="fREML",control=list(),
scale=0,gamma=1,knots=NULL,sp=NULL,min.sp=NULL,paraPen=NULL,
chunk.size=10000,rho=0,AR.start=NULL,sparse=FALSE,cluster=NULL,gc.level=1,
use.chol=FALSE,samfrac=1,drop.unused.levels=TRUE,...)
chunk.size=10000,rho=0,AR.start=NULL,sparse=FALSE,cluster=NULL,
gc.level=1,use.chol=FALSE,samfrac=1,drop.unused.levels=TRUE,...)
}
%- maybe also `usage' for other objects documented here.
......@@ -103,7 +103,8 @@ if smooths share smoothing parameters).}
\code{\link{gam.models}} explains more.}
\item{chunk.size}{The model matrix is created in chunks of this size, rather than ever being formed whole.}
\item{chunk.size}{The model matrix is created in chunks of this size, rather than ever being formed whole.
Reset to \code{4*p} if \code{chunk.size < 4*p} where \code{p} is the number of coefficients.}
\item{rho}{An AR1 error model can be used for the residuals (based on dataframe order), of Gaussian-identity
link models. This is the AR1 correlation parameter.}
......
......@@ -102,6 +102,31 @@ The parameterization used by \code{\link{gam}} can be controlled via
\seealso{ \code{\link{gam.control}},
\code{\link{smooth.construct}}, \code{\link{Predict.matrix}} }
\examples{
## example of using smoothCon and PredictMat to set up a basis
## to use for regression and make predictions using the result
library(MASS) ## load for mcycle data.
## set up a smoother...
sm <- smoothCon(s(times,k=10),data=mcycle,knots=NULL)[[1]]
## use it to fit a regression spline model...
beta <- coef(lm(mcycle$accel~sm$X-1))
with(mcycle,plot(times,accel)) ## plot data
times <- seq(0,60,length=200) ## creat prediction times
## Get matrix mapping beta to spline prediction at 'times'
Xp <- PredictMat(sm,data.frame(times=times))
lines(times,Xp\%*\%beta) ## add smooth to plot
## Same again but using a penalized regression spline of
## rank 30....
sm <- smoothCon(s(times,k=30),data=mcycle,knots=NULL)[[1]]
E <- t(mroot(sm$S[[1]])) ## square root penalty
X <- rbind(sm$X,0.1*E) ## augmented model matrix
y <- c(mcycle$accel,rep(0,nrow(E))) ## augmented data
beta <- coef(lm(y~X-1)) ## fit penalized regression spline
Xp <- PredictMat(sm,data.frame(times=times)) ## prediction matrix
with(mcycle,plot(times,accel)) ## plot data
lines(times,Xp\%*\%beta) ## overlay smooth
}
\keyword{models} \keyword{smooth} \keyword{regression}%-- one or more ...
File mode changed from 100644 to 100755
......@@ -71,12 +71,14 @@ double diagABt(double *d,double *A,double *B,int *r,int *c)
r by c stored column-wise.
*/
{ int j;
double tr,*pa,*pb,*p1,*pd;
for (pa=A,pb=B,p1=pa + *r,pd=d;pa<p1;pa++,pb++,pd++) *pd = *pa * *pb;
for (j=1;j < *c;j++)
for (p1=d + *r,pd=d;pd<p1;pa++,pb++,pd++) *pd += *pa * *pb;
/* d now contains diag(AB') */
for (tr=0.0,pd=d,p1=d + *r;pd < p1;pd++) tr += *pd;
double tr=0.0,*pa,*pb,*p1,*pd;
if (*c>0) {
for (pa=A,pb=B,p1=pa + *r,pd=d;pa<p1;pa++,pb++,pd++) *pd = *pa * *pb;
for (j=1;j < *c;j++)
for (p1=d + *r,pd=d;pd<p1;pa++,pb++,pd++) *pd += *pa * *pb;
/* d now contains diag(AB') */
for (pd=d,p1=d + *r;pd < p1;pd++) tr += *pd;
}
return(tr);
}
......@@ -1108,7 +1110,9 @@ void get_ddetXWXpS(double *det1,double *det2,double *P,double *K,double *sp,
rSoff = (int *)R_chk_calloc((size_t)*M,sizeof(int));
rSoff[0] = 0;for (m=0;m < *M-1;m++) rSoff[m+1] = rSoff[m] + rSncol[m];
if (*M>0) {
rSoff[0] = 0;for (m=0;m < *M-1;m++) rSoff[m+1] = rSoff[m] + rSncol[m];
}
tid = 0;
#ifdef SUPPORT_OPENMP
#pragma omp parallel private(m,bt,ct,tid) num_threads(nthreads)
......
......@@ -577,7 +577,7 @@ void getXtWX(double *XtWX, double *X,double *w,int *r,int *c,double *work)
else for (j=0;j<=i;j++) XtWX[i * *c + j] = w2[j];
}
/* now fill in the other triangle */
*XtWX = xx;
if (*r * *c > 0) *XtWX = xx;
for (i=0;i< *c;i++) for (j=0;j<i;j++) XtWX[j * *c + i] = XtWX[i * *c + j];
}
......
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