Commit 4079654a authored by Dirk Eddelbuettel's avatar Dirk Eddelbuettel

Import Upstream version 1.8-6

parent 2ffccbf0
** denotes quite substantial/important changes
*** denotes really big changes
1.8-6
* Generalization of list formula handling to allow linear predictors to
share terms. e.g. gam(list(y1~s(x),y2~s(z),1+2~s(v)+w-1),family=mvn(d=2))
* New German translation thanks to Detlef Steuer.
* plot.gam now silently returns a list of plotting data, to help advanced
users (Fabian Scheipl) to produce custimized plot.
* bam can now set up an object suitable for fitting, but not actually do
the fit, following a suggestion by Fabian Scheipl. See arguments 'fit'
and 'G'.
1.8-5
* Korean translation added thanks to Chel Hee Lee.
* scale parameter handling in edf in logLik.gam made consistent with glm
(affects AIC).
......
Package: mgcv
Version: 1.8-5
Version: 1.8-6
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
......@@ -15,7 +15,7 @@ Suggests: splines, parallel, survival, MASS
LazyLoad: yes
ByteCompile: yes
License: GPL (>= 2)
Packaged: 2015-03-02 21:05:39 UTC; sw283
Packaged: 2015-03-30 13:02:31 UTC; sw283
NeedsCompilation: yes
Repository: CRAN
Date/Publication: 2015-03-03 09:18:09
Date/Publication: 2015-03-31 12:46:38
dcdf699f4b790dfe04b100f3c980436c *ChangeLog
cadc1146dc050c7dc00162207980193b *DESCRIPTION
2f31dc4982b7fa7c7f08ff3fad94b88d *ChangeLog
2244e729b443669c81a2732d005dd044 *DESCRIPTION
eb723b61539feef013de476e68b5c50a *GPL-2
af08e19c4214ac86d453bd90e1ce9e4e *NAMESPACE
4d8460a0933930189b287f8f77bdfabf *R/bam.r
4d14a8a6cdbe2e25634829667d3251ec *NAMESPACE
bb90db9c74f22d202abc4b0fc562448e *R/bam.r
52ea7081e235df890d3e51198bcc3092 *R/coxph.r
8e0796908b72ff7d4d5f57740050a488 *R/efam.r
a53c7592bb4c10a961f58e6537c87355 *R/efam.r
8cde8b6379034d9eb78caaf1b0d04cf8 *R/fast-REML.r
5df55d3cb82ac4669fc89909769d5a95 *R/gam.fit3.r
cbeb60506c33d463465efb5d1ada43f5 *R/gam.fit4.r
5b6cfb632772fb514428730b22066480 *R/gam.fit3.r
fc965ed3d5edf72c6f0335347c383dbb *R/gam.fit4.r
e63276b78a4ff2566af2e651bfc6aa9c *R/gam.sim.r
9a4ef9d8c7362256e314be0e65b5a96c *R/gamlss.r
633c048e36083646c7d409e99cf7b037 *R/gamlss.r
28ddc8624513070a62d0cc2aba7776a2 *R/gamm.r
09240f99d77e54848bf15d540c267022 *R/jagam.r
9a451b851f1e08fb7c7c60738556392d *R/mgcv.r
2f7ab8cc9f9cf2f6cb52fc72f9f69bfb *R/misc.r
6c53ea135a7ac956fbe8a3d16c59baaa *R/mvam.r
58413439e82b60d8942d070ca72f8e14 *R/plots.r
70eb07e6f5c4d61f6519863b8eebdfc0 *R/mgcv.r
bf5f3e43220b1e6db57c80caf66f1e54 *R/misc.r
66c24aed2f8fc9f6bce321794f8aff87 *R/mvam.r
f031621e351c06b12c328e4033e2f097 *R/plots.r
3e87a7ef880c8aa40ce57201f6ea5458 *R/smooth.r
68348617c5d3b0e2ea805c4360a8cdd4 *R/soap.r
76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
40874e3ced720a596750f499ded8a60a *data/columb.rda
074adce4131eced4cc71c67ad1d63c52 *inst/CITATION
5280bb00fa6595ff4a4b70d2e80114a8 *inst/po/de/LC_MESSAGES/R-mgcv.mo
1a78725a3c78f2b6f88e86f2a6efca1a *inst/po/de/LC_MESSAGES/R-mgcv.mo
c745969a1292eb3d49dfd0d0c2c997d4 *inst/po/de/LC_MESSAGES/mgcv.mo
9a49cc8d28e67874337cbe11d8632138 *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
447116b9728f9a82da5ce5d821850917 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
2d98abb49e987bb8ebdf9be5cdbbb68b *inst/po/fr/LC_MESSAGES/R-mgcv.mo
983aad6e8961a53c15c240d2ca586e35 *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
7aba5f2423e057971c95cdb30804af20 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
35f0d553e3181b070fed958a24dbf945 *inst/po/fr/LC_MESSAGES/R-mgcv.mo
643671acdb430cb0790f43559addbb0d *inst/po/fr/LC_MESSAGES/mgcv.mo
ea6f2b5c19372fda131171a03c983d6c *inst/po/ko/LC_MESSAGES/R-mgcv.mo
8a47c286e60655824529dead54f8368d *inst/po/ko/LC_MESSAGES/R-mgcv.mo
c1b1475e5fef49fe49929d2796ff87b6 *inst/po/ko/LC_MESSAGES/mgcv.mo
f8a930e0893c980fc2e9cfe276107e0f *inst/po/pl/LC_MESSAGES/R-mgcv.mo
cd7e6d1282796c089c320fbff388047f *inst/po/pl/LC_MESSAGES/R-mgcv.mo
715e52c0debf9848bbda15e94f5e7315 *inst/po/pl/LC_MESSAGES/mgcv.mo
251140c0e4b621fd578dae4ce8a28c8f *inst/po/po/LC_MESSAGES/R-mgcv.mo
6224353bb8556c6ae4106352e244a3ab *inst/po/po/LC_MESSAGES/R-mgcv.mo
4e7e8faae00111b7db61daab04358202 *inst/po/po/LC_MESSAGES/mgcv.mo
c574fe1ca9d55a9818d308906f16d16e *man/Beta.Rd
5bf12ddc0dab9daae72271b96a15c539 *man/Predict.matrix.Rd
......@@ -41,7 +41,7 @@ e9b0a2e31b130cf2cb38618ec50d1919 *man/Predict.matrix.soap.film.Rd
f12468625253dd3603907de233762fd6 *man/Rrank.Rd
f6cadda5999c35800fd65c08d6812f7b *man/Tweedie.Rd
80f8763baa4987579e2aa56073a9e94e *man/anova.gam.Rd
ef84108378322022dfe6ca5d5c9c917a *man/bam.Rd
de9b8ca4369a4fc7ba69c68c1cdb98c2 *man/bam.Rd
71bde8b8caa24a36529ce7e0ac3165d8 *man/bam.update.Rd
a2beb811b1093c5e82ef32d7de1f7d32 *man/cSplineDes.Rd
9b4d616d1b6c4a46ca77d16cded3f806 *man/choose.k.Rd
......@@ -54,10 +54,10 @@ db8f5192f7cce5e2ed07fb5bd97b580f *man/family.mgcv.Rd
44ad0563add1c560027d502ce41483f5 *man/fix.family.link.Rd
4d4eea9ad2b78e765a96b2a0065725c1 *man/fixDependence.Rd
c7561a0f2e114247955e212aeabc6ae9 *man/formXtViX.Rd
d87ff6287d7e343ea24d2053f4724b35 *man/formula.gam.Rd
c6fd48d86959402982c75566876baa16 *man/formula.gam.Rd
4da4d585b329769eb44f0c7a6e7dd554 *man/fs.test.Rd
6f405acde2d7b6f464cf45f5395113ba *man/full.score.Rd
a54ea44bb9641389bcbd21ce0e578d8f *man/gam.Rd
86678ba1579e9deb8ac3337af60a9c20 *man/gam.Rd
fe61dd0efab9e920c17335faf3d5764c *man/gam.check.Rd
a65bc22f606e45d185bc375fbf5698a1 *man/gam.control.Rd
44db24b66ce63bc16d2c8bc3f5b42ac5 *man/gam.convergence.Rd
......@@ -93,7 +93,7 @@ c00bcfe2d0b44b2ea955f3934421807c *man/interpret.gam.Rd
00ccf213c31910cd14f1df65a300eb33 *man/model.matrix.gam.Rd
bc9b89db7e7ff246749551c16f5f1f07 *man/mono.con.Rd
d33914a328f645af13f5a42914ca0f35 *man/mroot.Rd
2423b710038a5d53d7d28dc307f33642 *man/mvn.Rd
d0c2b73407680aac50a17753346d814f *man/mvn.Rd
9e95f6b9eec476188de830beb9ff3e0b *man/negbin.Rd
8b766a6ad848b0f1ca469e381ded0169 *man/new.name.Rd
7d8f62e182f7c428cc9d46ddd4d97d43 *man/notExp.Rd
......@@ -105,7 +105,7 @@ c64b44bc684cbadfe77e7f6dc833ddb4 *man/pcls.Rd
8c0f8575b427f30316b639a326193aeb *man/pdTens.Rd
1721f1b266d9e14827e8226e2cb74a81 *man/pen.edf.Rd
edf57071572275a8443b2f0b66d44424 *man/place.knots.Rd
c488918033a996ce97943bf20ff8ac54 *man/plot.gam.Rd
04a9af3e294bd14ad590aacb6d26332d *man/plot.gam.Rd
c27a6b886929b1dc83bf4b90cae848f9 *man/polys.plot.Rd
085dff270e5255315cb244368faf79df *man/predict.bam.Rd
900647986ab1365d437368cc28b62658 *man/predict.gam.Rd
......@@ -115,6 +115,7 @@ f77ca1471881d2f93c74864d076c0a0e *man/rTweedie.Rd
5ff3bd5034e8f2afa3174c0a2d989d97 *man/random.effects.Rd
c523210ae95cb9aaa0aaa1c37da1a4c5 *man/residuals.gam.Rd
f6f1333be2587ffef5970905b13468ea *man/rig.Rd
9f6f46f5c5da080bc82f9aa4685d364a *man/rmvn.Rd
7258dfc0149fff020054117fd2ee6bd8 *man/s.Rd
c438eb3e41cb1f51e72283380145d210 *man/scat.Rd
2199afd400a1821b74a24e2578bc0224 *man/single.index.Rd
......@@ -146,18 +147,18 @@ a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
2c5f6815e609f2cdb56b0067a183f915 *man/ziP.Rd
ae8388103d8b1be39f55f426b205b576 *man/ziplss.Rd
a7c1582252be4e0013e4275ffd3aac4c *po/R-de.po
4550c3e06d76e6e8a9b22fb8bbd76eb3 *po/R-de.po
0bdfcf98961b0d52b60f806dc1dba77e *po/R-en@quot.po
eb3707a62b0e3d406ddf5c15bc40c673 *po/R-fr.po
c7317bc4fb746b3ce2a9bcdf9d928716 *po/R-ko.po
5e7691372a93b11f3c0efa48ec9dae2f *po/R-mgcv.pot
62fa356ec68d7cc8a9ec1b3915b0db13 *po/R-pl.po
a5335d43184be0f9da9df45cfadc3e33 *po/R-po.po
4e65e93fef4d034a399f90421e8f323a *po/R-fr.po
73cdaf7a5a69f0b7cbfe411cd0c468b6 *po/R-ko.po
7eb472ce4c2d1dc30b4dd1091c0e88df *po/R-mgcv.pot
7b07899266c3acf3d2a625850d7cd6ef *po/R-pl.po
5b91cecd0b9e52154185d380a273c623 *po/R-po.po
ccf3140169b68ec7aff4b15e6a97e5db *po/de.po
93f72334356fe6f05a64e567efd35c8e *po/en@quot.po
fb829b82760779929951d49fe29ed2e5 *po/fr.po
dc1ef92ff4454734c3a24876e299b760 *po/ko.po
5fee1431cf41b342cde9e5eb6dd27cab *po/mgcv.pot
8ad4757e026d1841c8f43eb97072c06e *po/mgcv.pot
dfd4eec9edc7d1ab6354d47b6e2bd42f *po/pl.po
b56dac4547037ea45e2c8f9bce7aa9ef *po/po.po
03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
......
......@@ -38,7 +38,7 @@ export(anova.gam, bam, bam.update, betar, cox.ph,concurvity, cSplineDes,
Predict.matrix.random.effect,
Predict.matrix.t2.smooth,
qq.gam,
residuals.gam,rig,rTweedie,
residuals.gam,rig,rTweedie,rmvn,
Rrank,s,scat,
sim2jam,
slanczos,
......
This diff is collapsed.
......@@ -950,7 +950,8 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
theta <- exp(theta) ## note log theta supplied
muth <- mu*theta
## yth <- y*theta
2* wt * (-lgamma(theta) +lgamma(muth) + lgamma(theta - muth) - muth*log(y/(1-y)) - theta*log(1-y) + log(y*(1-y)))
2* wt * (-lgamma(theta) +lgamma(muth) + lgamma(theta - muth) - muth*(log(y)-log1p(-y)) -
theta*log1p(-y) + log(y) + log1p(-y))
}
Dd <- function(y, mu, theta, wt, level=0) {
......@@ -970,7 +971,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
psi2.onemuth <- psigamma(onemuth,2)
psi3.muth <- psigamma(muth,3)
psi3.onemuth <- psigamma(onemuth,3)
log.yoney <- log(y)-log(oney)
log.yoney <- log(y)-log1p(-y)
r <- list()
## get the quantities needed for IRLS.
## Dmu2eta2 is deriv of D w.r.t mu twice and eta twice,
......@@ -979,7 +980,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
r$Dmu2 <- 2 * wt * theta^2*(psi1.muth+psi1.onemuth)
r$EDmu2 <- r$Dmu2
if (level>0) { ## quantities needed for first derivatives
r$Dth <- 2 * wt *theta*(-mu*log.yoney - log(oney)+ mu*psi0.muth+onemu*psi0.onemuth -psi0.th)
r$Dth <- 2 * wt *theta*(-mu*log.yoney - log1p(-y)+ mu*psi0.muth+onemu*psi0.onemuth -psi0.th)
r$Dmuth <- r$Dmu + 2 * wt * theta^2*(mu*psi1.muth -onemu*psi1.onemuth)
r$Dmu3 <- 2 * wt *theta^3 * (psi2.muth - psi2.onemuth)
r$Dmu2th <- 2* r$Dmu2 + 2 * wt * theta^3* (mu*psi2.muth + onemu*psi2.onemuth)
......@@ -1001,7 +1002,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
theta <- exp(theta)
muth <- mu*theta
term <- -lgamma(theta)+lgamma(muth)+lgamma(theta-muth)-(muth-1)*log(y)-
(theta-muth-1)*log(1-y) ## `-' log likelihood for each observation
(theta-muth-1)*log1p(-y) ## `-' log likelihood for each observation
2 * sum(term * wt)
}
......@@ -1042,7 +1043,7 @@ betar <- function (theta = NULL, link = "logit",eps=.Machine$double.eps*100) {
mu[!ind] <- (a + b*expeta[!ind])/(1+expeta[!ind])
l <- dbeta(y,phi*mu,phi*(1-mu),log=TRUE)
if (deriv) {
g <- phi * log(y) - phi * log(1-y) - phi * digamma(mu*phi) + phi * digamma((1-mu)*phi)
g <- phi * log(y) - phi * log1p(-y) - phi * digamma(mu*phi) + phi * digamma((1-mu)*phi)
h <- -phi^2*(trigamma(mu*phi)+trigamma((1-mu)*phi))
dmueta2 <- dmueta1 <- eta
dmueta1 <- expeta*(b-a)/(1+expeta)^2
......@@ -1235,7 +1236,7 @@ scat <- function (theta = NULL, link = "identity") {
dev.resids <- function(y, mu, wt,theta=NULL) {
if (is.null(theta)) theta <- get(".Theta")
nu <- exp(theta[1])+2; sig <- exp(theta[2])
wt * (nu + 1)*log(1+(1/nu)*((y-mu)/sig)^2)
wt * (nu + 1)*log1p((1/nu)*((y-mu)/sig)^2)
}
Dd <- function(y, mu, theta, wt, level=0) {
......@@ -1312,7 +1313,7 @@ scat <- function (theta = NULL, link = "identity") {
if (is.null(theta)) theta <- get(".Theta")
nu <- exp(theta[1])+2; sig <- exp(theta[2])
term <- -lgamma((nu+1)/2)+ lgamma(nu/2) + log(sig*(pi*nu)^.5) +
(nu+1)*log(1 + ((y-mu)/sig)^2/nu)/2 ## `-'log likelihood for each observation
(nu+1)*log1p(((y-mu)/sig)^2/nu)/2 ## `-'log likelihood for each observation
2 * sum(term * wt)
}
......
......@@ -2012,14 +2012,14 @@ fix.family.link.extended.family <- function(fam) {
} ## cauchit
if (link == "cloglog") {
## g = log(-log(1-mu)), g' = -1/(log(1-mu)*(1-mu))
fam$g2g <- function(mu) { l1m <- log(1-mu)
fam$g2g <- function(mu) { l1m <- log1p(-mu)
-l1m - 1
}
fam$g3g <- function(mu) { l1m <- log(1-mu)
fam$g3g <- function(mu) { l1m <- log1p(-mu)
l1m*(2*l1m + 3) + 2
}
fam$g4g <- function(mu){
l1m <- log(1-mu)
l1m <- log1p(-mu)
-l1m*(l1m*(6*l1m+11)+12)-6
}
return(fam)
......@@ -2096,15 +2096,15 @@ fix.family.link.family <- function(fam)
}
if (link == "cloglog") {
## g = log(-log(1-mu)), g' = -1/(log(1-mu)*(1-mu))
fam$d2link <- function(mu) { l1m <- log(1-mu)
fam$d2link <- function(mu) { l1m <- log1p(-mu)
-1/((1 - mu)^2*l1m) *(1+ 1/l1m)
}
fam$d3link <- function(mu) { l1m <- log(1-mu)
fam$d3link <- function(mu) { l1m <- log1p(-mu)
mu3 <- (1-mu)^3
(-2 - 3*l1m - 2*l1m^2)/mu3/l1m^3
}
fam$d4link <- function(mu){
l1m <- log(1-mu)
l1m <- log1p(-mu)
mu4 <- (1-mu)^4
( - 12 - 11 * l1m - 6 * l1m^2 - 6/l1m )/mu4 /l1m^3
}
......
......@@ -694,7 +694,6 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
if (!is.null(start0)) start <- start0
coef <- as.numeric(start)
if (is.null(weights)) weights <- rep.int(1, nobs)
if (is.null(offset)) offset <- rep.int(0, nobs)
......@@ -821,14 +820,18 @@ gam.fit5 <- function(x,y,lsp,Sl,weights=NULL,offset=NULL,deriv=2,family,
St <- St[-drop,-drop]
x <- x[,-drop] ## dropping columns from model matrix
if (!is.null(lpi)) { ## need to adjust column indexes as well
k <- 0
ii <- (1:q)[!bdrop];ij <- rep(NA,q)
ij[ii] <- 1:length(ii) ## col i of old model matrix is col ij[i] of new
#k <- 0
for (i in 1:length(lpi)) {
kk <- sum(lpi[[i]]%in%drop==FALSE) ## how many left undropped?
lpi[[i]] <- 1:kk + k ## new index - note strong assumptions on structure here
k <- k + kk
#kk <- sum(lpi[[i]]%in%drop==FALSE) ## how many left undropped?
#lpi[[i]] <- 1:kk + k ## new index - note strong assumptions on structure here
#k <- k + kk
lpi[[i]] <- ij[lpi[[i]][!(lpi[[i]]%in%drop)]] # drop and shuffle up
}
} ## lpi adjustment done
attr(x,"lpi") <- lpi
attr(x,"drop") <- drop ## useful if family has precomputed something from x
ll <- family$ll(y,x,coef,weights,family,deriv=1)
ll0 <- ll$l - (t(coef)%*%St%*%coef)/2
}
......@@ -1126,7 +1129,7 @@ deriv.check5 <- function(x, y, sp,
{ if (!deriv%in%c(1,2)) stop("deriv should be 1 or 2")
if (control$epsilon>1e-9) control$epsilon <- 1e-9
## first obtain the fit corresponding to sp...
b <- gam.fit5(x=x,y=y,lsp=sp,Sl=Sl,weights=weights,offset=offset,deriv=2,
b <- gam.fit5(x=x,y=y,lsp=sp,Sl=Sl,weights=weights,offset=offset,deriv=deriv,
family=family,control=control,Mp=Mp,start=start)
## now get the derivatives of the likelihood w.r.t. coefs...
ll <- family$ll(y=y,X=x,coef=b$coefficients,wt=weights,family=family,
......
......@@ -175,7 +175,7 @@ gamlss.etamu <- function(l1,l2,l3=NULL,l4=NULL,ig1,g2,g3=NULL,g4=NULL,i2,i3=NULL
} # gamlss.etamu
gamlss.gH <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) {
gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) {
## X[,jj[[i]]] is the ith model matrix.
## lj contains jth derivatives of the likelihood for each datum,
## columns are w.r.t. different combinations of parameters.
......@@ -229,7 +229,7 @@ gamlss.gH <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=N
v <- v + l3[,i3[i,i,q]] * d1eta[ind,l]
ind <- ind + n
}
d1H[jj[[j]],l] <- colSums(X[,jj[[i]],drop=FALSE]*(v*X[,jj[[i]],drop=FALSE]))
d1H[jj[[i]],l] <- colSums(X[,jj[[i]],drop=FALSE]*(v*X[,jj[[i]],drop=FALSE]))
}
}
} ## if deriv==1
......@@ -304,6 +304,142 @@ gamlss.gH <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=N
}
} ## if deriv>2
list(lb=lb,lbb=lbb,d1H=d1H,d2H=d2H,trHid2H=trHid2H)
} ## end of gamlss.gH0
gamlss.gH <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=NULL,D=NULL) {
## X[,jj[[i]]] is the ith model matrix.
## lj contains jth derivatives of the likelihood for each datum,
## columns are w.r.t. different combinations of parameters.
## ij is the symmetric array indexer for the jth order derivs...
## e.g. l4[,i4[i,j,l,m]] contains derivatives with
## respect to parameters indexed by i,j,l,m
## d1b and d2b are first and second derivatives of beta w.r.t. sps.
## fh is a factorization of the penalized hessian, while D contains the corresponding
## Diagonal pre-conditioning weights.
## deriv: 0 - just grad and Hess
## 1 - diagonal of first deriv of Hess
## 2 - first deriv of Hess
## 3 - everything.
K <- length(jj)
p <- ncol(X);n <- nrow(X)
trHid2H <- d1H <- d2H <- NULL ## defaults
## the gradient...
lb <- rep(0,p)
for (i in 1:K) { ## first derivative loop
lb[jj[[i]]] <- lb[jj[[i]]] + colSums(l1[,i]*X[,jj[[i]],drop=FALSE]) ## !
}
## the Hessian...
lbb <- matrix(0,p,p)
for (i in 1:K) for (j in i:K) {
A <- t(X[,jj[[i]],drop=FALSE])%*%(l2[,i2[i,j]]*X[,jj[[j]],drop=FALSE])
lbb[jj[[i]],jj[[j]]] <- lbb[jj[[i]],jj[[j]]] + A
if (j>i) lbb[jj[[j]],jj[[i]]] <- lbb[jj[[j]],jj[[i]]] + t(A)
}
if (deriv>0) {
## the first derivative of the Hessian, using d1b
## the first derivates of the coefficients wrt the sps
m <- ncol(d1b) ## number of smoothing parameters
## stack the derivatives of the various linear predictors on top
## of each other...
d1eta <- matrix(0,n*K,m)
ind <- 1:n
for (i in 1:K) {
d1eta[ind,] <- X[,jj[[i]],drop=FALSE]%*%d1b[jj[[i]],]
ind <- ind + n
}
}
if (deriv==1) {
d1H <- matrix(0,p,m) ## only store diagonals of d1H
for (l in 1:m) {
for (i in 1:K) {
v <- rep(0,n);ind <- 1:n
for (q in 1:K) {
v <- v + l3[,i3[i,i,q]] * d1eta[ind,l]
ind <- ind + n
}
d1H[jj[[i]],l] <- d1H[jj[[i]],l] + colSums(X[,jj[[i]],drop=FALSE]*(v*X[,jj[[i]],drop=FALSE]))
}
}
} ## if deriv==1
if (deriv>1) {
d1H <- list()
for (l in 1:m) {
d1H[[l]] <- matrix(0,p,p)
for (i in 1:K) for (j in i:K) {
v <- rep(0,n);ind <- 1:n
for (q in 1:K) {
v <- v + l3[,i3[i,j,q]] * d1eta[ind,l]
ind <- ind + n
}
## d1H[[l]][jj[[j]],jj[[i]]] <-
A <- t(X[,jj[[i]],drop=FALSE])%*%(v*X[,jj[[j]],drop=FALSE])
d1H[[l]][jj[[i]],jj[[j]]] <- d1H[[l]][jj[[i]],jj[[j]]] + A
if (j>i) d1H[[l]][jj[[j]],jj[[i]]] <- d1H[[l]][jj[[j]],jj[[i]]] + t(A)
}
}
} ## if deriv>1
if (deriv>2) {
## need tr(Hp^{-1} d^2H/drho_k drho_j)
## First form the expanded model matrix...
VX <- Xe <- matrix(0,K*n,ncol(X))
ind <- 1:n
for (i in 1:K) {
Xe[ind,jj[[i]]] <- X[,jj[[i]]]
ind <- ind + n
}
## Now form Hp^{-1} Xe'...
if (is.list(fh)) { ## then the supplied factor is an eigen-decomposition
d <- fh$values;d[d>0] <- 1/d[d>0];d[d<=0] <- 0
Xe <- t(D*((fh$vectors%*%(d*t(fh$vectors)))%*%(D*t(Xe))))
} else { ## the supplied factor is a choleski factor
ipiv <- piv <- attr(fh,"pivot");ipiv[piv] <- 1:p
Xe <- t(D*(backsolve(fh,forwardsolve(t(fh),(D*t(Xe))[piv,]))[ipiv,]))
}
## now compute the required trace terms
d2eta <- matrix(0,n*K,ncol(d2b))
ind <- 1:n
for (i in 1:K) {
d2eta[ind,] <- X[,jj[[i]],drop=FALSE]%*%d2b[jj[[i]],]
ind <- ind + n
}
trHid2H <- rep(0,ncol(d2b))
kk <- 0 ## counter for second derivatives
for (k in 1:m) for (l in k:m) { ## looping over smoothing parameters...
kk <- kk + 1
for (i in 1:K) for (j in 1:K) {
v <- rep(0,n);ind <- 1:n
for (q in 1:K) { ## accumulate the diagonal matrix for X_i'diag(v)X_j
v <- v + d2eta[ind,kk]*l3[,i3[i,j,q]]
ins <- 1:n
for (s in 1:K) {
v <- v + d1eta[ind,k]*d1eta[ins,l]*l4[,i4[i,j,q,s]]
ins <- ins + n
}
ind <- ind + n
}
if (i==j) {
rind <- 1:n + (i-1)*n
VX[rind,jj[[i]]] <- v * X[,jj[[i]]]
} else {
rind1 <- 1:n + (i-1)*n
rind2 <- 1:n + (j-1)*n
VX[rind2,jj[[i]]] <- v * X[,jj[[i]]]
VX[rind1,jj[[j]]] <- v * X[,jj[[j]]]
}
}
trHid2H[kk] <- sum(Xe*VX)
}
} ## if deriv>2
list(lb=lb,lbb=lbb,d1H=d1H,d2H=d2H,trHid2H=trHid2H)
} ## end of gamlss.gH
......
This diff is collapsed.
......@@ -2,6 +2,22 @@
## Many of the following are simple wrappers for C functions, used largely
## for testing purposes
rmvn <- function(n,mu,V) {
## generate multivariate normal deviates. e.g.
## V <- matrix(c(2,1,1,2),2,2); mu <- c(1,1);n <- 1000;z <- rmvn(n,mu,V);crossprod(sweep(z,2,colMeans(z)))/n
p <- ncol(V)
R <- mroot(V,rank=ncol(V)) ## RR' = V
if (is.matrix(mu)) {
if (ncol(mu)!=p||nrow(mu)!=n) stop("mu dimensions wrong")
z <- matrix(rnorm(p*n),n,p)%*%t(R) + mu
} else {
if (length(mu)!=p) stop("mu dimensions wrong")
z <- t(R%*% matrix(rnorm(p*n),p,n) + mu)
if (n==1) z <- as.numeric(z)
}
z
} ## rmvn
mgcv.omp <- function() {
## does open MP appear to be available?
oo <- .C(C_mgcv_omp,a=as.integer(-1))
......
## (c) Simon N. Wood (2013, 2014) mvn model extended family.
## (c) Simon N. Wood (2013-2015) mvn model extended family.
## Released under GPL2 ...
lpi.expand <- function(X,trailing=TRUE) {
## takes a model matrix X, with "lpi" attribute, and produces
## full redundant version in which each column block is the full
## model matrix for one linear predictor, which may involve
## repeating columns between blocks.
## See mvn family (ll) for prototypic application
lpi <- attr(X,"lpi")
if (!attr(lpi,"overlap")) return(X) ## nothing to do
ip <- unlist(lpi)
if (trailing&&max(ip)<ncol(X)) { ## include any unindexed trailing blocks
ii <- (max(ip)+1):ncol(X)
X <- cbind(X[,ip],X[,ii])
} else X <- X[,ip]
k <- 0
for (i in 1:length(lpi)) {
lpi[[i]] <- 1:length(lpi[[i]]) + k
k <- k + length(lpi[[i]])
}
attr(X,"lpi") <- lpi
X
} ## lpi.expand
lpi.contract <- function(x,lpi,type="rc",trailing=TRUE) {
## takes a vector or matrix x, and applies an lpi contraction to it
## if x is a matrix then type can be "r", "c" or "rc" for row, col
## or row, column contraction.
## See mvn family (ll) for prototypic application
ip <- unlist(lpi)
if (trailing) { ## copy across any un-indexed trailing blocks
lip <- length(ip) ## last row/col indexed in x
## get the length, nt, of any unindexed trailing block
nt <- 0
if (is.matrix(x)) {
if (type=="r"&&nrow(x)>lip) nt <- nrow(x)-lip else
if (ncol(x)>lip) nt <- ncol(x) - lip
} else if (length(x)>lip) nt <- length(x) - lip
if (nt>0) { ## there is a trailing block - index it in lpi
lpi[[length(lpi)+1]] <- 1:nt + max(ip)
ip <- unlist(lpi)
}
}
p <- max(ip) ## dimension of result
if (is.matrix(x)) {
if (type=="c"||type=="rc") { ## column contraction
k <- 0
z <- matrix(0,nrow(x),p)
for (i in 1:length(lpi)) {
ii <- 1:length(lpi[[i]]) + k
k <- k + length(ii)
z[,lpi[[i]]] <- z[,lpi[[i]]] + x[,ii]
}
if (type=="rc") x <- z
}
if (type=="r"||type=="rc") { ## row contraction
z <- matrix(0,p,ncol(x))
k <- 0
for (i in 1:length(lpi)) {
ii <- 1:length(lpi[[i]]) + k
k <- k + length(ii)
z[lpi[[i]],] <- z[lpi[[i]],] + x[ii,]
}
}
} else { ## vector
z <- rep(0,p);k <- 0
for (i in 1:length(lpi)) {
ii <- 1:length(lpi[[i]]) + k
k <- k + length(ii)
z[lpi[[i]]] <- z[lpi[[i]]] + x[ii]
}
}
z
} ## lpi.contract
mvn <- function(d=2) {
## Extended family object for multivariate normal additive model.
if (d<2) stop("mvn requires 2 or more dimensional data")
......@@ -9,7 +82,7 @@ mvn <- function(d=2) {
stats[[i]] <- make.link("identity")
}
env <- new.env(parent = .GlobalEnv)
##env <- new.env(parent = .GlobalEnv)
validmu <- function(mu) all(is.finite(mu))
## initialization has to add in the extra parameters of
......@@ -21,7 +94,7 @@ mvn <- function(d=2) {
## finds initial coefficients
ydim <- ncol(G$y) ## dimension of response
nbeta <- ncol(G$X)
ntheta <- ydim*(ydim+1)/2 ## numer of cov matrix factor params
ntheta <- ydim*(ydim+1)/2 ## number of cov matrix factor params
lpi <- attr(G$X,"lpi")
XX <- crossprod(G$X)
G$X <- cbind(G$X,matrix(0,nrow(G$X),ntheta)) ## add dummy columns to G$X
......@@ -36,8 +109,10 @@ mvn <- function(d=2) {
## now get initial parameters and store in family...
for (k in 1:ydim) {
sin <- G$off %in% lpi[[k]]
Sk <- G$S[sin]
um <- magic(G$y[,k],G$X[,lpi[[k]]],rep(-1,sum(sin)),G$S[sin],G$off[sin]-lpi[[k]][1]+1,nt=control$nthreads)
#Sk <- G$S[sin]
um <- magic(G$y[,k],G$X[,lpi[[k]]],rep(-1,sum(sin)),G$S[sin],
match(G$off[sin],lpi[[k]]), ## turn G$off global indices into indices for this predictor
nt=control$nthreads)
G$family$ibeta[lpi[[k]]] <- um$b
G$family$ibeta[nbeta+1] <- -.5*log(um$scale) ## initial log root precision
nbeta <- nbeta + ydim - k + 1
......@@ -64,8 +139,6 @@ mvn <- function(d=2) {
initialize <- expression({
## called in gam.fit5 and initial.spg
## Ideally fit separate models to each component and
## extract initial coefs, s.p.s and variances this way
n <- rep(1, nobs)
if (is.null(start)) start <- family$ibeta
## need to re-parameterize XX is non-standard
......@@ -84,7 +157,7 @@ mvn <- function(d=2) {
##rd <- qf <- NULL ## these functions currently undefined for
ll <- function(y,X,coef,wt,family,deriv=0,d1b=0,d2b=0,Hp=NULL,rank=0,fh=NULL,D=NULL) {
ll <- function(y,X,coef,wt,family,deriv=0,d1b=NULL,d2b=NULL,Hp=NULL,rank=0,fh=NULL,D=NULL) {
## function defining the Multivariate Normal model log lik.
## Calls C code "mvn_ll"
## deriv codes: 0 - eval; 1 - grad and Hessian
......@@ -97,11 +170,31 @@ mvn <- function(d=2) {
## D is the diagonal pre-conditioning matrix used to obtain Hp
## if Hr is the raw Hp then Hp = D*t(D*Hr)
lpi <- attr(X,"lpi") ## lpi[[k]] is index of model matrix columns for dim k
m <- length(lpi) ## number of dimensions of MVN
overlap <- attr(lpi,"overlap") ## do dimensions share terms?
drop <- attr(X,"drop")
if (!is.null(drop)) { ## the optimizer has dropped some parameters
## it will have adjusted lpi automatically, but XX is mvn specific
attr(X,"XX") <- attr(X,"XX")[-drop,-drop]
}
m <- length(lpi) ## number of dimensions of MVN
if (overlap) { ## linear predictors share terms - expand to redundant representation
ip <- unlist(lpi)
XX <- attr(X,"XX")[ip,ip]
X <- lpi.expand(X)