version 1.8-17

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: mgcv
-Version: 1.8-16
+Version: 1.8-17
 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
@@ -16,6 +16,6 @@ LazyLoad: yes
 ByteCompile: yes
 License: GPL (>= 2)
 NeedsCompilation: yes
-Packaged: 2016-11-07 10:22:00 UTC; sw283
+Packaged: 2017-02-06 11:03:57 UTC; sw283
 Repository: CRAN
-Date/Publication: 2016-11-07 19:28:16
+Date/Publication: 2017-02-08 21:30:16

MD5

@@ -1,36 +1,36 @@
-a38e3e6c9f0ae98bdcf3d14127347fde *ChangeLog
-e051228d6327e0b3d640616d4fa510c3 *DESCRIPTION
+fdc45fe0bc09d3a9dea994040574ba9d *ChangeLog
+1671dbe5a90499d24659b13d04370628 *DESCRIPTION
 eb723b61539feef013de476e68b5c50a *GPL-2
-595f7fa74dd678ae4be1f6968e174555 *NAMESPACE
+fe1f2bde9713ba97eeb59d928310f628 *NAMESPACE
 1aee82e17a9feb6166dbe107659029c9 *R/bam.r
 87933a9e2845f871334d96f537ee11e9 *R/coxph.r
-a25f5a55d9d8b4ada56ac0df3ec99a18 *R/efam.r
-4b370253beb1eda19d389132ade30212 *R/fast-REML.r
-cb165a0a11de9686910fb27b0d6450ec *R/gam.fit3.r
-7734ba62022328a50b771e721df98a72 *R/gam.fit4.r
+97c5e0261e26e0aedd8127f9ed41f817 *R/efam.r
+deaf0c6a0e39d42a9bad487772eb5904 *R/fast-REML.r
+dbf376ceeec4b343d15daa90fe5a81fb *R/gam.fit3.r
+9f73b1e3e9781d8ee26a46c1a2127c73 *R/gam.fit4.r
 1b620b840ca67c9759639bd6e53824e3 *R/gam.sim.r
-66257913c556f135657b7c12b6ed733d *R/gamlss.r
+a6b6a5e70c7889370d33876e08396a09 *R/gamlss.r
 cceac26b3d01513f8f87eacae91ddae0 *R/gamm.r
 10facb791e4cfd123d183f05660119c6 *R/jagam.r
-8cebce4d68d4c45d661073756ca7b202 *R/mgcv.r
+088976e7d2ee3df696d09b3fcdfdd563 *R/mgcv.r
 2feca0dc9d354de7bc707c67a5891364 *R/misc.r
 03772998ab05887f2eeae10dd6efe983 *R/mvam.r
 5d07d6e4f75a64b6d04cda9e23278d70 *R/plots.r
-cf3543d1d4b7fe6fb8576f681b345722 *R/smooth.r
+9a4db620ad080f012b56134c1aba10bb *R/smooth.r
 666d7fd36fda68b928993d5388b0d7fc *R/soap.r
 76cc875719bf0ef9eab45ea5bfeccda6 *R/sparse.r
 e468195a83fab90da8e760c2c3884bd3 *data/columb.polys.rda
 40874e3ced720a596750f499ded8a60a *data/columb.rda
-074adce4131eced4cc71c67ad1d63c52 *inst/CITATION
-7279083a0b1eb3505c3e8105c7c5d2e7 *inst/po/de/LC_MESSAGES/R-mgcv.mo
+e74e6d792e3dc69f7a20c68c16097145 *inst/CITATION
+342b1840b6a569b8d8d4a3355a0d2bfa *inst/po/de/LC_MESSAGES/R-mgcv.mo
 fe9d11e3087789da149e3688309df670 *inst/po/de/LC_MESSAGES/mgcv.mo
-14104d65bd4a0005bbc5a3932a37ae8f *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
-6b8ae16fa2affd2676d2a47c61885dd1 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
-a78bad32fe0b74e954a15ccd98f3a995 *inst/po/fr/LC_MESSAGES/R-mgcv.mo
+72ccd181324c1acd4671f69311ed494a *inst/po/en@quot/LC_MESSAGES/R-mgcv.mo
+dcb794aad0b69ab80cb6877510c3a925 *inst/po/en@quot/LC_MESSAGES/mgcv.mo
+72ffab2a2903e1efa30cc81424bce43d *inst/po/fr/LC_MESSAGES/R-mgcv.mo
 418bef2f2c1ed07bff6bbdb6884d2858 *inst/po/fr/LC_MESSAGES/mgcv.mo
