[r-cran-gam] 04/20: Imported Upstream version 1.08
Andreas Tille
tille at debian.org
Fri Dec 16 14:32:10 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-gam.
commit a622f76aa2e07bfd8277a15ee2c4b1c800d05a20
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 13:32:21 2016 +0100
Imported Upstream version 1.08
---
DESCRIPTION | 15 +--
MD5 | 30 +++---
R/gam.exact.R | 2 +-
R/gam.wlist.R | 1 +
R/{onLoad.R => onAttach.R} | 2 +-
R/predict.gam.R | 8 +-
inst/ratfor/backfit.r | 235 ++++++++++++++++++++++++++++++++++++++++++++
inst/ratfor/backlo.r | 170 ++++++++++++++++++++++++++++++++
inst/ratfor/linear.r | 78 +++++++--------
inst/ratfor/lo.r | 187 +++++++++++++++++++++++++++++++++++
inst/ratfor/splsm.r | 239 +++++++++++++++++++++++++++++++++++++++++++++
man/predict.gam.Rd | 2 +-
src/Makevars.win | 2 +-
src/backfit.f | 10 +-
src/backlo.f | 12 +--
src/linear.f | 82 ++++++++--------
src/lo.f | 12 +--
src/splsm.f | 14 +--
18 files changed, 970 insertions(+), 131 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 2035559..805b7f4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,17 +1,18 @@
Package: gam
Type: Package
Title: Generalized Additive Models
-Date: 2011-12-05
-Version: 1.06.2
+Date: 2013-04-23
+Version: 1.08
Author: Trevor Hastie
Description: Functions for fitting and working with generalized
- additive models, as described in chapter 7 of "Statistical Models in
- S" (Chambers and Hastie (eds), 1991), and "Generalized Additive
- Models" (Hastie and Tibshirani, 1990).
+ additive models, as described in chapter 7 of "Statistical
+ Models in S" (Chambers and Hastie (eds), 1991), and
+ "Generalized Additive Models" (Hastie and Tibshirani, 1990).
Maintainer: Trevor Hastie <hastie at stanford.edu>
Depends: stats, splines
Suggests: akima
License: GPL-2
-Packaged: 2011-12-06 01:44:56 UTC; hastie
+Packaged: 2013-04-24 02:02:27 UTC; hastie
+NeedsCompilation: yes
Repository: CRAN
-Date/Publication: 2011-12-06 07:00:29
+Date/Publication: 2013-04-24 09:46:10
diff --git a/MD5 b/MD5
index 6b99772..6610c5a 100644
--- a/MD5
+++ b/MD5
@@ -1,4 +1,4 @@
-88414eb45ea6971765a99732c9f3ba98 *DESCRIPTION
+0d66a823f3136a9c9fbdbdff74965935 *DESCRIPTION
af77f82fb0aa5e383808c5f36aa47066 *INDEX
105cd87b9c243fbe33aabf84cf2b3271 *NAMESPACE
3f56ee7eddd13ec792f7c2150b1e1eca *R/all.wam.R
@@ -12,7 +12,7 @@ c67ddb150e807c4a0ef77acf132c1021 *R/as.anova.R
4ab8260f7666ffd65e9ef915fb245531 *R/deviance.lm.R
60657ce607548e744477b5b157776116 *R/gam.R
6f9d6d8a11d7b20d791233f1ed3445dc *R/gam.control.R
-4be1a4ec9218947bab2f6f708acbb03d *R/gam.exact.R
+2c633d79441282175358f1a0c148bd50 *R/gam.exact.R
52ac1d379ebc2189d982e3f86e62c8af *R/gam.fit.R
bd96ac15ba688c863f09f0e632a411f2 *R/gam.lo.R
946d0c85f06753f768700803d197c247 *R/gam.match.R
@@ -22,7 +22,7 @@ bd96ac15ba688c863f09f0e632a411f2 *R/gam.lo.R
a6bc1a9e490a60aaa1d5d8d47b872c00 *R/gam.scope.R
6ccfc23d8ade2bae10532a2f0d6ad683 *R/gam.slist.R
1f8d4b4e6776250d0e4080e64a1273e9 *R/gam.sp.R
-c83f70ebeee69ffd30cb96bcabd2c2c7 *R/gam.wlist.R
+865f790f7ee2bcf180a2741fba26f7ac *R/gam.wlist.R
3145f500d12b60398d2bf36eb20973c9 *R/gamlist.R
6bc9b975aa99176b8d021818760af91c *R/gplot.R
f8026ded9300952fe6a9810e8bfeecdf *R/gplot.default.R
@@ -35,11 +35,11 @@ f8026ded9300952fe6a9810e8bfeecdf *R/gplot.default.R
f8e577af7d28bdf5252ef489f4cb824d *R/lo.wam.R
901999a84e9e5db5552b706a99fcc303 *R/na.gam.replace.R
a82abbfbf9d226c04290dd721e1d8b51 *R/newdata.predict.gam.R
-0b14c3fb8b9f3899d54b63cf090491ec *R/onLoad.R
+39cfb671132d4928c0a77a79edb1967d *R/onAttach.R
9229d04542aa65b686e6a44fc5e5eb52 *R/plot.gam.R
07745ab717cea5503a96ccfbd3d66529 *R/plot.preplot.gam.R
b0b4e3ea29649f25919a9b3b9a01886f *R/polylo.R
-8e6ec4cd97cfd239ae7ef53d3e39be85 *R/predict.gam.R
+23706ce18cdf6313374bbf40e1161a48 *R/predict.gam.R
4314caf8a4bc9da21bd847ea1322146c *R/preplot.gam.R
87477db085b0209bab1a5fc82cca43be *R/print.gam.R
ca3a618d4376069f1d43a2aabaea5e4c *R/print.gamex.R
@@ -54,7 +54,11 @@ addd48040c8374e5a9ed5d587cfe77e7 *R/step.gam.R
ba66638e3de17b868b4d98dffe95009d *data/gam.data.RData
83529cbff37939aff8d96d32d6458f12 *data/gam.newdata.RData
4df6bee1aa3bee402b8f1083dda04b04 *data/kyphosis.RData
-b1c5751278caefdae9c643f6dd3daded *inst/ratfor/linear.r
+a42dba5f95d8760e06a78389d780b170 *inst/ratfor/backfit.r
+a009bd4d2232ec7cc19b9d7d10280ffb *inst/ratfor/backlo.r
+1ca924fd62063613d16c1cf607abb6b2 *inst/ratfor/linear.r
+4e0b184dc647e3abac7fb7f023ed4a69 *inst/ratfor/lo.r
+90ae42a7500cb1b749601e2d4f5656cc *inst/ratfor/splsm.r
ddba04cdcf5d2044919c08bcf71c7786 *man/anova.gam.Rd
d39a839642a1ee2b71d93df9a41c62c5 *man/gam-internal.Rd
ffbb793ebf7414ea69df18b384e40534 *man/gam.Rd
@@ -65,18 +69,18 @@ dd3553bd8578858873bd1384b000273a *man/kyphosis.Rd
800542a98f81e40c0930605f288c9ca4 *man/lo.Rd
5068084693ec99de54e0531a53d4e637 *man/na.gam.replace.Rd
a79905d34d4d93cb9939329b7cc8f507 *man/plot.gam.Rd
-fe97676fb4413c259621109a5d78edac *man/predict.gam.Rd
+7de080b0a4e684b2a2b4e26e17fe38d8 *man/predict.gam.Rd
7f2b4a226c564121816e2594f7fd61c3 *man/s.Rd
03718b63cb73f8be1d573ecca2c9a48f *man/step.gam.Rd
f80de2889856cba0512e5377926af1aa *src/Makevars
-2fa4c7011c2bc0f7449ae151d5cc44ae *src/Makevars.win
-65a8e52762837be96db5d33cf8f859a7 *src/backfit.f
-819e82fd2b11f2acdda88efbf0a9fba0 *src/backlo.f
+1c81cbb8cbdd19a2e909dc5f802ddb07 *src/Makevars.win
+45061d1eb26bda41c4c458126dd303da *src/backfit.f
+7205467d48f1b2f05d6b830aa34afac1 *src/backlo.f
10f91c907532cc70e76074ca47d14716 *src/bsplvd.f
2cf888d1f6a6ca37a7211c7c6ab86c0e *src/bvalue.f
40c1f6b57d6c3a03a1d972f009c5e5df *src/bvalus.f
-f8d84e4eab8700aa419a7aefd1a8ddd4 *src/linear.f
-85c9d000ac0b8cc1dc61a8bfa40f3ed2 *src/lo.f
+114d8313198201d0f01482982e5e8de7 *src/linear.f
+6b4fde12653de078154a99510fed6d64 *src/lo.f
388f5201bdc80e4c25a575f4810915d6 *src/loessc.c
425c0bc459b1dee59d08e5cd6c28d89c *src/loessf.f
6ada085e39fd48968fcd79368c0ccddd *src/modreg.h
@@ -84,6 +88,6 @@ e1c0a8ef61f04b6239ff2ea6874be92d *src/qsbart.f
876032799d52ef84fb846af58939a318 *src/sbart.c
82da999d24034505301e31c78d1e58cc *src/sgram.f
d277bb97eef775673f5fa2da911d81de *src/sinerp.f
-362186247368f8550e59c0bf1cde8eb5 *src/splsm.f
+65db5d599cdd5a5d81110cff681d374a *src/splsm.f
39a0d0a34a95130be92b075a49feac36 *src/sslvrg.f
4f6275039d4731d4f2920fc3de5f61e7 *src/stxwx.f
diff --git a/R/gam.exact.R b/R/gam.exact.R
index df2d231..d0fa89d 100644
--- a/R/gam.exact.R
+++ b/R/gam.exact.R
@@ -111,7 +111,7 @@ function(gam.obj)
this.var <- paste("x",k,sep="")
upd.form <- update(as.formula(form),paste(this.var,"~. -",this.var))
- XX[,k] <- gam(formula=upd.form,data=mydat,family=gaussian,weight=w,
+ XX[,k] <- gam(formula=upd.form,data=mydat,family=gaussian,weights=w,
control=eval(Control))$fitted
}
}
diff --git a/R/gam.wlist.R b/R/gam.wlist.R
index 1a71a07..13a6f30 100644
--- a/R/gam.wlist.R
+++ b/R/gam.wlist.R
@@ -1,2 +1,3 @@
"gam.wlist" <-
c("s","lo")
+
diff --git a/R/onLoad.R b/R/onAttach.R
similarity index 73%
rename from R/onLoad.R
rename to R/onAttach.R
index f637676..64584b2 100644
--- a/R/onLoad.R
+++ b/R/onAttach.R
@@ -1,3 +1,3 @@
-.onLoad=function(libname,pkgname){
+.onAttach=function(libname,pkgname){
packageStartupMessage("Loaded gam ", as.character(packageDescription("gam")[["Version"]]),"\n")
}
diff --git a/R/predict.gam.R b/R/predict.gam.R
index 6f8c764..f335084 100644
--- a/R/predict.gam.R
+++ b/R/predict.gam.R
@@ -30,12 +30,13 @@
terms = {
out <- NextMethod("predict")
TT <- dimnames(s <- object$smooth)[[2]]
+ TT=intersect(terms,TT)##added to protect subsets
out$fit[, TT] <- out$fit[,
- TT] + s
+ TT] + s[,TT]
TS <- out$residual.scale^2
out$se.fit[, TT] <- sqrt(out$
se.fit[, TT]^2 + TS *
- object$var)
+ object$var[,TT])
out
}
)
@@ -43,7 +44,8 @@
terms = {
out <- NextMethod("predict")
TT <- dimnames(s <- object$smooth)[[2]]
- out[, TT] <- out[, TT] + s
+ TT=intersect(terms,TT)##added to protect subsets
+ out[, TT] <- out[, TT] + s[,TT]
out
}
,
diff --git a/inst/ratfor/backfit.r b/inst/ratfor/backfit.r
new file mode 100644
index 0000000..d6a0e91
--- /dev/null
+++ b/inst/ratfor/backfit.r
@@ -0,0 +1,235 @@
+subroutine bakfit(x,npetc,y,w,which,spar,dof,match,nef,
+ etal,s,eta,beta,var,tol,
+ qr,qraux,qpivot,effect,work)
+#integer npetc(7)
+#1:n
+#2:p
+#3:q
+#4:ifvar
+#5:nit
+#6:maxit
+#7:qrank
+#subroutine bakfit(x,n,p,y,w,q,which,spar,dof,match,nef,
+# etal,s,eta,beta,var,ifvar,tol,nit,maxit,
+# qr,qraux,qrank,qpivot,work)
+#This subroutine fits an additive spline fit to y
+#All arguments are either double precision or integer
+# bakfit uses the modified backfitting algorithm described in Buja, Hastie
+# and Tibshirani, Annals of Statistics, 1989. It calls splsm, and some
+# linpack based routines
+# This was written by Trevor Hastie in 1990
+# It has been modified from the S3 version by Trevor Hastie
+# in March 2005, to accommodate the modified sbart routine in R
+# Note that spar has changed, and we change it here to conform with
+# the smooth.spline routine in R
+#INPUT
+#
+#x double dim n by p ; x variables, includes constant
+#n integer number of rows in x
+#p integer number of columns of x
+#y double length n ; y variable for smoothing
+#w double length n ; prior weights for smoothing, > 0
+#q integer number of nonlinear terms
+#which integer length q indices of columns of x for nonlinear fits
+#spar double length q spars for smoothing; see below
+#dof double length q dof for/from smoothing; see below
+#match integer n by q matrix of match'es; see below
+#nef integer q vector of nef's; see below
+#s double n by q nonlinear part of the smooth functions
+# used as starting values. the linear part is
+# irrelevant
+#ifvar logical should the variance information be computed
+#tol double tolerance for backfitting convergence; 0.0005 is good
+#maxit integer maximum number of iterations; 15 is good
+#qr double n by p weighted qr decomposition of x
+#qraux double p belongs with qr
+#qrank integer rank of x ; if qrank=0, then bakfit computes qr and qraux
+#qpivot integer p the columns of qr are rearranged according to pivot
+#effec double n effect vector
+#work double
+# Let nk=max(nef)+2, then
+# work should be (10+2*4)*nk+5*nef+5*n+15 +q double
+#BELOW
+#the following comments come from documentation for splsm
+# they apply to each element of spar,dof match etc
+#spar double smoothing parameter -1.5 <spar<1.5; default is 1
+#dof double equivalent degrees of freedom
+# if dof is 0, spar is used
+# if 0< dof <1, dof = 1
+# if dof >=1, dof is used
+# note: dof does not use the constant term
+#match integer length n -- in S language x[i] == sort(unique(x)[match[i]]
+# match is produced by subroutine namat
+#nef number of unique elements in x; so match has values between 1 and nef+1
+# missing data are given the match number nef+1
+#work double workspace of length (10+2*4)*(nef+2)+5*nef+n+15
+#
+#OUTPUT
+#
+#x,y,w,n,p,which,q,maxit,match,nef are untouched
+#spar for each element of spar:
+# if spar was 0 and dof was 0, then spar is that spar
+# that minimized gcv
+# if spar was 0 and dof > 0, then spar is that which achieves dof
+#dof the dof of the fitted smooth. Note: even if dof was given
+# as 4, it will be returned as say 3.995 which is what
+# spar produces
+#etal double length n linear component of the fit
+#s double n by q nonlinear part of the smooth functions
+#eta double length n fitted values
+#beta double length p linear coefficients
+# So, the centered fitted functions are:
+# b(j)*(x(i,j)-mean(x(.,j)) +s(i,j)
+# where j is an element of which
+#var double n by q
+# if ifvar was .true.
+# the unscaled variance elements for the NONLINEAR
+# and UNIQUE part of s, in the order of sort(unique(x))
+# var is lev(i)/w(i) -h(i)/w where h(i) is the hat element from
+# the simple weighted least squares fit. This is used in gamcov
+#
+#nit number of iterations used
+#qr etc the qr is returned
+
+implicit double precision(a-h,o-z)
+logical ifvar
+integer npetc(7),iter
+integer n,p,q,which(*),match(*),nef(*),nit,maxit,qrank,qpivot(*)
+double precision x(*),y(*),w(*),spar(*),dof(*),
+ etal(*),s(*),eta(*),beta(*),var(*),tol,
+ qr(*),qraux(*),effect(*),work(*)
+n=npetc(1)
+p=npetc(2)
+q=npetc(3)
+ifvar=.false.
+if(npetc(4)==1)ifvar=.true.
+maxit=npetc(6)
+qrank=npetc(7)
+do i=1,q{work(i)=dof(i)}
+call backf1(x,n,p,y,w,q,which,spar,dof,match,nef,
+ etal,s,eta,beta,var,ifvar,tol,nit,maxit,
+ qr,qraux,qrank,qpivot,effect,work(q+1),work(q+n+1),
+ work(q+2*n+1),work(q+3*n+1),work(q+4*n+1))
+npetc(7)=qrank
+return
+end
+
+subroutine backf1(x,n,p,y,w,q,which,spar,dof,match,nef,
+ etal,s,eta,beta,var,ifvar,tol,nit,maxit,
+ qr,qraux,qrank,qpivot,effect,z,old,sqwt,sqwti,work)
+implicit double precision(a-h,o-z)
+logical ifvar
+integer n,p,q,which(q),match(n,q),nef(q),nit,maxit,qrank,qpivot(p)
+double precision x(n,p),y(n),w(n),spar(q),dof(q),
+ etal(n),s(n,q),eta(n),beta(p),var(n,q),tol,
+ qr(n,p),qraux(p),effect(n),work(*)
+double precision z(*),old(*),dwrss,ratio
+double precision sqwt(n),sqwti(n)
+logical anyzwt
+double precision deltaf, normf,onedm7
+integer job,info
+onedm7=1d-7
+job=1101;info=1
+if(q==0)maxit=1
+ratio=1d0
+# fix up sqy's for weighted problems.
+anyzwt=.false.
+do i=1,n{
+ if(w(i)>0d0){
+ sqwt(i)=dsqrt(w(i))
+ sqwti(i)=1d0/sqwt(i)
+ }
+ else{
+ sqwt(i)=0d0
+ sqwti(i)=0d0
+ anyzwt=.true.
+ }
+ }
+# if qrank > 0 then qr etc contain the qr decomposition
+# else bakfit computes it.
+if(qrank==0){
+ do i=1,n{
+ do j=1,p{
+ qr(i,j)=x(i,j)*sqwt(i)
+ }
+ }
+ do j=1,p{qpivot(j)=j}
+ call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7)
+ }
+do i=1,n{
+ eta(i)=0d0
+ for(j=1;j<=q;j=j+1){
+ eta(i)=eta(i)+s(i,j)
+ }
+ }
+nit=0
+while ((ratio > tol )&(nit < maxit)){
+ # first the linear fit
+ deltaf=0d0
+ nit=nit+1
+ do i=1,n{
+ z(i)=(y(i)-eta(i))*sqwt(i)
+ old(i)=etal(i)
+ }
+# call dqrsl1(qr,dq,qraux,qrank,sqz,one,work(1),etal,two,three)
+#job=1101 -- computes fits, effects and beta
+ call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta,
+ work(1),etal,job,info)
+
+# now unsqrt the fits
+#Note: we dont have to fix up the zero weights till the end, since their fits
+#are always immaterial to the computation
+ do i=1,n{
+ etal(i)=etal(i)*sqwti(i)
+ }
+ # now a single non-linear backfitting loop
+ for(k=1;k<=q;k=k+1){
+ j=which(k)
+ do i=1,n{
+ old(i)=s(i,k)
+ z(i)=y(i)-etal(i)-eta(i)+old(i)
+ }
+ # this uses spar to set smoothing after iteration 1
+ if(nit>1){dof(k)=0d0}
+ call splsm(x(1,j),z,w,n,match(1,k),nef(k),spar(k),
+ dof(k),s(1,k),s0,var(1,k),ifvar,work)
+ do i=1,n{
+ eta(i)=eta(i)+s(i,k)-old(i)
+ etal(i)=etal(i)+s0
+ }
+ deltaf=deltaf+dwrss(n,old,s(1,k),w)
+ }
+ normf=0d0
+ do i=1,n{
+ normf=normf+w(i)*eta(i)*eta(i)
+ }
+ if(normf>0d0){
+ ratio=dsqrt(deltaf/normf)
+ }
+ else {ratio = 0d0}
+# call DBLEPR("ratio",-1,ratio,1)
+ }
+#now package up the results
+do j=1,p {work(j)=beta(j)}
+do j=1,p {beta(qpivot(j))=work(j)}
+if(anyzwt){
+ do i=1,n {
+ if(w(i) <= 0d0){
+ etal(i)=0d0
+ do j=1,p{
+ etal(i)=etal(i)+beta(j)*x(i,j)
+ }
+ }
+ }
+ }
+
+do i=1,n
+ eta(i)=eta(i)+etal(i)
+
+do j=1,q {
+ call unpck(n,nef(j),match(1,j),var(1,j),old)
+ do i=1,n {var(i,j)=old(i)}
+ }
+
+return
+end
diff --git a/inst/ratfor/backlo.r b/inst/ratfor/backlo.r
new file mode 100644
index 0000000..2df87f2
--- /dev/null
+++ b/inst/ratfor/backlo.r
@@ -0,0 +1,170 @@
+subroutine baklo(x,y,w,npetc,wddnfl,spatol,match,
+ etal,s,eta,beta,var,dof,
+ qr,qraux,qpivot,effect,iv,v,iwork,work)
+implicit double precision(a-h,o-z)
+integer n,p,q,nit,maxit,qrank
+integer npetc(7),wddnfl(*),match(*),qpivot(*),iv(*),iwork(*)
+##integer which(q),dwhich(q),degree(q),nef(q),liv(q),lv(q),nvmax(q)
+double precision x(*),y(*),w(*),spatol(*),
+ etal(*),s(*),eta(*),beta(*),var(*),dof(*),
+ qr(*),qraux(*),v(*),effect(*),work(*)
+#work size: 4*n + sum( nef(k)*(pj+dj+4)+5+3*pj ) +5*n
+# = 9*n + sum( nef(k)*(pj+dj+4)+5+3*pj )
+
+n=npetc(1)
+p=npetc(2)
+q=npetc(3)
+maxit=npetc(5)
+qrank=npetc(6)
+call baklo0(x,n,p,y,w,q,wddnfl(1),wddnfl(q+1),wddnfl(2*q+1),
+ spatol(1),wddnfl(3*q+1),dof,match,wddnfl(4*q+1),
+ etal,s,eta,beta,var,spatol(q+1),
+ nit,maxit,qr,qraux,qrank,qpivot,effect,
+ work(1),work(n+1),work(2*n+1),work(3*n+1),
+ iv,wddnfl(5*q+1),wddnfl(6*q+1),v,wddnfl(7*q+1),
+ iwork(1),work(4*n+1))
+npetc(4)=nit
+npetc(6)=qrank
+return
+end
+
+subroutine baklo0(x,n,p,y,w,q,which,dwhich,pwhich,span,degree,dof,match,nef,
+ etal,s,eta,beta,var,tol,nit,maxit,
+ qr,qraux,qrank,qpivot,effect,z,old,sqwt,sqwti,
+ iv,liv,lv,v,nvmax,iwork,work)
+implicit double precision(a-h,o-z)
+integer n,p,q,which(q),dwhich(q),pwhich(q),degree(q),match(n,q),nef(q),nit,
+ maxit,qrank,qpivot(p),iv(*),liv(q),lv(q),nvmax(q),iwork(q)
+double precision x(n,p),y(n),w(n),span(q),dof(q),
+ etal(n),s(n,q),eta(n),beta(p),var(n,q),tol,
+ qr(n,p),qraux(p),v(*),effect(n),work(*)
+#work should be sum( nef(k)*(pj+dj+4)+5+3*pj ) +5*n
+double precision z(*),old(*),dwrss,ratio
+double precision sqwt(n),sqwti(n)
+logical anyzwt
+double precision deltaf, normf,onedm7
+integer job,info,slv,sliv,iw,j,dj,pj
+onedm7=1d-7
+job=1101;info=1
+if(q==0)maxit=1
+ratio=1d0
+# fix up sqy's for weighted problems.
+anyzwt=.false.
+do i=1,n{
+ if(w(i)>0d0){
+ sqwt(i)=dsqrt(w(i))
+ sqwti(i)=1d0/sqwt(i)
+ }
+ else{
+ sqwt(i)=0d0
+ sqwti(i)=0d0
+ anyzwt=.true.
+ }
+ }
+# if qrank > 0 then qr etc contain the qr decomposition
+# else baklo computes it.
+if(qrank==0){
+ do i=1,n{
+ do j=1,p{
+ qr(i,j)=x(i,j)*sqwt(i)
+ }
+ }
+ do j=1,p{qpivot(j)=j}
+ call dqrdca(qr,n,n,p,qraux,qpivot,work,qrank,onedm7)
+ }
+do i=1,n{
+ eta(i)=0d0
+ for(j=1;j<=q;j=j+1){
+ eta(i)=eta(i)+s(i,j)
+ }
+ }
+nit=0
+while ((ratio > tol )&(nit < maxit)){
+ # first the linear fit
+ deltaf=0d0
+ nit=nit+1
+ do i=1,n{
+ z(i)=(y(i)-eta(i))*sqwt(i)
+ old(i)=etal(i)
+ }
+# call dqrsl1(qr,dq,qraux,qrank,sqz,one,work(1),etal,two,three)
+#job=1101 -- computes fits, effects and beta
+ call dqrsl(qr,n,n,qrank,qraux,z,work(1),effect(1),beta,
+ work(1),etal,job,info)
+
+# now unsqrt the fits
+#Note: we dont have to fix up the zero weights till the end, since their fits
+#are always immaterial to the computation
+ do i=1,n{
+ etal(i)=etal(i)*sqwti(i)
+ }
+ # now a single non-linear backfitting loop
+ sliv=1
+ slv=1
+ iw=5*n+1
+ for(k=1;k<=q;k=k+1){
+
+ j=which(k)
+ dj=dwhich(k)
+ pj=pwhich(k)
+ do i=1,n{
+ old(i)=s(i,k)
+ z(i)=y(i)-etal(i)-eta(i)+old(i)
+ }
+call lo1(x(1,j),z,w,n,dj,pj,nvmax(k),span(k),degree(k),match(1,k),
+ nef(k),nit,dof(k),s(1,k),var(1,k),work(iw),
+# xin,win
+ work(iw+pj+1),work(iw+nef(k)*dj+pj+1),
+# sqwin,sqwini,
+ work(iw+nef(k)*(dj+1)+pj+2),work(iw + nef(k)*(dj+2)+pj+2),
+# xqr,qrank,
+ work(iw+nef(k)*(dj+3)+pj+2),work(iw+nef(k)*(pj+dj+4)+pj+2),
+# qpivot,qraux,
+# work(iw+nef(k)*(pj+dj+4)+pj+3),work(iw+nef(k)*(pj+dj+4)+4+2*pj),
+ iwork(1),work(iw+nef(k)*(pj+dj+4)+4+2*pj),
+ iv(sliv),liv(k),lv(k),v(slv),
+ work(1) )
+#work should be sum( nef(k)*(pj+dj+4)+5+3*pj ) +5*n
+# In the call above I give lo1 pieces of work to use for storing
+# the qr decomposition, and it gets the same undisturbed portion
+# each time. The fact that it is given a double work word for qrank
+# is irrelevant but convenient; it still stores the integer qrank there.
+# I do this because there is a partition like this for each lo() term
+# in the model, and the number of them is variable
+ sliv=sliv+liv(k)
+ slv=slv+lv(k)
+ iw=iw+nef(k)*(pj+dj+4)+5+3*pj
+ do i=1,n{
+ eta(i)=eta(i)+s(i,k)-old(i)
+ }
+ deltaf=deltaf+dwrss(n,old,s(1,k),w)
+ }
+ normf=0d0
+ do i=1,n{
+ normf=normf+w(i)*eta(i)*eta(i)
+ }
+ if(normf>0d0){
+ ratio=dsqrt(deltaf/normf)
+ }
+ else {ratio = 0d0}
+ }
+#now package up the results
+do j=1,p {work(j)=beta(j)}
+do j=1,p {beta(qpivot(j))=work(j)}
+if(anyzwt){
+ do i=1,n {
+ if(w(i) <= 0d0){
+ etal(i)=0d0
+ do j=1,p{
+ etal(i)=etal(i)+beta(j)*x(i,j)
+ }
+ }
+ }
+ }
+
+do i=1,n
+ eta(i)=eta(i)+etal(i)
+
+
+return
+end
diff --git a/inst/ratfor/linear.r b/inst/ratfor/linear.r
index fa755a9..de4a18f 100644
--- a/inst/ratfor/linear.r
+++ b/inst/ratfor/linear.r
@@ -1,7 +1,7 @@
subroutine dqrls(x,dx,pivot,qraux,y,dy,beta,res,qt,tol,scrtch,rank)
-integer pivot(1),dx(2),dy(2),rank
-double precision x(1), qraux(1), y(1), beta(1),res(1),qt(1),tol(1),
- scrtch(1)
+integer pivot(*),dx(2),dy(2),rank
+double precision x(*), qraux(*), y(*), beta(*),res(*),qt(*),tol(*),
+ scrtch(*)
integer n,p,q,kn,kp,k,info
@@ -19,7 +19,7 @@ end
#apply the qr decomposition to do various jobs
subroutine dqrsl1(qr,dq,qra,rank,y,k,qy,qb,job,info)
-double precision qr(1),qra(1),y(1),qy(1),qb(1); integer dq(2),job,k,rank
+double precision qr(*),qra(*),y(*),qy(*),qb(*); integer dq(2),job,k,rank
integer n,kn,kb,j
double precision ourqty(1), ourqy(1), ourb(1), ourrsd(1), ourxb(1)
ourqty(1) = 0d0
@@ -62,8 +62,8 @@ return
end
subroutine dqr(x,dx,pivot,qraux,tol,scrtch,rank)
-integer pivot(1),dx(2),rank
-double precision x(1), qraux(1), tol(1), scrtch(1)
+integer pivot(*),dx(2),rank
+double precision x(*), qraux(*), tol(*), scrtch(*)
integer n,p
@@ -77,8 +77,8 @@ end
# ordering and rank estimation
subroutine dqrdca(x,ldx,n,p,qraux,jpvt,work,rank,eps)
integer ldx,n,p,rank
-integer jpvt(1)
-double precision x(ldx,1),qraux(1),work(1),eps
+integer jpvt(*)
+double precision x(ldx,*),qraux(*),work(*),eps
integer j,jj,jp,l,lup,curpvt
double precision dnrm2,tt
double precision ddot,nrmxl,t,ww
@@ -325,8 +325,8 @@ return
end
subroutine chol(a,p,work,jpvt,job,info)
-integer p,jpvt(1),job,info(1)
-double precision a(p,1),work(1)
+integer p,jpvt(*),job,info(*)
+double precision a(p,*),work(*)
integer i,j
for(j =2; j<=p; j = j+1)
for(i=1; i<j; i = i+1)
@@ -340,7 +340,7 @@ end
#x is a real symmetric matrix
subroutine crs(x,dmx,matz,w,z,fv1,fv2,ierr)
-double precision x(1),w(1),z(1),fv1(1),fv2(1)
+double precision x(*),w(*),z(*),fv1(*),fv2(*)
integer dmx(2),nx,nv,ierr,matz
nx=dmx(1)
nv=dmx(2)
@@ -349,9 +349,9 @@ return
end
subroutine dqrls2(x,dx,pivot,qraux,y,dy,beta,res,qt,scrtch,eps)
-integer pivot(1),dx(2),dy(2)
-double precision x(1), qraux(1), y(1), beta(1),res(1),qt(1),
- scrtch(1),eps
+integer pivot(*),dx(2),dy(2)
+double precision x(*), qraux(*), y(*), beta(*),res(*),qt(*),
+ scrtch(*),eps
integer n,p,q,kn,kp,k,info,rank
@@ -368,7 +368,7 @@ return
end
subroutine dsvdc1(x,dmx,job,work,e,s,u,v,info)
-double precision x(1),work(1),s(1),e(1),u(1),v(1)
+double precision x(*),work(*),s(*),e(*),u(*),v(*)
integer dmx(2),nx,nv,job,info
nx=dmx(1)
nv=dmx(2)
@@ -1617,7 +1617,7 @@ end
subroutine dmatp(x,dx,y,dy,z)
integer dx(2),dy(2)
-double precision x(1), y(1),z(1),ddot
+double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j
@@ -1634,7 +1634,7 @@ end
subroutine dmatpt(x,dx,y,dy,z)
integer dx(2),dy(2)
-double precision x(1), y(1),z(1),ddot
+double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j,ii
@@ -1652,9 +1652,9 @@ end
subroutine matpm(x,dx,mmx,mx,y,dy,mmy,my,z)
integer dx(2),dy(2)
-integer mmx(1), mmy(1)
-integer mx(1), my(1)
-double precision x(1), y(1),z(1),ddot
+integer mmx(*), mmy(*)
+integer mx(*), my(*)
+double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j
@@ -1674,9 +1674,9 @@ end
subroutine matptm(x,dx,mmx,mx,y,dy,mmy,my,z)
integer dx(2),dy(2)
-integer mmx(1), mmy(1)
-integer mx(1), my(1)
-double precision x(1), y(1),z(1),ddot
+integer mmx(*), mmy(*)
+integer mx(*), my(*)
+double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j
call colmis(mmx,dx(1),dx(2),mx)
@@ -1697,7 +1697,7 @@ end
subroutine rowmis(m,n,p,vec)
integer n,p
-integer m(n,p); integer vec(1)
+integer m(n,p); integer vec(*)
do i = 1,n {
vec(i)=0
# vec(i)=.false.
@@ -1710,7 +1710,7 @@ end
subroutine colmis(m,n,p,vec)
integer n,p
-integer m(n,p); integer vec(1)
+integer m(n,p); integer vec(*)
do j = 1,p {
vec(j)=0
do i = 1,n {
@@ -1721,7 +1721,7 @@ return
end
subroutine daxpy(n,da,dx,incx,dy,incy)
-double precision dx(1),dy(1),da
+double precision dx(*),dy(*),da
integer i,incx,incy,m,mp1,n
if (n>0)
if (da!=0.0d0)
@@ -1760,7 +1760,7 @@ end
subroutine dcopy(n,dx,incx,dy,incy)
-double precision dx(1),dy(1)
+double precision dx(*),dy(*)
integer i,incx,incy,ix,iy,m,mp1,n
if (n>0)
if (incx!=1||incy!=1) {
@@ -1801,7 +1801,7 @@ end
double precision function ddot(n,dx,incx,dy,incy)
-double precision dx(1),dy(1),dtemp
+double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
ddot = 0.0d0
dtemp = 0.0d0
@@ -1840,7 +1840,7 @@ end
double precision function dnrm2(n,dx,incx)
integer nst
-double precision dx(1),cutlo,cuthi,hitest,sum,xmax,zero,one
+double precision dx(*),cutlo,cuthi,hitest,sum,xmax,zero,one
data zero,one/0.0d0,1.0d0/
data cutlo,cuthi/8.232d-11,1.304d19/
if (n<=0)
@@ -1905,7 +1905,7 @@ end
subroutine dscal(n,da,dx,incx)
-double precision da,dx(1)
+double precision da,dx(*)
integer i,incx,m,mp1,n,nincx
if (n>0)
if (incx!=1) {
@@ -1936,7 +1936,7 @@ end
subroutine dswap(n,dx,incx,dy,incy)
-double precision dx(1),dy(1),dtemp
+double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
if (n>0)
if (incx!=1||incy!=1) {
@@ -2000,8 +2000,8 @@ end
subroutine rtod(dx,dy,n)
-real dx(1)
-double precision dy(1)
+real dx(*)
+double precision dy(*)
integer i,m,mp1,n
if (n>0) {
m = mod(n,7)
@@ -2028,8 +2028,8 @@ end
subroutine dtor(dx,dy,n)
-double precision dx(1)
-real dy(1)
+double precision dx(*)
+real dy(*)
integer i,m,mp1,n
if (n>0) {
m = mod(n,7)
@@ -2056,7 +2056,7 @@ end
subroutine drot(n,dx,incx,dy,incy,c,s)
-double precision dx(1),dy(1),dtemp,c,s
+double precision dx(*),dy(*),dtemp,c,s
integer i,incx,incy,ix,iy,n
if (n>0)
if (incx==1&&incy==1)
@@ -2116,7 +2116,7 @@ end
subroutine dqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
integer ldx,n,k,job,info
-double precision x(ldx,1),qraux(1),y(1),qy(1),qty(1),b(1),rsd(1),xb(1)
+double precision x(ldx,*),qraux(*),y(*),qy(*),qty(*),b(*),rsd(*),xb(*)
integer i,j,jj,ju,kp1
double precision ddot,t,temp
logical cb,cqy,cqty,cr,cxb
@@ -2217,7 +2217,7 @@ end
subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
integer ldx,n,p,ldu,ldv,job,info
-double precision x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1)
+double precision x(ldx,*),s(*),e(*),u(ldu,*),v(ldv,*),work(*)
integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit,mm,mm1,mp1,nct,nctp1,ncu,nrt,nrtp1
double precision ddot,t
double precision b,c,cs,el,emm1,f,g,dnrm2,scale,shift,sl,sm,sn,smm1,t1,test,ztest
@@ -2503,7 +2503,7 @@ end
subroutine dtrsl(t,ldt,n,b,job,info)
integer ldt,n,job,info
-double precision t(ldt,1),b(1)
+double precision t(ldt,*),b(*)
double precision ddot,temp
integer which,j,jj
# check for zero diagonal elements.
diff --git a/inst/ratfor/lo.r b/inst/ratfor/lo.r
new file mode 100644
index 0000000..b73d52e
--- /dev/null
+++ b/inst/ratfor/lo.r
@@ -0,0 +1,187 @@
+subroutine lo0(x,y,w,n,d,p,nvmax,span,degree,match,nef,dof,s,var,
+ beta,iv,liv,lv,v,iwork,work)
+integer n,d,p,nvmax,degree,match(*),nef,liv,lv,iv(liv),iwork(*)
+double precision x(n,d),y(n),w(n),span,dof,s(n),var(n),v(lv),work(*)
+double precision beta(p+1)
+#work should be nef*(p+d+8) + 2*p + n +8
+integer qrank
+ call lo1(x,y,w,n,d,p,nvmax,span,degree,match,nef,0,dof,s,var,beta,
+# xin,win,sqwin,sqwini,
+ work(1),work(nef*d+1),work(nef*(d+1)+2),work(nef*(d+2)+2),
+# xqr,qrank,qpivot,qraux,
+ work(nef*(d+3)+2),qrank,iwork(1),work(nef*(p+d+4)+3+p),
+ iv,liv,lv,v,
+ work(nef*(p+d+4)+4+2*p) )
+return
+end
+
+subroutine lo1(x,y,w,n,d,p,nvmax,span,degree,match,nef,nit,dof,s,var,beta,
+ xin,win,sqwin,sqwini,xqr,qrank,qpivot,qraux,
+ iv,liv,lv,v,
+ work)
+integer n,d,p,nvmax,degree,match(*),nef,nit,qrank,qpivot(p+1)
+integer iv(liv),liv,lv
+double precision x(n,d),y(n),w(n),span,dof,s(n),var(n),beta(p+1),
+ xin(nef,d),win(nef+1),sqwin(nef),sqwini(nef),xqr(nef,p+1),
+ qraux(p+1),v(lv),
+ work(*)
+#work should have size n +4*(nef+1)
+call lo2(x,y,w,n,d,p,nvmax,span,degree,match,nef,nit,dof,s,var,beta,
+ xin,win,sqwin,sqwini,xqr,qrank,qpivot,qraux,
+ iv,liv,lv,v,
+ work(1),work(nef+2),work(2*nef+3),work(3*nef+4))
+return
+end
+
+
+
+subroutine lo2(x,y,w,n,d,p,nvmax,span,degree,match,nef,nit,dof,s,var,beta,
+ xin,win,sqwin,sqwini,xqr,qrank,qpivot,qraux,
+ iv,liv,lv,v,
+ levout,sout,yin,work)
+integer n,d,p,nvmax,degree,match(*),nef,nit,qrank,qpivot(p+1)
+integer iv(liv),liv,lv
+double precision x(n,d),y(n),w(n),span,dof,s(n),var(n),beta(p+1),
+ xin(nef,d),win(nef+1),sqwin(nef),sqwini(nef),xqr(nef,p+1),
+ qraux(p+1),v(lv),
+ levout(nef+1), sout(nef+1),yin(nef+1),work(*)
+#work should be length n
+double precision junk, onedm7
+integer job, info
+logical setLf, ifvar
+job=110;info=1
+ifvar=.true.
+onedm7=1d-7
+if(nit<=1){
+ call pck(n,nef,match,w,win)
+ do i=1,nef{
+ if(win(i)>0d0){
+ sqwin(i)=dsqrt(win(i))
+ sqwini(i)=1d0/sqwin(i)
+ }
+ else{
+ sqwin(i)=1d-5
+ sqwini(i)=1d5
+ }
+ }
+ do i=1,n{
+ k=match(i)
+ if(k<=nef){
+ do j=1,d
+ xin(k,j)=x(i,j)
+ for(j=d+1;j<=p;j=j+1)
+ xqr(k,j+1)=x(i,j)
+
+ }
+ }
+ do i=1,nef{
+ xqr(i,1)=sqwin(i)
+ do j=1,d
+ xqr(i,j+1)=xin(i,j)*sqwin(i)
+ for(j=d+2;j<=p+1;j=j+1)
+ xqr(i,j)=xqr(i,j)*sqwin(i)
+ }
+ for(j=1;j<=p+1;j=j+1)
+ qpivot(j)=j
+ call dqrdca(xqr,nef,nef,p+1,qraux,qpivot,work,qrank,onedm7)
+ setLf = (nit==1)
+ call lowesd(106,iv,liv,lv,v,d,nef,span,degree,nvmax,setLf)
+ v(2)=span/5d0
+ }
+do i=1,n
+ work(i)=y(i)*w(i)
+call pck(n,nef,match,work,yin)
+ do i=1,nef
+ yin(i)=yin(i)*sqwini(i)*sqwini(i)
+
+
+if(nit<=1)call lowesb(xin,yin,win,levout,ifvar,iv,liv,lv,v)
+ else call lowesr(yin,iv,liv,lv,v)
+call lowese(iv,liv,lv,v,nef,xin,sout)
+
+#now remove the parametric piece
+do i=1,nef
+ sout(i)=sout(i)*sqwin(i)
+call dqrsl(xqr,nef,nef,qrank,qraux,sout,work(1),work(1),beta,
+ sout,work(1),job,info)
+#####dqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
+do i=1,nef
+ sout(i)=sout(i)*sqwini(i)
+
+#now clean up
+
+if(nit<=1){
+#get rid of the parametric component of the leverage
+ job=10000
+ for(j=1;j<=p+1;j=j+1){
+ do i=1,nef
+ work(i)=0d0
+ work(j)=1d0
+ call dqrsl(xqr,nef,nef,qrank,qraux,work,var,junk,junk,
+ junk,junk,job,info)
+ do i=1,nef
+ levout(i)=levout(i) - var(i)**2
+ }
+ dof=0d0
+ do i=1,nef {
+ if(win(i)>0d0) {
+ levout(i)=levout(i)/win(i)
+ }
+ else {levout(i)=0d0}
+ }
+ do i=1,nef {dof=dof+levout(i)*win(i)}
+ call unpck(n,nef,match,levout,var)
+ for(j=1;j<=p+1;j=j+1){work(j)=beta(j)}
+ for(j=1;j<=p+1;j=j+1){beta(qpivot(j))=work(j)}
+ }
+call unpck(n,nef,match,sout,s)
+return
+end
+subroutine pck(n,p,match,x,xbar)
+integer match(n),p,n
+double precision x(n),xbar(n)
+ do i=1,p
+ xbar(i)=0d0
+ do i=1,n
+ xbar(match(i))=xbar(match(i))+x(i)
+return
+end
+
+subroutine suff(n,p,match,x,y,w,xbar,ybar,wbar,work)
+integer match(n),p,n
+double precision x(n),xbar(n),y(n),ybar(n),w(n),wbar(n),work(n)
+call pck(n,p,match,w,wbar)
+do i=1,n
+ xbar(match(i))=x(i)
+do i=1,n
+ work(i)=y(i)*w(i)
+call pck(n,p,match,work,ybar)
+ do i=1,p{
+ if(wbar(i)>0d0)
+ ybar(i)=ybar(i)/wbar(i)
+ else ybar(i)=0d0
+
+ }
+return
+end
+subroutine unpck(n,p,match,xbar,x)
+integer match(n),p,n
+double precision x(n),xbar(p+1)
+if(p<n)xbar(p+1)=0d0
+do i = 1,n
+ x(i)=xbar(match(i))
+return
+end
+double precision function dwrss(n,y,eta,w)
+integer n
+double precision y(n),w(n),wtot,wsum,work,eta(n)
+wsum=0d0
+wtot=0d0
+do i = 1,n{
+ work=y(i)-eta(i)
+ wsum=wsum+w(i)*work*work
+ wtot=wtot+w(i)
+}
+if (wtot > 0d0) {dwrss=wsum/wtot} else {dwrss=0d0}
+return
+end
diff --git a/inst/ratfor/splsm.r b/inst/ratfor/splsm.r
new file mode 100644
index 0000000..bb5d241
--- /dev/null
+++ b/inst/ratfor/splsm.r
@@ -0,0 +1,239 @@
+subroutine sknotl(x,n,knot,k)
+implicit double precision(a-h,o-z)
+double precision x(n),knot(n+6),a1,a2,a3,a4
+integer n,k,ndk,j
+
+
+ # Allocate knots acording to the cutoffs given below
+
+
+ # Cutoff constants
+
+ a1 = log(50d0)/log(2d0) ; a2 = log(100d0)/log(2d0)
+ a3 = log(140d0)/log(2d0) ; a4 = log(200d0)/log(2d0)
+
+ # Cutoff Criteria
+
+ if(n<50) { ndk = n }
+ else if (n>=50 & n<200) { ndk = 2.**(a1+(a2-a1)*(n-50.)/150.) }
+ else if (n>=200 & n<800) { ndk = 2.**(a2+(a3-a2)*(n-200.)/600.) }
+ else if (n>=800 & n<3200) { ndk = 2.**(a3+(a4-a3)*(n-800.)/2400.) }
+ else if (n>=3200) { ndk = 200. + float(n-3200)**.2 }
+
+ k = ndk + 6
+
+
+ # Allocate Knots ( note no account is taken of any weighting vector )
+
+ do j=1,3 { knot(j) = x(1) }
+ do j=1,ndk { knot(j+3) = x( 1 + (j-1)*(n-1)/(ndk-1) ) }
+ do j=1,3 { knot(ndk+3+j) = x(n) }
+
+return
+end
+subroutine splsm(x,y,w,n,match,nef,spar,dof,smo,s0,cov,ifcov,work)
+#This subroutine performs a smoothing spline fit
+# This was written by Trevor Hastie in 1990
+# It has been modified from the S3 version by Trevor Hastie
+# in July 2004, to accommodate the modified sbart routine in R
+# and also to accommodate only the gam bakfit routine.
+# Note that spar has changed, and we change it here to conform with
+# the smooth.spline routine in R
+#All arguments are either double precision or integer
+#INPUT
+#
+#x double length n ; x variable for smoothing
+#y double length n ; y variable for smoothing
+#w double length n ; weights for smoothing, > 0
+#n integer length above
+#match integer length n -- in S language x[i] == sort(unique(x)[match[i]]
+# match is produced by subroutine namat
+#nef number of unique elements in x; so match has values between 1 and nef+1
+# missing data are given the match number nef+1
+#spar double smoothing parameter -1.5 <spar<1.5; default is 1
+#dof double equivalent degrees of freedom
+# if dof is 0, spar is used
+# if 0< dof <1, dof = 1
+# if dof >=1, dof is used
+# note: dof does not use the constant term
+#ifcov integer if 1, the unscaled variance information is computed
+#work double workspace of length (10+2*4)*(nef+2)+5*nef+n+15
+#
+#OUTPUT
+#
+#x,y,w,n,match,nef are untouched
+#spar if dof > 1, then spar is that which achieves dof
+#dof the dof of the fitted smooth. Note: even if dof was given
+# as 4, it will be returned as say 3.995 which is what
+# spar produces
+#smo double length n the fitted values, with weighted average 0
+#s0 double weighted mean of y
+#cov double length nef the unscaled variance elements for the NONLINEAR
+# and UNIQUE part of smo, in the order of sort(unique(x))
+# cov is lev(i)/w(i) -h(i)/w where h(i) is the hat element from
+# the simple weighted least squares fit. This is passed on
+# to bakfit and used in gamcov
+#
+# splsm calls (eventually after some memory management dummy calls)
+# sbart, the spline routine of Finbarr O'Sullivan, slightly modified
+# by Trevor Hastie, 8/2/89
+
+implicit double precision(a-h,o-z)
+double precision x(*),y(*),w(*),spar,dof,smo(*),s0,cov(*),work(*)
+integer n,match(*),nef
+integer ifcov
+# work should be (10+2*ld4)*nk+5*nef+n+15 double
+# ld4 =4 nk<= nef+2
+call splsm1(x,y,w,n,match,nef,spar,dof,smo,s0,cov,ifcov,
+# xin(nef+1),yin(nef+1), win(nef+1), knot(n+6),
+ work(1), work(nef+2),work(2*nef+3),work(3*nef+4),
+ work(3*nef+n+10))
+return
+end
+
+subroutine splsm1(x,y,w,n,match,nef,spar,dof,smo,s0,lev,ifcov,
+ xin,yin,win,knot,
+ work)
+implicit double precision(a-h,o-z)
+double precision x(*),y(*),w(*),spar,dof,smo(*),s0,lev(*),work(*)
+integer n,match(*),nef
+integer ifcov
+double precision xin(nef+1),yin(nef+1),win(nef+1),knot(nef+6)
+integer nk,ldnk,ld4,k
+double precision xmin,xrange
+call suff(n,nef,match,x,y,w,xin,yin,win,work(1))
+xmin=xin(1)
+xrange=xin(nef)-xin(1)
+do i=1,nef {xin(i)=(xin(i)-xmin)/xrange}
+call sknotl(xin,nef,knot,k)
+nk=k-4
+ld4=4
+ldnk=1 # p21p nd ldnk is not used
+call splsm2(x,y,w,n,match,nef,spar,dof,smo,s0,lev,ifcov,
+ xin,yin,win,knot,
+# coef(nk),sout(nef+1), levout(nef+1), xwy(nk),
+# hs0(nk), hs1(nk), hs2(nk),
+# hs3(nk),
+# sg0(nk), sg1(nk), sg2(nk),
+# sg3(nk),
+# abd(ld4,nk), p1ip(ld4,nk),
+# p2ip(ldnk,nk)
+ work(1), work(nk+1), work(nk+nef+2),work(nk+2*nef+3),
+ work(2*nk+2*nef+3),work(3*nk+2*nef+3),work(4*nk+2*nef+3),
+ work(5*nk+2*nef+3),
+ work(6*nk+2*nef+3),work(7*nk+2*nef+3),work(8*nk+2*nef+3),
+ work(9*nk+2*nef+3),
+ work(10*nk+2*nef+3),work((10+ld4)*nk+2*nef+3),
+ work((10+2*ld4)*nk+2*nef+3),
+ ld4,ldnk,nk)
+
+return
+end
+subroutine splsm2(x,y,w,n,match,nef,spar,dof,smo,s0,lev,ifcov,
+ xin,yin,win,knot,
+ coef,sout,levout,xwy,
+ hs0,hs1,hs2,hs3,
+ sg0,sg1,sg2,sg3,
+ abd,p1ip,p2ip,ld4,ldnk,nk)
+implicit double precision(a-h,o-z)
+double precision x(*),y(*),w(*),spar,dof,smo(*),s0,lev(*)
+integer n,match(*),nef
+integer nk,ldnk,ld4
+integer ifcov
+double precision xin(nef+1),yin(nef+1),win(nef+1),knot(nk+4)
+double precision coef(nk),sout(nef+1),levout(nef+1),xwy(nk),
+ hs0(nk),hs1(nk),hs2(nk),hs3(nk),
+ sg0(nk),sg1(nk),sg2(nk),sg3(nk),
+ abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,*)
+# local variables
+integer ispar,icrit,isetup,ier
+double precision lspar,uspar,tol,penalt,
+ sumwin,dofoff,crit,xbar,dsum,xsbar
+double precision yssw, eps
+integer maxit
+# yssw is an additional parameter introduced in R version of sbart
+double precision wmean
+crit=0d0
+# Note we only allow limited options here
+if(dof <= 0d0){
+# use spar
+ ispar=1
+ icrit=3
+ dofoff=0d0
+ }
+else{
+ if( dof < 1d0 )dof=1d0
+ ispar=0
+ icrit=3
+ dofoff=dof+1d0
+}
+#Here we set some default parameters similar to the smooth.spline in R
+isetup=0
+ier=1
+penalt=1d0
+lspar= -1.5
+uspar= 1.5
+tol=1d-4
+eps=2d-8
+maxit=200
+do i=1,nef
+ sout(i)=yin(i)*yin(i)
+sumwin=0d0
+do i=1,nef
+ sumwin=sumwin+win(i)
+yssw=wmean(nef,sout,win)
+s0=wmean(n,y,w)
+# which should be equal to wmean(nef,yin,win)
+yssw=yssw*(sumwin-s0*s0)
+
+call sbart(penalt,dofoff,xin,yin,win,yssw,nef,knot,nk,
+ coef,sout,levout,crit,
+ icrit,spar,ispar,maxit,
+ lspar,uspar,tol,eps,
+ isetup,
+ xwy,
+ hs0,hs1,hs2,hs3,
+ sg0,sg1,sg2,sg3,
+ abd,p1ip,p2ip,ld4,ldnk,ier)
+#return
+#now clean up
+do i=1,nef {
+ win(i)=win(i)*win(i) #we sqrted them in sbart
+ }
+sbar=wmean(nef,sout,win)
+xbar=wmean(nef,xin,win)
+
+
+#now remove the linear leverage from the leverage for the smooths
+# will be altered at this stage
+do i=1,nef {lev(i)=(xin(i)-xbar)*sout(i) }
+xsbar=wmean(nef,lev,win)
+do i=1,nef {lev(i)=(xin(i)-xbar)**2 }
+dsum=wmean(nef,lev,win)
+do i=1,nef {
+ if(win(i)>0d0) {
+ lev(i)=levout(i)/win(i)-1d0/sumwin -lev(i)/(sumwin*dsum)
+ }
+ else {lev(i)=0d0}
+ }
+dof=0d0
+do i=1,nef {dof=dof+lev(i)*win(i)}
+dof=dof+1d0
+do i=1,nef
+ sout(i)=sout(i)-sbar -(xin(i)-xbar)*xsbar/dsum
+call unpck(n,nef,match,sout,smo)
+return
+end
+
+double precision function wmean(n,y,w)
+integer n
+double precision y(n),w(n),wtot,wsum
+wtot=0d0
+wsum=0d0
+do i=1,n{
+ wsum=wsum+y(i)*w(i)
+ wtot=wtot+w(i)
+}
+if(wtot > 0d0) {wmean=wsum/wtot} else {wmean=0d0}
+return
+end
diff --git a/man/predict.gam.Rd b/man/predict.gam.Rd
index 4c85d58..399c2c6 100644
--- a/man/predict.gam.Rd
+++ b/man/predict.gam.Rd
@@ -4,7 +4,7 @@
\description{Obtains predictions and optionally estimates standard errors of
those predictions from a fitted generalized additive model object.}
\usage{
-predict.gam(object, newdata, type, dispersion, se.fit = FALSE,na.action, terms,\dots)
+\method{predict}{gam}(object, newdata, type, dispersion, se.fit = FALSE,na.action, terms,\dots)
}
\arguments{
\item{object}{
diff --git a/src/Makevars.win b/src/Makevars.win
index 34a28f2..ea686ad 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -1 +1 @@
-PKG_LIBS = $(BLAS_LIBS) $(FLIBS)
+PKG_LIBS = $(BLAS_LIBS) $(FLIBS)
diff --git a/src/backfit.f b/src/backfit.f
index 4d027be..9293b85 100644
--- a/src/backfit.f
+++ b/src/backfit.f
@@ -4,9 +4,9 @@ C Output from Public domain Ratfor, version 1.0
implicit double precision(a-h,o-z)
logical ifvar
integer npetc(7),iter
- integer n,p,q,which(1),match(1),nef(1),nit,maxit,qrank,qpivot(1)
- double precision x(1),y(1),w(1),spar(1),dof(1), etal(1),s(1),eta(1
- *),beta(1),var(1),tol, qr(1),qraux(1),effect(1),work(1)
+ integer n,p,q,which(*),match(*),nef(*),nit,maxit,qrank,qpivot(*)
+ double precision x(*),y(*),w(*),spar(*),dof(*), etal(*),s(*),eta(*
+ *),beta(*),var(*),tol, qr(*),qraux(*),effect(*),work(*)
n=npetc(1)
p=npetc(2)
q=npetc(3)
@@ -33,8 +33,8 @@ C Output from Public domain Ratfor, version 1.0
logical ifvar
integer n,p,q,which(q),match(n,q),nef(q),nit,maxit,qrank,qpivot(p)
double precision x(n,p),y(n),w(n),spar(q),dof(q), etal(n),s(n,q),e
- *ta(n),beta(p),var(n,q),tol, qr(n,p),qraux(p),effect(n),work(1)
- double precision z(1),old(1),dwrss,ratio
+ *ta(n),beta(p),var(n,q),tol, qr(n,p),qraux(p),effect(n),work(*)
+ double precision z(*),old(*),dwrss,ratio
double precision sqwt(n),sqwti(n)
logical anyzwt
double precision deltaf, normf,onedm7
diff --git a/src/backlo.f b/src/backlo.f
index 5a516c1..d5ed52a 100644
--- a/src/backlo.f
+++ b/src/backlo.f
@@ -3,9 +3,9 @@ C Output from Public domain Ratfor, version 1.0
*var,dof, qr,qraux,qpivot,effect,iv,v,iwork,work)
implicit double precision(a-h,o-z)
integer n,p,q,nit,maxit,qrank
- integer npetc(7),wddnfl(1),match(1),qpivot(1),iv(1),iwork(1)
- double precision x(1),y(1),w(1),spatol(1), etal(1),s(1),eta(1),bet
- *a(1),var(1),dof(1), qr(1),qraux(1),v(1),effect(1),work(1)
+ integer npetc(7),wddnfl(*),match(*),qpivot(*),iv(*),iwork(*)
+ double precision x(*),y(*),w(*),spatol(*), etal(*),s(*),eta(*),bet
+ *a(*),var(*),dof(*), qr(*),qraux(*),v(*),effect(*),work(*)
n=npetc(1)
p=npetc(2)
q=npetc(3)
@@ -25,12 +25,12 @@ C Output from Public domain Ratfor, version 1.0
*t,effect,z,old,sqwt,sqwti, iv,liv,lv,v,nvmax,iwork,work)
implicit double precision(a-h,o-z)
integer n,p,q,which(q),dwhich(q),pwhich(q),degree(q),match(n,q),ne
- *f(q),nit, maxit,qrank,qpivot(p),iv(1),liv(q),lv(q),nvmax(q),iwork(
+ *f(q),nit, maxit,qrank,qpivot(p),iv(*),liv(q),lv(q),nvmax(q),iwork(
*q)
double precision x(n,p),y(n),w(n),span(q),dof(q), etal(n),s(n,q),e
- *ta(n),beta(p),var(n,q),tol, qr(n,p),qraux(p),v(1),effect(n),work(1
+ *ta(n),beta(p),var(n,q),tol, qr(n,p),qraux(p),v(*),effect(n),work(*
*)
- double precision z(1),old(1),dwrss,ratio
+ double precision z(*),old(*),dwrss,ratio
double precision sqwt(n),sqwti(n)
logical anyzwt
double precision deltaf, normf,onedm7
diff --git a/src/linear.f b/src/linear.f
index d927aca..afc4736 100644
--- a/src/linear.f
+++ b/src/linear.f
@@ -1,9 +1,9 @@
-C Output from Public domain Ratfor, version 1.01
+C Output from Public domain Ratfor, version 1.0
subroutine dqrls(x,dx,pivot,qraux,y,dy,beta,res,qt,tol,scrtch,rank
*)
- integer pivot(1),dx(2),dy(2),rank
- double precision x(1), qraux(1), y(1), beta(1),res(1),qt(1),tol(1)
- *, scrtch(1)
+ integer pivot(*),dx(2),dy(2),rank
+ double precision x(*), qraux(*), y(*), beta(*),res(*),qt(*),tol(*)
+ *, scrtch(*)
integer n,p,q,kn,kp,k,info
n=dx(1)
p=dx(2)
@@ -25,7 +25,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dqrsl1(qr,dq,qra,rank,y,k,qy,qb,job,info)
- double precision qr(1),qra(1),y(1),qy(1),qb(1)
+ double precision qr(*),qra(*),y(*),qy(*),qb(*)
integer dq(2),job,k,rank
integer n,kn,kb,j
double precision ourqty(1), ourqy(1), ourb(1), ourrsd(1), ourxb(1)
@@ -104,8 +104,8 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dqr(x,dx,pivot,qraux,tol,scrtch,rank)
- integer pivot(1),dx(2),rank
- double precision x(1), qraux(1), tol(1), scrtch(1)
+ integer pivot(*),dx(2),rank
+ double precision x(*), qraux(*), tol(*), scrtch(*)
integer n,p
n=dx(1)
p=dx(2)
@@ -114,8 +114,8 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine dqrdca(x,ldx,n,p,qraux,jpvt,work,rank,eps)
integer ldx,n,p,rank
- integer jpvt(1)
- double precision x(ldx,1),qraux(1),work(1),eps
+ integer jpvt(*)
+ double precision x(ldx,*),qraux(*),work(*),eps
integer j,jj,jp,l,lup,curpvt
double precision dnrm2,tt
double precision ddot,nrmxl,t,ww
@@ -394,8 +394,8 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine chol(a,p,work,jpvt,job,info)
- integer p,jpvt(1),job,info(1)
- double precision a(p,1),work(1)
+ integer p,jpvt(*),job,info(*)
+ double precision a(p,*),work(*)
integer i,j
j =2
23124 if(.not.(j.le.p))goto 23126
@@ -426,7 +426,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine crs(x,dmx,matz,w,z,fv1,fv2,ierr)
- double precision x(1),w(1),z(1),fv1(1),fv2(1)
+ double precision x(*),w(*),z(*),fv1(*),fv2(*)
integer dmx(2),nx,nv,ierr,matz
nx=dmx(1)
nv=dmx(2)
@@ -434,9 +434,9 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dqrls2(x,dx,pivot,qraux,y,dy,beta,res,qt,scrtch,eps)
- integer pivot(1),dx(2),dy(2)
- double precision x(1), qraux(1), y(1), beta(1),res(1),qt(1), scrtc
- *h(1),eps
+ integer pivot(*),dx(2),dy(2)
+ double precision x(*), qraux(*), y(*), beta(*),res(*),qt(*), scrtc
+ *h(*),eps
integer n,p,q,kn,kp,k,info,rank
n=dx(1)
p=dx(2)
@@ -456,7 +456,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dsvdc1(x,dmx,job,work,e,s,u,v,info)
- double precision x(1),work(1),s(1),e(1),u(1),v(1)
+ double precision x(*),work(*),s(*),e(*),u(*),v(*)
integer dmx(2),nx,nv,job,info
nx=dmx(1)
nv=dmx(2)
@@ -1783,7 +1783,7 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine dmatp(x,dx,y,dy,z)
integer dx(2),dy(2)
- double precision x(1), y(1),z(1),ddot
+ double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j
n=dx(1)
p=dx(2)
@@ -1805,7 +1805,7 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine dmatpt(x,dx,y,dy,z)
integer dx(2),dy(2)
- double precision x(1), y(1),z(1),ddot
+ double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j,ii
n=dx(1)
p=dx(2)
@@ -1829,9 +1829,9 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine matpm(x,dx,mmx,mx,y,dy,mmy,my,z)
integer dx(2),dy(2)
- integer mmx(1), mmy(1)
- integer mx(1), my(1)
- double precision x(1), y(1),z(1),ddot
+ integer mmx(*), mmy(*)
+ integer mx(*), my(*)
+ double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j
n=dx(1)
p=dx(2)
@@ -1857,9 +1857,9 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine matptm(x,dx,mmx,mx,y,dy,mmy,my,z)
integer dx(2),dy(2)
- integer mmx(1), mmy(1)
- integer mx(1), my(1)
- double precision x(1), y(1),z(1),ddot
+ integer mmx(*), mmy(*)
+ integer mx(*), my(*)
+ double precision x(*), y(*),z(*),ddot
integer n,p,q,i,j
call colmis(mmx,dx(1),dx(2),mx)
call colmis(mmy,dy(1),dy(2),my)
@@ -1888,7 +1888,7 @@ C Output from Public domain Ratfor, version 1.01
subroutine rowmis(m,n,p,vec)
integer n,p
integer m(n,p)
- integer vec(1)
+ integer vec(*)
do23659 i = 1,n
vec(i)=0
do23661 j = 1,p
@@ -1904,7 +1904,7 @@ C Output from Public domain Ratfor, version 1.01
subroutine colmis(m,n,p,vec)
integer n,p
integer m(n,p)
- integer vec(1)
+ integer vec(*)
do23665 j = 1,p
vec(j)=0
do23667 i = 1,n
@@ -1918,7 +1918,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine daxpy(n,da,dx,incx,dy,incy)
- double precision dx(1),dy(1),da
+ double precision dx(*),dy(*),da
integer i,incx,incy,m,mp1,n
if(n.gt.0)then
if(da.ne.0.0d0)then
@@ -1962,7 +1962,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dcopy(n,dx,incx,dy,incy)
- double precision dx(1),dy(1)
+ double precision dx(*),dy(*)
integer i,incx,incy,ix,iy,m,mp1,n
if(n.gt.0)then
if(incx.ne.1.or.incy.ne.1)then
@@ -2007,7 +2007,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
double precision function ddot(n,dx,incx,dy,incy)
- double precision dx(1),dy(1),dtemp
+ double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
ddot = 0.0d0
dtemp = 0.0d0
@@ -2052,7 +2052,7 @@ C Output from Public domain Ratfor, version 1.01
end
double precision function dnrm2(n,dx,incx)
integer nst
- double precision dx(1),cutlo,cuthi,hitest,sum,xmax,zero,one
+ double precision dx(*),cutlo,cuthi,hitest,sum,xmax,zero,one
data zero,one/0.0d0,1.0d0/
data cutlo,cuthi/8.232d-11,1.304d19/
if(n.le.0)then
@@ -2128,7 +2128,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dscal(n,da,dx,incx)
- double precision da,dx(1)
+ double precision da,dx(*)
integer i,incx,m,mp1,n,nincx
if(n.gt.0)then
if(incx.ne.1)then
@@ -2162,7 +2162,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dswap(n,dx,incx,dy,incy)
- double precision dx(1),dy(1),dtemp
+ double precision dx(*),dy(*),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
if(n.gt.0)then
if(incx.ne.1.or.incy.ne.1)then
@@ -2230,8 +2230,8 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine rtod(dx,dy,n)
- real dx(1)
- double precision dy(1)
+ real dx(*)
+ double precision dy(*)
integer i,m,mp1,n
if(n.gt.0)then
m = mod(n,7)
@@ -2259,8 +2259,8 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine dtor(dx,dy,n)
- double precision dx(1)
- real dy(1)
+ double precision dx(*)
+ real dy(*)
integer i,m,mp1,n
if(n.gt.0)then
m = mod(n,7)
@@ -2288,7 +2288,7 @@ C Output from Public domain Ratfor, version 1.01
return
end
subroutine drot(n,dx,incx,dy,incy,c,s)
- double precision dx(1),dy(1),dtemp,c,s
+ double precision dx(*),dy(*),dtemp,c,s
integer i,incx,incy,ix,iy,n
if(n.gt.0)then
if(incx.eq.1.and.incy.eq.1)then
@@ -2349,8 +2349,8 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine dqrsl(x,ldx,n,k,qraux,y,qy,qty,b,rsd,xb,job,info)
integer ldx,n,k,job,info
- double precision x(ldx,1),qraux(1),y(1),qy(1),qty(1),b(1),rsd(1),x
- *b(1)
+ double precision x(ldx,*),qraux(*),y(*),qy(*),qty(*),b(*),rsd(*),x
+ *b(*)
integer i,j,jj,ju,kp1
double precision ddot,t,temp
logical cb,cqy,cqty,cr,cxb
@@ -2475,7 +2475,7 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine dsvdc(x,ldx,n,p,s,e,u,ldu,v,ldv,work,job,info)
integer ldx,n,p,ldu,ldv,job,info
- double precision x(ldx,1),s(1),e(1),u(ldu,1),v(ldv,1),work(1)
+ double precision x(ldx,*),s(*),e(*),u(ldu,*),v(ldv,*),work(*)
integer i,iter,j,jobu,k,kase,kk,l,ll,lls,lm1,lp1,ls,lu,m,maxit,mm,
*mm1,mp1,nct,nctp1,ncu,nrt,nrtp1
double precision ddot,t
@@ -2853,7 +2853,7 @@ C Output from Public domain Ratfor, version 1.01
end
subroutine dtrsl(t,ldt,n,b,job,info)
integer ldt,n,job,info
- double precision t(ldt,1),b(1)
+ double precision t(ldt,*),b(*)
double precision ddot,temp
integer which,j,jj
do24063 info = 1,n
diff --git a/src/lo.f b/src/lo.f
index 7e124ff..824c910 100644
--- a/src/lo.f
+++ b/src/lo.f
@@ -1,9 +1,9 @@
C Output from Public domain Ratfor, version 1.0
subroutine lo0(x,y,w,n,d,p,nvmax,span,degree,match,nef,dof,s,var,
*beta,iv,liv,lv,v,iwork,work)
- integer n,d,p,nvmax,degree,match(1),nef,liv,lv,iv(liv),iwork(1)
+ integer n,d,p,nvmax,degree,match(*),nef,liv,lv,iv(liv),iwork(*)
double precision x(n,d),y(n),w(n),span,dof,s(n),var(n),v(lv),work(
- *1)
+ **)
double precision beta(p+1)
integer qrank
call lo1(x,y,w,n,d,p,nvmax,span,degree,match,nef,0,dof,s,var,beta,
@@ -15,11 +15,11 @@ C Output from Public domain Ratfor, version 1.0
subroutine lo1(x,y,w,n,d,p,nvmax,span,degree,match,nef,nit,dof,s,v
*ar,beta, xin,win,sqwin,sqwini,xqr,qrank,qpivot,qraux, iv,liv,lv,v,
* work)
- integer n,d,p,nvmax,degree,match(1),nef,nit,qrank,qpivot(p+1)
+ integer n,d,p,nvmax,degree,match(*),nef,nit,qrank,qpivot(p+1)
integer iv(liv),liv,lv
double precision x(n,d),y(n),w(n),span,dof,s(n),var(n),beta(p+1),
*xin(nef,d),win(nef+1),sqwin(nef),sqwini(nef),xqr(nef,p+1), qraux(p
- *+1),v(lv), work(1)
+ *+1),v(lv), work(*)
call lo2(x,y,w,n,d,p,nvmax,span,degree,match,nef,nit,dof,s,var,bet
*a, xin,win,sqwin,sqwini,xqr,qrank,qpivot,qraux, iv,liv,lv,v, work(
*1),work(nef+2),work(2*nef+3),work(3*nef+4))
@@ -28,11 +28,11 @@ C Output from Public domain Ratfor, version 1.0
subroutine lo2(x,y,w,n,d,p,nvmax,span,degree,match,nef,nit,dof,s,v
*ar,beta, xin,win,sqwin,sqwini,xqr,qrank,qpivot,qraux, iv,liv,lv,v,
* levout,sout,yin,work)
- integer n,d,p,nvmax,degree,match(1),nef,nit,qrank,qpivot(p+1)
+ integer n,d,p,nvmax,degree,match(*),nef,nit,qrank,qpivot(p+1)
integer iv(liv),liv,lv
double precision x(n,d),y(n),w(n),span,dof,s(n),var(n),beta(p+1),
*xin(nef,d),win(nef+1),sqwin(nef),sqwini(nef),xqr(nef,p+1), qraux(p
- *+1),v(lv), levout(nef+1), sout(nef+1),yin(nef+1),work(1)
+ *+1),v(lv), levout(nef+1), sout(nef+1),yin(nef+1),work(*)
double precision junk, onedm7
integer job, info
logical setlf, ifvar
diff --git a/src/splsm.f b/src/splsm.f
index 5d2a480..eb9ceb0 100644
--- a/src/splsm.f
+++ b/src/splsm.f
@@ -43,8 +43,8 @@ C Output from Public domain Ratfor, version 1.0
end
subroutine splsm(x,y,w,n,match,nef,spar,dof,smo,s0,cov,ifcov,work)
implicit double precision(a-h,o-z)
- double precision x(1),y(1),w(1),spar,dof,smo(1),s0,cov(1),work(1)
- integer n,match(1),nef
+ double precision x(*),y(*),w(*),spar,dof,smo(*),s0,cov(*),work(*)
+ integer n,match(*),nef
integer ifcov
call splsm1(x,y,w,n,match,nef,spar,dof,smo,s0,cov,ifcov, work(1),
*work(nef+2),work(2*nef+3),work(3*nef+4), work(3*nef+n+10))
@@ -53,8 +53,8 @@ C Output from Public domain Ratfor, version 1.0
subroutine splsm1(x,y,w,n,match,nef,spar,dof,smo,s0,lev,ifcov, xin
*,yin,win,knot, work)
implicit double precision(a-h,o-z)
- double precision x(1),y(1),w(1),spar,dof,smo(1),s0,lev(1),work(1)
- integer n,match(1),nef
+ double precision x(*),y(*),w(*),spar,dof,smo(*),s0,lev(*),work(*)
+ integer n,match(*),nef
integer ifcov
double precision xin(nef+1),yin(nef+1),win(nef+1),knot(nef+6)
integer nk,ldnk,ld4,k
@@ -82,14 +82,14 @@ C Output from Public domain Ratfor, version 1.0
*,yin,win,knot, coef,sout,levout,xwy, hs0,hs1,hs2,hs3, sg0,sg1,sg2,
*sg3, abd,p1ip,p2ip,ld4,ldnk,nk)
implicit double precision(a-h,o-z)
- double precision x(1),y(1),w(1),spar,dof,smo(1),s0,lev(1)
- integer n,match(1),nef
+ double precision x(*),y(*),w(*),spar,dof,smo(*),s0,lev(*)
+ integer n,match(*),nef
integer nk,ldnk,ld4
integer ifcov
double precision xin(nef+1),yin(nef+1),win(nef+1),knot(nk+4)
double precision coef(nk),sout(nef+1),levout(nef+1),xwy(nk), hs0(n
*k),hs1(nk),hs2(nk),hs3(nk), sg0(nk),sg1(nk),sg2(nk),sg3(nk), abd(l
- *d4,nk),p1ip(ld4,nk),p2ip(ldnk,1)
+ *d4,nk),p1ip(ld4,nk),p2ip(ldnk,*)
integer ispar,icrit,isetup,ier
double precision lspar,uspar,tol,penalt, sumwin,dofoff,crit,xbar,d
*sum,xsbar
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-gam.git
More information about the debian-science-commits
mailing list