Commit 8b8f464f authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-1

parent 86e4d7fd
** denotes quite substantial/important changes
*** denotes really big changes
1.8-1
* bam modified so that choleski based fitting works properly with rank
deficient model matrix (without regularization).
* fix of 1.8-0 bug - gam prior weights mishandled in computation of cov matrix.
Thanks Fabian Scheipl.
1.8-0
*** Cox Proportional Hazard family 'cox.ph' added as example of general
......
Package: mgcv
Version: 1.8-0
Version: 1.8-1
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-06-25 05:49:03 UTC; sw283
Packaged: 2014-07-05 09:43:28 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2014-06-25 09:05:33
Date/Publication: 2014-07-06 00:16:40
064eafe39b266bd8f4c0f070467de907 *ChangeLog
2fb3b015c2e5425ebc81cb02b1252129 *DESCRIPTION
bda4042be09d0aa5d24bd6426ff711e3 *ChangeLog
eda47726d381d2a2e3bf0a34d8046426 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
d1659aca631e70420c5743e07bc71935 *NAMESPACE
8adb8ddf38c97741cb7713cc21117997 *R/bam.r
752f0cb574845c1596bf656e9ff48ad8 *R/bam.r
9d5bb48cc0072b9ba638567986928a73 *R/coxph.r
d6f40f3d2acd3e30f3a59cbc85da6bc3 *R/efam.r
2fb8891700c1a813dc387cdce8dd817b *R/fast-REML.r
db354e1dd127aca2e9eefc5494103050 *R/gam.fit3.r
7953f7d1ff8d5ba96778c84934620e64 *R/gam.fit3.r
6586de2772276237525c646362926f8b *R/gam.fit4.r
e63276b78a4ff2566af2e651bfc6aa9c *R/gam.sim.r
bd651b8eec80d83d6e0de5f47715a0de *R/gamlss.r
......
......@@ -29,29 +29,18 @@ rwMatrix <- function(stop,row,weight,X) {
chol2qr <- function(XX,Xy) {
## takes X'X and X'y and returns R and f
## equivalent to qr update. Uses simple
## regularization if XX not +ve def.
XXeps <- norm(XX)*.Machine$double.eps^.9
n <- ncol(XX)
for (i in 0:20) {
ok <- TRUE
R <- try(chol(XX+diag(n)*XXeps*i,pivot=TRUE),silent=TRUE) ## R'R = X'X
if (inherits(R,"try-error")) ok <- FALSE else {
ipiv <- piv <- attr(R,"pivot")
f <- try(forwardsolve(t(R),Xy[piv]),silent=TRUE)
if (inherits(f,"try-error")) ok <- FALSE
}
if (ok) {
ipiv[piv] <- 1:ncol(R)
R <- R[,ipiv]
break; ## success
}
}
if (i==20 && !ok) stop("Choleski based method failed, switch to QR")
## equivalent to qr update.
op <- options(warn = -1) ## otherwise warns if +ve semidef
R <- chol(XX,pivot=TRUE)
options(op)
p <- length(Xy)
ipiv <- piv <- attr(R,"pivot");ipiv[piv] <- 1:p
rank <- attr(R,"rank");ind <- 1:rank
if (rank<p) R[(rank+1):p,] <- 0 ## chol is buggy (R 3.1.0) with pivot=TRUE
f <- c(forwardsolve(t(R[ind,ind]),Xy[piv[ind]]),rep(0,p-rank))[ipiv]
R <- R[ipiv,ipiv]
list(R=R,f=f)
}
} ## chol2qr
qr.update <- function(Xn,yn,R=NULL,f=rep(0,0),y.norm2=0,use.chol=FALSE)
## Let X = QR and f = Q'y. This routine updates f and R
......
......@@ -595,7 +595,7 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
mu <- linkinv(eta)
}
trA <- oo$trA;
pearson <- sum((y-mu)^2/family$variance(mu)) ## Pearson statistic
pearson <- sum(weights*(y-mu)^2/family$variance(mu)) ## Pearson statistic
scale.est <- (pearson+dev.extra)/(n.true-trA)
......@@ -766,7 +766,8 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
gam.fit3.post.proc <- function(X,L,object) {
## get edf array and covariance matrices after a gam fit.
## X is original model matrix, L the mapping from working to full sp
Vb <- object$rV%*%t(object$rV)*object$scale.est ## Bayesian cov.
scale <- if (object$scale.estimated) object$scale.est else object$scale
Vb <- object$rV%*%t(object$rV)*scale ## Bayesian cov.
# PKt <- object$rV%*%t(object$K)
PKt <- .Call(C_mgcv_pmmult2,object$rV,object$K,0,1,object$control$nthreads)
# F <- PKt%*%(sqrt(object$weights)*X)
......@@ -797,7 +798,7 @@ gam.fit3.post.proc <- function(X,L,object) {
Vc <- crossprod(rV%*%t(object$db.drho))
Vc <- Vb + Vc ## Bayesian cov matrix with sp uncertainty
## finite sample size check on edf sanity...
edf2 <- rowSums(Vc*crossprod(R))/object$scale.est
edf2 <- rowSums(Vc*crossprod(R))/scale
if (sum(edf2)>sum(edf1)) edf2 <- edf1
} else edf2 <- Vc <- NULL
list(Vc=Vc,Vb=Vb,Ve=Ve,edf=edf,edf1=edf1,edf2=edf2,hat=hat,F=F,R=R)
......
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