-1553db580cf9e89f3876fe888f9e4e94 *inst/po/ko/LC_MESSAGES/R-mgcv.mo
+fbe9f7d7c30fbcebf97c28bf17c2ef50 *inst/po/ko/LC_MESSAGES/R-mgcv.mo
 e6196a86ad3a8e42df84242a861d362a *inst/po/ko/LC_MESSAGES/mgcv.mo
-6291a3775233a0315dd626e50fc89f2f *inst/po/pl/LC_MESSAGES/R-mgcv.mo
+5c8c8783edc20ee9196bb17bb86c420d *inst/po/pl/LC_MESSAGES/R-mgcv.mo
 07e822258166c032ff9f1a4e96025841 *inst/po/pl/LC_MESSAGES/mgcv.mo
 c574fe1ca9d55a9818d308906f16d16e *man/Beta.Rd
 5bf12ddc0dab9daae72271b96a15c539 *man/Predict.matrix.Rd
@@ -73,6 +73,8 @@ b2ff252537dd2155524b774b2435e66e *man/gam.vcomp.Rd
 eb8648cc6b3c9374b899abd2f2b12f7b *man/gam2objective.Rd
 717401fd7efa3b39d90418a5d1d0c216 *man/gamObject.Rd
 a2593990d6a87f7b783d0435062dec02 *man/gamSim.Rd
+e5d2541f32dab56972f58b0773eba50c *man/gamlss.etamu.Rd
+c7f140d128d1d1d76909499900faf49e *man/gamlss.gH.Rd
 351009f7345aaff40f3377e77f4ce7ad *man/gamm.Rd
 ec2d1b2aa87cc3f8424e07b6dc0340d5 *man/gaulss.Rd
 398a5c12285401c1d37a8edb58780bc3 *man/get.var.Rd
@@ -83,7 +85,7 @@ faf325db054dce698286f87e09ddc202 *man/gevlss.Rd
 2f222eeeb3d7bc42f93869bf8c2af58a *man/influence.gam.Rd
 39b9de9dbac7d9dc5c849e1a37def675 *man/initial.sp.Rd
 2a37ae59a9f9f5a0a58af45947eca524 *man/interpret.gam.Rd
-557a6b7a01d8c0ca7d4db3d9cf329ce6 *man/jagam.Rd
+29f94d99b9a3decd4d71618eb839a7f7 *man/jagam.Rd
 87d942b17bee05bb662270b894b183e6 *man/ldTweedie.Rd
 58e73ac26b93dc9d28bb27c8699e12cf *man/linear.functional.terms.Rd
 93035193b0faa32700e1421ce8c1e9f6 *man/logLik.gam.Rd
@@ -129,7 +131,7 @@ d515e51ec98d73af6166f7b31aeaba9b *man/scat.Rd
 6f03e337d54221bc167d531e25af1eea *man/slanczos.Rd
 8020154bd5c709d11f0e7cf043df886d *man/smooth.construct.Rd
 4a689eba97e4fed138dccb8cad13205e *man/smooth.construct.ad.smooth.spec.Rd
-9590099266d7f6fafd6bdd1e5c3f61da *man/smooth.construct.bs.smooth.spec.Rd
+24927a9c5ce97da988ac8ff3299c533a *man/smooth.construct.bs.smooth.spec.Rd
 76013feaf70d00976bba0154b6f2c946 *man/smooth.construct.cr.smooth.spec.Rd
 f5e6d0f5122f61c336827b3615482157 *man/smooth.construct.ds.smooth.spec.Rd
 db75c958cbfb561914a3291ab58b9786 *man/smooth.construct.fs.smooth.spec.Rd
@@ -141,7 +143,7 @@ b45d8e71bda4ceca8203dffea577e441 *man/smooth.construct.re.smooth.spec.Rd
 0bfe981f2c3e6ea5b8d5372076ccde53 *man/smooth.construct.sos.smooth.spec.Rd
 3cb4e59f915c8d64b90754eaeeb5a86f *man/smooth.construct.t2.smooth.spec.Rd
 8672633a1fad8df3cb1f53d7fa883620 *man/smooth.construct.tensor.smooth.spec.Rd
-3e6d88ef6a8ab21bd6f120120602dcf6 *man/smooth.construct.tp.smooth.spec.Rd
+a088e69cf148d07a78ce95a69759c95c *man/smooth.construct.tp.smooth.spec.Rd
 d4083ff900aa69fa07610e0af2a2987b *man/smooth.terms.Rd
 ca676e32c30ddaedbce9f706327cc55b *man/smooth2random.Rd
 844f9653d74441293d05a24dd3e2876a *man/smoothCon.Rd
@@ -153,35 +155,36 @@ a0b0988dba55cca5b4b970e035e3c749 *man/t2.Rd
 a27690f33b9a7bd56d9c1779c64896cc *man/te.Rd
 6eebb6ef90374ee09453d6da6449ed79 *man/tensor.prod.model.matrix.Rd
 f22f1cee0ff2b70628846d1d0f8e9a66 *man/trichol.Rd
-94154ff18af819a7bb83919ee10db0de *man/uniquecombs.Rd
+87e6b4437d00fab4fc814f4cefa3795c *man/trind.generator.Rd
+96c48dd705710f639d76a0d0cc3fb128 *man/uniquecombs.Rd
 a16b3a5a4d13c705dcab8d1cd1b3347e *man/vcov.gam.Rd
 281e73658c726997196727a99a4a1f9e *man/vis.gam.Rd
 07a73758156dfa580c6e92edd34b0654 *man/ziP.Rd
 8bc98d4cb86d851ea0970d68799522cb *man/ziplss.Rd
-92cee4501729e5715f0bb95b498f8425 *po/R-de.po
+f32c7932ef429aa3113007140da82221 *po/R-de.po
 0bdfcf98961b0d52b60f806dc1dba77e *po/R-en@quot.po
-1a3e7a9e3d2f18c8df0ccd49b66bf26d *po/R-fr.po
-6855e193db0dd078d534c913c21a0398 *po/R-ko.po
-461a961615a7c780d6ce5d9694b112d9 *po/R-mgcv.pot
-dcca72cd51cd749f5c54478e31abb706 *po/R-pl.po
+18713570046eecabea9478b9b864a2cf *po/R-fr.po
+d24a793282e6f43c9b67235e46102199 *po/R-ko.po
+753bc70ea9b8cb7c06f658a1eb81ce36 *po/R-mgcv.pot
+1bb61f6df62fc1787864f66933757701 *po/R-pl.po
 8c33d89a914170dbc9f7c5fe598d2135 *po/de.po
 93f72334356fe6f05a64e567efd35c8e *po/en@quot.po
 b3dfaf74ca2d76ca26eec986a14f5584 *po/fr.po
 9116fc041ab458e49b3e498f8c0ac0d9 *po/ko.po
-24f393ff94fa39c8e66250eb0d0e2fcd *po/mgcv.pot
+cd15f571b582fdd31e597083a5fe66ff *po/mgcv.pot
 ed7cb61912e4990cb0076d4cdcf11da8 *po/pl.po
 03972284b3400cf82cacd5d2dc4b8cb3 *src/Makevars
 342aa30c8f6f1e99ffa2576a6f29d7ce *src/coxph.c
-555f6948880bff0b6fa23eeb51609c1c *src/discrete.c
-9721ea07a126278eacd041c814161968 *src/gdi.c
+0d723ffa162b4cb0c2c0fa958ccb4edd *src/discrete.c
+dba99f7d7cc412dd9255f6307c8c7fa7 *src/gdi.c
 2436f9b328e80370ce2203dbf1dd813c *src/general.h
 acf1a2cff05d68856f2b6acee1e63cf7 *src/init.c
-5d9a48e07b438a7c3b32c94fe67ae30c *src/magic.c
-195fa3d343a58a362dfac0ce75b5eed7 *src/mat.c
+a9151b5236852eef9e6947590bfcf88a *src/magic.c
+16b3a72d8bf608026f45de356cf37e65 *src/mat.c
 0545dabf3524a110d616ea5e6373295d *src/matrix.c
 6b781cbd5b9cfee68ad30bb7ce31ef3a *src/matrix.h
 48cde0e19d5dd54b131ba66c777c0ec2 *src/mgcv.c
-3e6fb646794403a7ad4d375f286117f0 *src/mgcv.h
+473bc0f430ec18b53d4649c7276c2962 *src/mgcv.h
 97e3717e95a70b1470b4c3071e144d17 *src/misc.c
 465b8790ca2dfb6e8c5635cacabf5460 *src/mvn.c
 8f480dc455f9ff011c3e8f059efec2c5 *src/qp.c

NAMESPACE

@@ -8,7 +8,9 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
        fs.test,fs.boundary,gam, gam2derivative, 
        gam2objective,
        gamm, gam.check, gam.control,gam.fit3,
-       gam.fit, gam.outer,gam.vcomp, gamSim , 
+       gam.fit, 
+       gamlss.etamu,gamlss.gH,
+       gam.outer,gam.vcomp, gamSim , 
        gaulss,gam.side,get.var,gevlss,
        influence.gam, 
        in.out,inSide,interpret.gam,initial.sp,
@@ -68,7 +70,8 @@ export(anova.gam, bam, bam.update,bandchol, betar, cox.ph,concurvity,
        summary.gam,sp.vcov,
        spasm.construct,spasm.sp,spasm.smooth,
        t2,te,ti,tensor.prod.model.matrix,tensor.prod.penalties,
-       trichol,Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
+       trichol,trind.generator,
+       Tweedie,tw,uniquecombs, vcov.gam, vis.gam, ziP, ziplss)
 
 importFrom(grDevices,cm.colors,dev.interactive,devAskNewPage,gray,grey,heat.colors,terrain.colors,topo.colors)
 importFrom(graphics,abline,axis,axTicks,box,contour,hist,image,lines,

R/efam.r

@@ -439,7 +439,7 @@ ocat <- function(theta=NULL,link="identity",R=NULL) {
   preinitialize <- expression({
     ocat.ini <- function(R,y) {
     ## initialize theta from raw counts in each category
-      if (R<3) return
+      if (R<3) return()
       y <- c(1:R,y) ## make sure there is *something* in each class
       p <- cumsum(tabulate(y[is.finite(y)])/length(y[is.finite(y)]))
       eta <- if (p[1]==0) 5 else -1 - log(p[1]/(1-p[1])) ## mean of latent

R/fast-REML.r

@@ -941,7 +941,13 @@ fast.REML.fit <- function(Sl,X,y,rho,L=NULL,rho.0=NULL,log.phi=0,phi.fixed=TRUE,
   ## idea in following is only to exclude terms with zero first and second derivative
   ## from optimization, as it is only these that slow things down if included...  
   uconv.ind <- (abs(grad) > reml.scale*conv.tol*.1)|(abs(grad2)>reml.scale*conv.tol*.1)
-
+  ## if all appear converged at this stage, then there is probably something wrong,
+  ## but reset anyway to see if situation can be recovered. If we don't do this then
+  ## need to abort immediately, otherwise fails trying to eigen a 0 by 0 matrix
+  if (sum(uconv.ind)==0) {
+    warning("Possible divergence detected in fast.REML.fit",call.=FALSE,immediate.=TRUE)
+    uconv.ind <- rep(TRUE,length(grad)) 
+  }
   step.failed <- FALSE
   for (iter in 1:200) { ## the Newton loop
     ## Work only with unconverged (much quicker under indefiniteness)

R/gam.fit3.r

@@ -663,7 +663,8 @@ gam.fit3 <- function (x, y, sp, Eb,UrS=list(),
             pearson <- sum(weights*(y-mu)^2/family$variance(mu))
             scale.est <- (pearson+dev.extra)/(n.true-trA)
             if (control$scale.est%in%c("fletcher","Fletcher")) { ## Apply Fletcher (2012) correction
-              s.bar = mean(family$dvar(mu)*(y-mu)*sqrt(weights)/family$variance(mu))
+              ## note limited to 10 times Pearson...
+              s.bar = max(-.9,mean(family$dvar(mu)*(y-mu)*sqrt(weights)/family$variance(mu)))
               if (is.finite(s.bar)) scale.est <- scale.est/(1+s.bar)
             }
          } else { ## use the deviance estimator
@@ -1927,25 +1928,30 @@ bfgs <-  function(lsp,X,y,Eb,UrS,L,lsp0,offset,U1,Mp,family,weights,
      ### Derivative testing code. Not usually called and not part of BFGS...
      ok <- check.derivs
      while (ok) { ## derivative testing
-       deriv <- 1
+       #deriv <- 1
        ok <- FALSE ## set to TRUE to re-run (e.g. with different eps)
        deriv.check(x=X, y=y, sp=L%*%lsp+lsp0, Eb=Eb,UrS=UrS,
-         offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=deriv,
+         offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1,
          control=control,gamma=gamma,scale=scale,
          printWarn=FALSE,mustart=mustart,start=start,
          scoreType=scoreType,eps=eps,null.coef=null.coef,Sl=Sl,...)
-
-       fdH <- b$dH
-       fdb.dr <- b$db.drho*0
+       ## deal with fact that deriv might be 0...	 
+       bb <- if (deriv==1) b else gam.fit3(x=X, y=y, sp=L%*%lsp+lsp0,Eb=Eb,UrS=UrS,
+                    offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1,
+                    control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=prev$start,
+                    mustart=prev$mustart,scoreType=scoreType,null.coef=null.coef,
+                    pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
+       fdH <- bb$dH
+       fdb.dr <- bb$db.drho*0
        for (j in 1:length(lsp)) { ## check dH and db.drho
          lsp1 <- lsp;lsp1[j] <- lsp[j] + eps
          ba <- gam.fit3(x=X, y=y, sp=L%*%lsp1+lsp0,Eb=Eb,UrS=UrS,
-                    offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=deriv,
+                    offset = offset,U1=U1,Mp=Mp,family = family,weights=weights,deriv=1,
                     control=control,gamma=gamma,scale=scale,printWarn=FALSE,start=prev$start,
                     mustart=prev$mustart,scoreType=scoreType,null.coef=null.coef,
                     pearson.extra=pearson.extra,dev.extra=dev.extra,n.true=n.true,Sl=Sl,...)
-         fdH[[j]] <- (ba$H - b$H)/eps
-         fdb.dr[,j] <- (ba$coefficients - b$coefficients)/eps
+         fdH[[j]] <- (ba$H - bb$H)/eps
+         fdb.dr[,j] <- (ba$coefficients - bb$coefficients)/eps
        }
      } 
      ### end of derivative testing. BFGS code resumes...

R/gam.fit4.r

@@ -278,23 +278,34 @@ gam.fit4 <- function(x, y, sp, Eb,UrS=list(),
   ## need an initial `null deviance' to test for initial divergence...
   ## if (!is.null(start)) null.coef <- start - can be on edge of feasible - not good
   null.eta <- as.numeric(x%*%null.coef + as.numeric(offset))
-  old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef 
 
-  if (!is.null(start)) { ## check it's at least better than null.coef
-    pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
-    if (pdev>old.pdev) start <- mustart <- etastart <- NULL
-  }
+  #old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef 
+
+  #if (!is.null(start)) { ## check it's at least better than null.coef
+  #  pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
+  #  if (pdev>old.pdev) start <- mustart <- etastart <- NULL
+  #}
 
   ## call the families initialization code...
 
   if (is.null(mustart)) {
     eval(family$initialize)
+    mukeep <- NULL
   } else {
     mukeep <- mustart
     eval(family$initialize)
-    mustart <- mukeep
+    #mustart <- mukeep
   }
-
+
+  old.pdev <- sum(dev.resids(y, linkinv(null.eta), weights,theta)) + t(null.coef)%*%St%*%null.coef 
+
+  if (!is.null(start)) { ## check it's at least better than null.coef
+    pdev <- sum(dev.resids(y, linkinv(x%*%start+as.numeric(offset)), weights,theta)) + t(start)%*%St%*%start
+    if (pdev>old.pdev) start <- mukeep <- etastart <- NULL
+  }
+
+  if (!is.null(mukeep)) mustart <- mukeep
+
   ## and now finalize initialization of mu and eta...
 
   eta <- if (!is.null(etastart)) etastart

R/gamlss.r

@@ -1,4 +1,4 @@
-## (c) Simon N. Wood (2013,2014) distributed under GPL2
+## (c) Simon N. Wood (2013-2016) distributed under GPL2
 ## Code for the gamlss families.
 ## idea is that there are standard functions converting
 ## derivatives w.r.t. mu to derivatives w.r.t. eta, given 
@@ -186,8 +186,8 @@ gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=
 ## 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
+##        1 - diagonal of first deriv of Hess wrt sps
+##        2 - first deriv of Hess wrt sps
 ##        3 - everything.
   K <- length(jj)
   p <- ncol(X);n <- nrow(X)
@@ -304,7 +304,11 @@ gamlss.gH0 <- function(X,jj,l1,l2,i2,l3=0,i3=0,l4=0,i4=0,d1b=0,d2b=0,deriv=0,fh=
     }
   } ## if deriv>2
 
-  list(lb=lb,lbb=lbb,d1H=d1H,d2H=d2H,trHid2H=trHid2H)
+  list(lb=lb, ## grad w.r.t. coefs
+       lbb=lbb, ## hess w.r.t. coefs, H
+       d1H=d1H, ## grad of H wrt sps 
+       ## d2H=d2H,
+       trHid2H=trHid2H) ## tr(H^{-1}d^2H/dspsp)
 } ## end of gamlss.gH0
 
 

R/mgcv.r

@@ -44,9 +44,19 @@ rig <- function(n,mean,scale) {
 ## inverse guassian deviates generated by algorithm 5.7 of
 ## Gentle, 2003. scale = 1/lambda. 
   if (length(n)>1) n <- length(n)
-  y <- rnorm(n)^2
-  mu2 <- 0*y + mean^2 ## y is there to ensure mu2 is a vector
-  x <- mean + 0.5*scale*(mu2*y - mean*sqrt(4*mean*y/scale + mu2*y^2))
+  x <- y <- rnorm(n)^2
+  mys <- mean*scale*y
+  mu <- 0*y + mean ## y is there to ensure mu is a vector
+  mu2 <- mu^2;
+  ind <- mys < .Machine$double.eps^-.5 ## cut off for tail computation
+  x[ind] <- mu[ind]*(1 + 0.5*(mys[ind] - sqrt(mys[ind]*4+mys[ind]^2)))
+  x[!ind] <- mu[!ind]/mys[!ind] ## tail term (derived from Taylor of sqrt(1+eps) etc)
+  #my <- mean*y; sc <- 0*y + scale
+  #ind <- my > 1 ## cancellation error can be severe, without splitting
+  #x[!ind] <- mu[!ind]*(1 + 0.5*sc[!ind]*(my[!ind] - sqrt(4*my[!ind]/sc[!ind] + my[!ind]^2)))
+  ## really the sqrt in the next term should be expanded beyond first order and then
+  ## worked on - otherwise too many exact zeros?
+  #x[ind] <- pmax(0,mu[ind]*(1+my[ind]*.5*sc[ind]*(1-sqrt(1+ 4/(sc[ind]*my[ind])))))
   ind <- runif(n) > mean/(mean+x)
   x[ind] <- mu2[ind]/x[ind]
   x ## E(x) = mean; var(x) = scale*mean^3
@@ -1580,7 +1590,7 @@ estimate.gam <- function (G,method,optimizer,control,in.out,scale,gamma,start=NU
        #}
        ## make sure optimizer appropriate for available derivatives
        if (!is.null(G$family$available.derivs)) {
-         if (G$family$available.derivs==1) optimizer <- c("outer","bfgs")
+         if (G$family$available.deriv==1 && optimizer[1]!="efs")  optimizer <- c("outer","bfgs")
 	 if (G$family$available.derivs==0) optimizer <- "efs"
        }	 
     }
@@ -3863,7 +3873,7 @@ cooks.distance.gam <- function(model,...)
 sp.vcov <- function(x) {
 ## get cov matrix of smoothing parameters, if available
   if (!inherits(x,"gam")) stop("argument is not a gam object")
-  if (x$method%in%c("ML","P-ML","REML","P-REML")&&!is.null(x$outer.info$hess)) {
+  if (x$method%in%c("ML","P-ML","REML","P-REML","fREML")&&!is.null(x$outer.info$hess)) {
     return(solve(x$outer.info$hess))
   } else return(NULL)
 }
@@ -3873,7 +3883,7 @@ gam.vcomp <- function(x,rescale=TRUE,conf.lev=.95) {
 ## in a fitted `gam' object.
   if (!inherits(x,"gam")) stop("requires an object of class gam")
   if (!is.null(x$reml.scale)&&is.finite(x$reml.scale)) scale <- x$reml.scale else scale <- x$sig2
-  if (length(x$sp)==0) return
+  if (length(x$sp)==0) return()
   if (rescale) { ## undo any rescaling of S[[i]] that may have been done
     m <- length(x$smooth)
     if (is.null(x$paraPen)) {