GAIM is an Interactive Program for fitting Generalized Additive
Models.
GAIM is constantly under development by Trevor Hastie and Rob
Tibshirani, so if you have any problems with the program or suggested
improvements, please contact us. Trevor Hastie- AT&T Bell Labs,
2C-261, 600 Mountain Ave, Murray Hill, NJ07974. (201) 582-5647. Rob
Tibshirani- Dept. of Preventive Medicine and Biostatistics, U. of
Toronto, Toronto Canada M5S1A8, (416)-978-4642.
The program is written in MORTRAN, a FORTRAN preprocessor. The
distribution tape has both a FORTRAN and a MORTRAN version. However,
the FORTRAN produced by MORTRAN is rather messy.
ref: Hastie, T and Tibshirani, R (1986), Generalized Additive Models,
Statistical Science, 1, 297-318
Hastie, T. and Tibshirani, R. (1987), Generalized Additive Models;
some Applications, JASA, June 1987.
pause
LIMITATIONS
___________
The current limits are set at 1000 observations
and 20 variables. This can easily be changed (by you).
SUMMARY of CAPABILITIES
_______________________
GAIM can fit additive models to GAUSSIAN, BINOMIAL, and MULTINOMIAL
(proportional odds) data, as well as CENSORED survival type data and
MATCHED CASE-CONTROL data. For the binomial and multinomial, the data
can be grouped or ungrouped. Users can specifiy their own models by
writing their own versions of the link, deviance, inverse link, weight
and derivative routines. The models can range from all terms linear to
all terms "smooth", although typically one fits some of each. One can
generate dummy variables quite easily for the levels of a categorical
variable. Using the SPECIAL option, you can get confidence bands (+-2
sd) for the fitted functions, as well as variances and covariances for
the linear terms (the current implementation does not allow this for
the MULTINOMIAL, COX or MATCHED CASE-CONTROL Model)
The program can deal with MISSING data explicitly, and
implicitly one can fit interactions between discrete and smooth terms.
The program allows one to build models in a natural GLIM type way by
adding to or deleteing from a current model. After each fit, a page of
statistics and details are printed giving slopes, deviance etc. The
program allows one to print the fitted functions to a PLOT file, which
can then be independently plotted. At present, the subroutine DUMPS
creates a file for each plot, which can then be manipulated.
pause
DATA ENTRY
__________
The program expects 2 data files in a UNIX
directory, and will prompt for the name of this
directory.
Unit number 1 reads from a file called "name/data" which
has a N x P block of data: N observations on P variables in free
format. One observation, (p numbers) per line.
Unit number 11 reads
from a "name/datad" file (Data description). The first line
contains N and P. P lines follow, each line containing a variable
name.
In addition, this file (helpunix) should be in the gaim directory.
Unit 2 is the PLOT output, and UNIT 3 gives fitted values when these
are requested.
pause
some details...
ADD, DELETE, MODIFY or RUN
__________________________
These options are
rather obvious, and allow for model building.
ADD--displays a list of the variables NOT currently in the model, and
allows the user to select one. Upon selection, you are requested to
give a span which can be any real number. The following actions are
taken
SPAN ACTION
__________________________________________________________
s < 0 | A cubic Spline Smoother is used with lamda =-span s = 0 | The span is chosen by cross-validation
0<s<1 | 100s% of the data is used in each neighborhood
| except at the ends, which drop down to 100s/2%
s=1 | This fits a linear fit for this term.
s=2,3,..| Categorical variable with 2,3... categories. If there
| are less then s categories, it allows you to exit or
| continue. If there are more than s, it allows you to exit
| or combine the last additional categories.
The categorical option creates k-1 dummy variables if there are k
categories. The first category gets no dummy, and corresponds to the
constant term in the model. All smooth fits are standardized to have
mean zero---linear fits are not.
pause
DELETE-- This allows you to delete terms in the model.
MODIFY-- This allows you to change the thresh-hold for the BACKFITTING
algorithm (set at !RSSOLD-RSSNEW!/RSSOLD < THBAK=0.0001 ) or the LOCAL
SCORING algorithm (THGAM = 0.0001).
It also allow you to alter the spans for variables. Do not try to
change a continuous variable to a discrete one though.
RUN-- Fits the current model.
pause
MISSING DATA
____________
GAIM treats any observation having the value
99999 (5 nines) as missing (SORRY, blame Fortran). The procedure is
naturally geared to handle missing data; any fit (smooth, or linear) is
obtained via the backfitting algorithm, and is thus the result of a
simple smooth or linear fit. As such, the missing data is only excluded
when each individual function is estimated. i.e. Missing cases are NOT
removed. The contribution to the fit for a missing variable is the
average contribution for that variable. Cases with missing RESPONSES
are obviously removed.
MISSING DATA AND INTERACTIONS
_____________________________
For linear
terms, interactions between Linear and Dummy variables can be obtained
in the usual way--by multiplication. For non-parametric fits, this does
not work. Instead, one can use the missing value option to fudge it.
i.e.
________________________________________________________________
v1! 1, 2, 99999, 7, 6, 99999,999999,99999. !
v2! 99999, 99999, 9, 99999, 99999, 11, 8, 12 !
group! 1 1 2 1 1 2 2 2 !
_____!__________________________________________________________!
sets up two variables v1 and v2 according to wether group is 1 or 2.
pause
SPECIAL VARIABLES
_________________
Allowance is made for a WEIGHT
variable, and for the BINOMIAL model, one can specify N (r out of N).
For the COX model, one has to specify the censoring variable (0,1).
For the MATCHED CASE-CONTROL model, one codes the (SINGLE) case as 1,
the controls as 0, and then is asked for the matched set indicator. This is
an integer array with the numbers between 1 and K, where K is the number of sets.
SPECIAL OPTION
______________
This is used to get +-2 st deviation
curves for non-parametric fits and a covariance matrices for the
parameters of the linear fits. It uses a lot of CPU, since it
basically fits the current model N times. We are busy working on more
efficient approximate methods for calculating these quantities.
pause
SMOOTHER
________
The current version uses either a LOCAL LINEAR smoother,
or a CUBIC SPLINE Smoother.
The LOCAL LINEAR smoother has the following features:
Each neighbourhood contains 100s% of the total WEIGHT, with endpoint
modifications. Weight is defined initially, and also in the non-linear
models is defined implicitly. A preliminary scan thru the data
replaces all Y's tied in X by their average at that X. A post scan
ensures that all those fitted values are the same as well. See the
references for more details.
The CUBIC SPLINE smoother was kindly supplied by Finbarr O'Sullivan, UC
Berkeley. It uses the deBoor routines with some fancy Choleski work to
cross-validate in O(n) time. The present implementation requires a
fixed parameter, but this is easily changed. pause
USER SUPPLIED MODELS
____________________
The following are examples of
user supplied functions which operate in a similar way to the GLIM OWN
macros. They must be compiled and linked with the rest of the
program.
REAL FUNCTION DEVU(N,FITS,Y,W,NI) REAL FITS(1),Y(1),W(1),NI(1)
RSS=0.0 DO 1 I=1,N RSS=RSS+W(I)*(Y(I)-FITS(I))*2.0 1 CONTINUE
DEVU=RSS RETURN END
REAL FUNCTION LINKU(ETAHAT) LINKU=ETAHAT RETURN END
REAL FUNCTION LINVU(MUHAT) REAL MUHAT LINVU=MUHAT RETURN END
REAL FUNCTION DIRIVU(MUHAT) REAL MUHAT DIRIVU=1.0 RETURN END
REAL FUNCTION WEIGTU(W,MUHAT,NI) REAL W,MUHAT,NI WEIGTU=1.0
RETURN END
pause
PARTIAL RESIDUALS
Note: when partial residuals are requested, the GAIM prints out
x partial residual smooth(x) rank(x) fitted value y
pause
*************************************************************************
This Help file was created in a hurry, and will evolve in a natural
way.
*************************************************************************
stop
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine addtm
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
integer neffob(20)
common/neff/neffob
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
logical ok
ok=.false.
10011 continue
write(6,10020)
10020 format( 45h the following variables are not in the model)
icount=0
j=1
goto 10033
10031 j=j+1
10033 if((j).gt.(maxp))goto 10032
if(indat(j) .ne. 0)goto 10051
icount=icount+1
write(6,10060)j,(labels(k,j),k=1,10)
10060 format(i3, 5h ....,10a4)
10051 continue
goto 10031
10032 continue
if(icount .le. 0)goto 10081
10091 continue
write(6,10100)
10100 format( 33h index of new covariate (0= quit))
read(5,*)ir
if(ir .gt. 0)goto 10121
return
10121 continue
if(indat(ir) .ne. 0)goto 10141
goto 10092
10141 continue
write(6,10150)ir
10150 format(i3, 13h not allowed!)
goto 10091
10092 continue
p=p+1
indat(ir)=1
index(p)=ir
indmod(p)=1
modlab(p)=dumlab(1)
do 10161 i=1,n
x(i,p)=datam(i,ir)
10161 continue
10162 continue
cross(p)=.false.
write(6,10170)ir
10170 format( 19h span for covariate,i3, 18h in interval (0,1),/,
* 23h 0 --- cross-validated ,/, 13h 1 --- linear,/, 41h -
*--- cubic spline with lambda = - span ,/, 46h -999 cubic spl
*ine with auto span selection ,/, 32h k --- categorical with k
*levels)
read(5,*)spans(p)
if((spans(p) .ne. 0) .and. (spans(p) .ne. -999))goto 10191
cross(p)=.true.
indmod(p)=3
modlab(p)=dumlab(21)
10191 continue
if(spans(p) .le. 1.0)goto 10211
call getdum(ir)
goto 10221
10211 continue
if(spans(p) .ne. 1.0)goto 10241
indmod(p)=2
modlab(p)=dumlab(22)
10241 continue
call mysort(x(1,p),tag(1,p),1,n)
nef=n
10251 if(x(tag(nef,p),p) .ne. 99999) goto 10252
nef=nef-1
goto 10251
10252 continue
neffob(p)=nef
miss=n-nef
if(miss .le. 0)goto 10271
write(6,10280)ir,miss
10280 format (// 10h variable ,i2, 5h has ,i3, 21h missing observ
*ations,/, 56h they will receive a fitted value of 0 for this v
*ariable,/)
10271 continue
10221 continue
10201 continue
goto 10291
10081 continue
write(6,10300)
10300 format( 24h oops....they are all in)
return
10291 continue
10071 continue
goto 10011
10012 continue
return
end
subroutine getdum(ir)
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
integer neffob(20)
common/neff/neffob
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real catval(20),misans
call mysort(x(1,p),tag(1,p),1,n)
indcat=1
catval(1)=x(tag(1,p),p)
ii=2
goto 10313
10311 ii=ii+1
10313 if((ii).gt.(n))goto 10312
i=tag(ii,p)
if(w(i) .le. 0.0)goto 10331
if(x(i,p) .le. catval(indcat))goto 10351
indcat=indcat+1
if(indcat .le. 20)goto 10371
write(6,10380)
10380 format( 24h more than 20 categories)
p=p-1
indat(ir)=0
return
10371 continue
catval(indcat)=x(i,p)
10351 continue
10331 continue
goto 10311
10312 continue
icat=aint(spans(p))
ifmiss=0
if(catval(indcat) .ne. 99999)goto 10401
icat=icat+1
nef=n
do 10411 i=1,n
if(x(i,p) .ne. 99999)goto 10431
nef=nef-1
10431 continue
10411 continue
10412 continue
miss=n-nef
write(6,10440)ir,miss
10440 format (// 10h variable ,i2, 5h has ,i3, 21h missing observ
*ations,/, 45h do you want a separate category for them (y),/,
* 56h or do you want them to receive the average constant (n))
read(5,10450)misans
if(misans .ne. no)goto 10471
ifmiss=1
10471 continue
10401 continue
ndum=indcat-1
if(indcat .eq. icat)goto 10491
i1=indcat-ifmiss
i2=icat-ifmiss
write(6,10500)i1,i2
10500 format( 11h there are ,i3, 22h categories; you said ,i3,/,
* 44h do you wish to continue (yes) or go back to,
* 14h the menu (no))
read(5,10450)answer
10450 format(a1)
if(answer .ne. no)goto 10521
p=p-1
indat(ir)=0
return
goto 10531
10521 continue
if(indcat .le. icat)goto 10551
ndum=icat-1
write(6,10560)
10560 format( 39h there are more categories than dummies,/, 36h in
*itial categories will be combined)
10551 continue
10531 continue
10511 continue
10491 continue
i1=ndum-ifmiss
write(6,10570)i1,i1
10570 format( 10h creating ,i2, 17h dummy variables,/, 14h for
*the last ,i2, 11h categories)
indat(ir)=indcat+100
indc=indcat-ndum
j=p
goto 10583
10581 j=j+1
10583 if((j).gt.(p+ndum-1))goto 10582
if(j .le. 20)goto 10601
write(6,10610)
10610 format ( 35h oops, you`ve exceeded $p variables)
stop
10601 continue
indc=indc+1
cross(j)=.false.
index(j)=ir
indmod(j)=100+j
modlab(j)=dumlab(indc)
do 10621 i=1,n
x(i,j)=0.0
10621 continue
10622 continue
do 10631 i=1,n
if(datam(i,ir) .ne. catval(indc))goto 10651
x(i,j)=1.0
10651 continue
tag(i,j)=i
10631 continue
10632 continue
spans(j)=1.0
neffob(j)=n
goto 10581
10582 continue
p=p+ndum-1
if(catval(indcat) .ne. 99999 .or. ifmiss .ne. 0)goto 10671
modlab(p)=dumlab(23)
10671 continue
if(ifmiss .ne. 1)goto 10691
j=p-1
goto 10703
10701 j=j+(-1)
10703 if((-1)*((j)-(p+1-ndum)).gt.0)goto 10702
do 10711 i=1,n
x(i,j)=x(i,j)+99999*x(i,p)
10711 continue
10712 continue
call mysort(x(1,j),tag(1,j),1,n)
neffob(j)=nef
goto 10701
10702 continue
p=p-1
10691 continue
return
end
subroutine readd
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
character*60 dsname
common/dnm/dsname,ilen
logical ds
character*80 datan,datadn
character lett
write(6,10720)
10720 format( 29h Give the data set directory )
read(5,10730)dsname
10730 format(a60)
lett=dsname(1:1)
i=1
10741 if((lett.eq. 1h ) ) goto 10742
i=i+1
lett=dsname(i:i)
goto 10741
10742 continue
ilen=i
datan=dsname
datadn=dsname
datadn(i:i+6)= 6h/datad
datan(i:i+5)= 5h/data
inquire(file=datan, exist=ds)
if(.not.(.not.ds))goto 10761
write(6,10770)datan
10770 format(a80, 27h does not exist! exiting...)
stop
10761 continue
inquire(file=datadn, exist=ds)
if(.not.(.not.ds))goto 10791
write(6,10800)datadn
10800 format(a80, 27h does not exist! exiting...)
stop
10791 continue
open(1,file=datan)
rewind 1
open(11,file=datadn)
rewind 11
do 10811 i=1,20
indat(i)=-1
10811 continue
10812 continue
read(11,*)n,maxp
if(maxp .le. 20)goto 10831
write(6,10840)
10840 format( 19h too many variables)
stop
10831 continue
do 10851 i=1,maxp
read(11,10860)(labels(j,i),j=1,10)
indat(i)=0
10851 continue
10852 continue
10860 format(10a4)
i=1
goto 10873
10871 i=i+1
10873 if((i).gt.(n))goto 10872
read(1,*,end=10880)(datam(i,j),j=1,maxp)
goto 10871
10872 continue
return
10880 continue
write(6,10890)
10890 format( 23h you ve run out of data)
stop
end
subroutine deltm
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
integer neffob(20)
common/neff/neffob
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
logical ok
ok=.false.
10901 continue
if(p .ne. 0)goto 10921
write(6,10930)
10930 format( 11h null model)
return
10921 continue
write(6,10940)
10940 format( 41h the following variables are in the model)
j=1
goto 10953
10951 j=j+1
10953 if((j).gt.(p))goto 10952
write(6,10960)index(j),modlab(j),(labels(k,index(j)),k=1,10)
10960 format (i3, 5h ....,11a4)
goto 10951
10952 continue
write(6,10970)
10970 format( 43h index of covariate to be deleted (0= quit))
read(5,*)ir
if(ir .ne. 0)goto 10991
return
10991 continue
if((indat(ir) .ne. 1) .and. (indat(ir) .le. 100))goto 11011
indat(ir)=0
j=1
goto 11023
11021 j=j+1
11023 if((j).gt.(p))goto 11022
if(ir .ne. index(j))goto 11041
goto 11022
11041 continue
goto 11021
11022 continue
if(p .le. 1)goto 11061
ij=1
jj=j
11071 if(index(jj+1) .ne. ir .or. (jj+1) .gt. p) goto 11072
jj=jj+1
ij=ij+1
goto 11071
11072 continue
k=j
goto 11083
11081 k=k+1
11083 if((k).gt.(p-ij))goto 11082
index(k)=index(k+ij)
cross(k)=cross(k+ij)
spans(k)=spans(k+ij)
indmod(k)=indmod(k+ij)
modlab(k)=modlab(k+ij)
neffob(k)=neffob(k+ij)
i=1
goto 11093
11091 i=i+1
11093 if((i).gt.(n))goto 11092
tag(i,k)=tag(i,k+ij)
x(i,k)=x(i,k+ij)
goto 11091
11092 continue
goto 11081
11082 continue
p=p-ij
goto 11101
11061 continue
p=0
return
11101 continue
11051 continue
goto 11111
11011 continue
write(6,11120)ir
11120 format( 10h variable ,i3, 20h is not in the model)
11111 continue
11001 continue
goto 10901
10902 continue
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine bakfit(n,p,y,x,w,tag,spans,dofs,slopes,cross, etahat,s
*,s0,rss,thresh,nit,verbos,quiet,itmax)
integer n,p,tag(1000,1),nit,itmax
logical verbos,quiet,cross(1)
integer neffob(20)
common/neff/neffob
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
integer bacnit
common /bacn/bacnit
real y(1),x(1000,1),w(1),spans(1),dofs(1),slopes(1), etahat(1),s(
*1000,1),s0,rss,thresh
real z(1000),scrat(1000)
if(p .ne. 1)goto 10021
itmax=1
10021 continue
rss=devg(n,etahat,y,w,junk)
rssold=10*rss+10
s1=0.0
sn=0.0
se=0.0
i=1
goto 10033
10031 i=i+1
10033 if((i).gt.(n))goto 10032
se=se+w(i)*etahat(i)
s1=s1+w(i)*y(i)
sn=sn+w(i)
goto 10031
10032 continue
s1=s1/sn
se=se/sn
do 10041 i=1,n
etahat(i)=etahat(i)-se+s1
10041 continue
10042 continue
s0=s0-se+s1
nit=0
10051 if(abs((rssold-rss)/rssold) .le. thresh .or. nit .ge. itmax) goto
*10052
rssold=rss
j=1
goto 10063
10061 j=j+1
10063 if((j).gt.(p))goto 10062
i=1
goto 10073
10071 i=i+1
10073 if((i).gt.(n))goto 10072
etahat(i)=etahat(i)-s(i,j)
z(i)=y(i)-etahat(i)
goto 10071
10072 continue
call smooth(x(1,j),z,w,tag(1,j),spans(j),dofs(j),n,neffob(j), cro
*ss(j),s(1,j),s1,rss,quiet,scrat)
slopes(j)=scrat(1)
do 10081 i=1,n
etahat(i)=etahat(i)+s(i,j)+s1
10081 continue
10082 continue
s0=s0+s1
if(.not.(verbos))goto 10101
write(6,10110)j,rss
10110 format(10x, 26hinner loop (backfit) var= , i3, 5h rss=,f16.5
*)
10101 continue
goto 10061
10062 continue
nit=nit+1
if(.not.(.not.quiet))goto 10131
write(6,10140)nit,rss
10140 format(/10x, 26houter loop (backfit) nit= , i3, 5h rss=,f16.
*5,/)
10131 continue
goto 10051
10052 continue
bacnit=nit
return
end
subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv )
calls bsplvb
calculates value and deriv.s of all b-splines which do not vanish at x
c
c****** i n p u t ******
c t the knot array, of length left+k (at least)
c k the order of the b-splines to be evaluated
c x the point at which these values are sought
c left an integer indicating the left endpoint of the interval of
c interest. the k b-splines whose support contains the interval
c (t(left), t(left+1))
c are to be considered.
c a s s u m p t i o n - - - it is assumed that
c t(left) .lt. t(left+1)
c division by zero will result otherwise (in b s p l v b ).
c also, the output is as advertised only if
c t(left) .le. x .le. t(left+1) .
c nderiv an integer indicating that values of b-splines and their
c derivatives up to but not including the nderiv-th are asked
c for. ( nderiv is replaced internally by the integer in (1,k)
c closest to it.)
c
c****** w o r k a r e a ******
c a an array of order (k,k), to contain b-coeff.s of the derivat-
c ives of a certain order of the k b-splines of interest.
c
c****** o u t p u t ******
c dbiatx an array of order (k,nderiv). its entry (i,m) contains
c value of (m-1)st derivative of (left-k+i)-th b-spline of
c order k for knot sequence t , i=m,...,k; m=1,...,nderiv.
c
c****** m e t h o d ******
c values at x of all the relevant b-splines of order k,k-1,...,
c k+1-nderiv are generated via bsplvb and stored temporarily
c in dbiatx . then, the b-coeffs of the required derivatives of the
c b-splines of interest are generated by differencing, each from the
c preceding one of lower order, and combined with the values of b-
c splines of corresponding order in dbiatx to produce the desired
c values.
c
integer k,left,nderiv, i,ideriv,il,j,jlow,jp1mid,kp1,kp1mm
* ,ldummy,m,mhigh
real a(k,k),dbiatx(k,nderiv),t(1),x, factor,fkp1mm,sum
mhigh = max0(min0(nderiv,k),1)
c mhigh is usually equal to nderiv.
kp1 = k+1
call bsplvb(t,kp1-mhigh,1,x,left,dbiatx)
if (mhigh .eq. 1) go to 99
c the first column of dbiatx always contains the b-spline values
c for the current order. these are stored in column k+1-current
c order before bsplvb is called to put values for the next
c higher order on top of it.
ideriv = mhigh
do 15 m=2,mhigh
jp1mid = 1
do 11 j=ideriv,k
dbiatx(j,ideriv) = dbiatx(jp1mid,1)
11 jp1mid = jp1mid + 1
ideriv = ideriv - 1
call bsplvb(t,kp1-ideriv,2,x,left,dbiatx)
15 continue
c
c at this point, b(left-k+i, k+1-j)(x) is in dbiatx(i,j) for
c i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the
c first column of dbiatx is already in final form. to obtain cor-
c responding derivatives of b-splines in subsequent columns, gene-
c rate their b-repr. by differencing, then evaluate at x.
c
jlow = 1
do 20 i=1,k
do 19 j=jlow,k
19 a(j,i) = 0.
jlow = i
20 a(i,i) = 1.
c at this point, a(.,j) contains the b-coeffs for the j-th of the
c k b-splines of interest here.
c
do 40 m=2,mhigh
kp1mm = kp1 - m
fkp1mm = float(kp1mm)
il = left
i = k
c
c for j=1,...,k, construct b-coeffs of (m-1)st derivative of
c b-splines from those for preceding derivative by differencing
c and store again in a(.,j) . the fact that a(i,j) = 0 for
c i .lt. j is used.sed.
do 25 ldummy=1,kp1mm
factor = fkp1mm/(t(il+kp1mm) - t(il))
c the assumption that t(left).lt.t(left+1) makes denominator
c in factor nonzero.
do 24 j=1,i
24 a(i,j) = (a(i,j) - a(i-1,j))*factor
il = il - 1
25 i = i - 1
c
c for i=1,...,k, combine b-coeffs a(.,i) with b-spline values
c stored in dbiatx(.,m) to get value of (m-1)st derivative of
c i-th b-spline (of interest here) at x , and store in
c dbiatx(i,m). storage of this value over the value of a b-spline
c of order m there is safe since the remaining b-spline derivat-
c ive of the same order do not use this value due to the fact
c that a(j,i) = 0 for j .lt. i .
30 do 40 i=1,k
sum = 0.
jlow = max0(i,m)
do 35 j=jlow,k
35 sum = a(j,i)*dbiatx(j,m) + sum
40 dbiatx(i,m) = sum
99 return
end
subroutine bsplvb ( t, jhigh, index, x, left, biatx )
calculates the value of all possibly nonzero b-splines at x of order
c
c jout = max( jhigh , (j+1)*(index-1) )
c
c with knot sequence t .
c
c****** i n p u t ******
c t.....knot sequence, of length left + jout , assumed to be nonde-
c creasing. a s s u m p t i o n . . . .
c t(left) .lt. t(left + 1) .
c d i v i s i o n b y z e r o will result if t(left) = t(left+1)
c jhigh,
c index.....integers which determine the order jout = max(jhigh,
c (j+1)*(index-1)) of the b-splines whose values at x are to
c be returned. index is used to avoid recalculations when seve-
c ral columns of the triangular array of b-spline values are nee-
c ded (e.g., in bvalue or in bsplvd ). precisely,
c if index = 1 ,
c the calculation starts from scratch and the entire triangular
c array of b-spline values of orders 1,2,...,jhigh is generated
c order by order , i.e., column by column .
c if index = 2 ,
c only the b-spline values of order j+1, j+2, ..., jout are ge-
c nerated, the assumption being that biatx , j , deltal , deltar
c are, on entry, as they were on exit at the previous call.
c in particular, if jhigh = 0, then jout = j+1, i.e., just
c the next column of b-spline values is generated.
c
c w a r n i n g . . . the restriction jout .le. jmax (= 20) is im-
c posed arbitrarily by the dimension statement for deltal and
c deltar below, but is n o w h e r e c h e c k e d for .
c
c x.....the point at which the b-splines are to be evaluated.
c left.....an integer chosen (usually) so that
c t(left) .le. x .le. t(left+1) .
c
c****** o u t p u t ******
c biatx.....array of length jout , with biatx(i) containing the val-
c ue at x of the polynomial of order jout which agrees with
c the b-spline b(left-jout+i,jout,t) on the interval (t(left),
c t(left+1)) .
c
c****** m e t h o d ******
c the recurrence relation
c
c x - t(i) t(i+j+1) - x
c b(i,j+1)(x) = -----------b(i,j)(x) + ---------------b(i+1,j)(x)
c t(i+j)-t(i) t(i+j+1)-t(i+1)
c
c is used (repeatedly) to generate the (j+1)-vector b(left-j,j+1)(x),
c ...,b(left,j+1)(x) from the j-vector b(left-j+1,j)(x),...,
c b(left,j)(x), storing the new values in biatx over the old. the
c facts that
c b(i,1) = 1 if t(i) .le. x .lt. t(i+1)
c and that
c b(i,j)(x) = 0 unless t(i) .le. x .lt. t(i+j)
c are used. the particular organization of the calculations follows al-
c gorithm (8) in chapter x of the text.
c
parameter(jmax = 20)
integer index,jhigh,left, i,j,jp1
real biatx(jhigh),t(1),x, deltal(jmax),deltar(jmax),saved,term
c dimension biatx(jout), t(left+jout)
current fortran standard makes it impossible to specify the length of
c t and of biatx precisely without the introduction of otherwise
c superfluous additional arguments.
data j/1/
c save j,deltal,deltar (valid in fortran 77)
c
go to (10,20), index
10 j = 1
biatx(1) = 1.
if (j .ge. jhigh) go to 99
c
20 jp1 = j + 1
deltar(j) = t(left+j) - x
deltal(j) = x - t(left+1-j)
saved = 0.
do 26 i=1,j
term = biatx(i)/(deltar(i) + deltal(jp1-i))
biatx(i) = saved + deltar(i)*term
26 saved = deltal(jp1-i)*term
biatx(jp1) = saved
j = jp1
if (j .lt. jhigh) go to 20
c
99 return
end
function bvalue ( t, bcoef, n, k, x, jderiv )
calls interv
c
calculates value at x of jderiv-th derivative of spline from b-repr.
c the spline is taken to be continuous from the right.
c
c****** i n p u t ******
c t, bcoef, n, k......forms the b-representation of the spline f to
c be evaluated. specifically,
c t.....knot sequence, of length n+k, assumed nondecreasing.
c bcoef.....b-coefficient sequence, of length n .
c n.....length of bcoef and dimension of s(k,t),
c a s s u m e d positive .
c k.....order of the spline .
c
c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed
c arbitrarily by the dimension statement for aj, dm, dm below,
c but is n o w h e r e c h e c k e d for.
c
c x.....the point at which to evaluate .
c jderiv.....integer giving the order of the derivative to be evaluated
c a s s u m e d to be zero or positive.
c
c****** o u t p u t ******
c bvalue.....the value of the (jderiv)-th derivative of f at x .
c
c****** m e t h o d ******
c the nontrivial knot interval (t(i),t(i+1)) containing x is lo-
c cated with the aid of interv . the k b-coeffs of f relevant for
c this interval are then obtained from bcoef (or taken to be zero if
c not explicitly available) and are then differenced jderiv times to
c obtain the b-coeffs of (d**jderiv)f relevant for that interval.
c precisely, with j = jderiv, we have from x.(12) of the text that
c
c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) )
c
c where
c / bcoef(.), , j .eq. 0
c /
c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1)
c / ----------------------------- , j .gt. 0
c / (t(.+k-j) - t(.))/(k-j)
c
c then, we use repeatedly the fact that
c
c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) )
c with
c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
c a(.,x) = ---------------------------------------
c (x - t(.)) + (t(.+m-1) - x)
c
c to write (d**j)f(x) eventually as a linear combination of b-splines
c of order 1 , and the coefficient for b(i,1,t)(x) must then
c be the desired number (d**j)f(x). (see x.(17)-(19) of text).
c
parameter(kmax = 20)
integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,km1,mflag,nmi
real bcoef(n),t(1),x, aj(kmax),dm(kmax),dp(kmax),fkmj
c dimension t(n+k)
current fortran standard makes it impossible to specify the length of
c t precisely without the introduction of otherwise superfluous
c additional arguments.
bvalue = 0.
if (jderiv .ge. k) go to 99
c
c *** find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and
c t(i) .le. x .lt. t(i+1) . if no such i can be found, x lies
c outside the support of the spline f and bvalue = 0.
c (the asymmetry in this choice of i makes f rightcontinuous)
if( (x.ne.t(n+1)) .or. (t(n+1).ne.t(n+k)) ) go to 700
i = n
go to 701
700 call interv ( t, n+k, x, i, mflag )
if (mflag .ne. 0) go to 99
701 continue
c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
km1 = k - 1
if (km1 .gt. 0) go to 1
bvalue = bcoef(i)
go to 99
c
c *** store the k b-spline coefficients relevant for the knot interval
c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dm(j) = x - t(i+1-j),
c dp(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
c from input to zero. set any t.s not obtainable equal to t(1) or
c to t(n+k) appropriately.
1 jcmin = 1
imk = i - k
if (imk .ge. 0) go to 8
jcmin = 1 - imk
do 5 j=1,i
5 dm(j) = x - t(i+1-j)
do 6 j=i,km1
aj(k-j) = 0.
6 dm(j) = dm(i)
go to 10
8 do 9 j=1,km1
9 dm(j) = x - t(i+1-j)
c
10 jcmax = k
nmi = n - i
if (nmi .ge. 0) go to 18
jcmax = k + nmi
do 15 j=1,jcmax
15 dp(j) = t(i+j) - x
do 16 j=jcmax,km1
aj(j+1) = 0.
16 dp(j) = dp(jcmax)
go to 20
18 do 19 j=1,km1
19 dp(j) = t(i+j) - x
c
20 do 21 jc=jcmin,jcmax
21 aj(jc) = bcoef(imk + jc)
c
c *** difference the coefficients jderiv times.
if (jderiv .eq. 0) go to 30
do 23 j=1,jderiv
kmj = k-j
fkmj = float(kmj)
ilo = kmj
do 23 jj=1,kmj
aj(jj) = ((aj(jj+1) - aj(jj))/(dm(ilo) + dp(jj)))*fkmj
23 ilo = ilo - 1
c
c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative,
c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
30 if (jderiv .eq. km1) go to 39
jdrvp1 = jderiv + 1
do 33 j=jdrvp1,km1
kmj = k-j
ilo = kmj
do 33 jj=1,kmj
aj(jj) = (aj(jj+1)*dm(ilo) + aj(jj)*dp(jj))/(dm(ilo)+dp(jj))
33 ilo = ilo - 1
39 bvalue = aj(1)
c
99 return
end
subroutine interv ( xt, lxt, x, left, mflag )
computes left = max( i ; 1 .le. i .le. lxt .and. xt(i) .le. x ) .
c
c****** i n p u t ******
c xt.....a real sequence, of length lxt , assumed to be nondecreasing
c lxt.....number of terms in the sequence xt .
c x.....the point whose location with respect to the sequence xt is
c to be determined.
c
c****** o u t p u t ******
c left, mflag.....both integers, whose value is
c
c 1 -1 if x .lt. xt(1)
c i 0 if xt(i) .le. x .lt. xt(i+1)
c lxt 1 if xt(lxt) .le. x
c
c in particular, mflag = 0 is the 'usual' case. mflag .ne. 0
c indicates that x lies outside the halfopen interval
c xt(1) .le. y .lt. xt(lxt) . the asymmetric treatment of the
c interval is due to the decision to make all pp functions cont-
c inuous from the right.
c
c****** m e t h o d ******
c the program is designed to be efficient in the common situation that
c it is called repeatedly, with x taken from an increasing or decrea-
c sing sequence. this will happen, e.g., when a pp function is to be
c graphed. the first guess for left is therefore taken to be the val-
c ue returned at the previous call and stored in the l o c a l varia-
c ble ilo . a first check ascertains that ilo .lt. lxt (this is nec-
c essary since the present call may have nothing to do with the previ-
c ous call). then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left =
c ilo and are done after just three comparisons.
c otherwise, we repeatedly double the difference istep = ihi - ilo
c while also moving ilo and ihi in the direction of x , until
c xt(ilo) .le. x .lt. xt(ihi) ,
c after which we use bisection to get, in addition, ilo+1 = ihi .
c left = ilo is then returned.
c
integer left,lxt,mflag, ihi,ilo,istep,middle
real x,xt(lxt)
data ilo /1/
c save ilo (a valid fortran statement in the new 1977 standard)
ihi = ilo + 1
if (ihi .lt. lxt) go to 20
if (x .ge. xt(lxt)) go to 110
if (lxt .le. 1) go to 90
ilo = lxt - 1
ihi = lxt
c
20 if (x .ge. xt(ihi)) go to 40
if (x .ge. xt(ilo)) go to 100
c
c **** now x .lt. xt(ilo) . decrease ilo to capture x .
30 istep = 1
31 ihi = ilo
ilo = ihi - istep
if (ilo .le. 1) go to 35
if (x .ge. xt(ilo)) go to 50
istep = istep*2
go to 31
35 ilo = 1
if (x .lt. xt(1)) go to 90
go to 50
c **** now x .ge. xt(ihi) . increase ihi to capture x .
40 istep = 1
41 ilo = ihi
ihi = ilo + istep
if (ihi .ge. lxt) go to 45
if (x .lt. xt(ihi)) go to 50
istep = istep*2
go to 41
45 if (x .ge. xt(lxt)) go to 110
ihi = lxt
c
c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
50 middle = (ilo + ihi)/2
if (middle .eq. ilo) go to 100
c note. it is assumed that middle = ilo in case ihi = ilo+1 .
if (x .lt. xt(middle)) go to 53
ilo = middle
go to 50
53 ihi = middle
go to 50
c**** set output and return.
90 mflag = -1
left = 1
return
100 mflag = 0
left = ilo
return
110 mflag = 1
left = lxt
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine dumpm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real pim(15),pm(15)
real ytot(15)
write(6,10010)(thetak(j),j=1,mcat)
10010 format(5f10.4)
i=1
goto 10023
10021 i=i+1
10023 if((i).gt.(n))goto 10022
call linkm(mcat,thetak,etahat(i),pim,pm)
ytot(1)=ym(i,1)
j=2
goto 10033
10031 j=j+1
10033 if((j).gt.(mcat-1))goto 10032
ytot(j)=ytot(j-1)+ym(i,j)
goto 10031
10032 continue
mcatm1=mcat-1
write(4,10040)i,(pim(j),j=1,mcatm1)
write(4,10040)i,(ytot(j),j=1,mcatm1)
goto 10021
10022 continue
10040 format(i5,10f7.4)
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine dumpp(n,p,x,lambda)
integer n,p
real x(506,p),lambda(506,2)
real y(6)
real sgrid(6,15,15),lam1(15),lam2(15)
integer table(15,15),lrlam1(2,15), lrlam2(2,15)
common /cont/ sgrid,lam1,lam2,table,lrlam1,lrlam2
igrids=15
i=p+2
write(14,*)n,i,igrids
do 10011 i=1,n
write(14,10020)(lambda(i,j),j=1,2),(x(i,j),j=1,p)
10011 continue
10012 continue
10020 format(8f10.4)
i=1
goto 10033
10031 i=i+1
10033 if((i).gt.(igrids))goto 10032
jgrid=lrlam1(2,i)-lrlam1(1,i)+1
write(14,*)jgrid
if(lrlam1(1,i) .le. lrlam1(2,i))goto 10051
goto 10031
10051 continue
j=lrlam1(1,i)
goto 10063
10061 j=j+1
10063 if((j).gt.(lrlam1(2,i)))goto 10062
k=1
goto 10073
10071 k=k+1
10073 if((k).gt.(p))goto 10072
y(k)=sgrid(k,i,j)
goto 10071
10072 continue
write(14,10020)lam1(i),lam2(j),(y(ij),ij=1,p)
goto 10061
10062 continue
goto 10031
10032 continue
j=1
goto 10083
10081 j=j+1
10083 if((j).gt.(igrids))goto 10082
jgrid=lrlam2(2,j)-lrlam2(1,j)+1
write(14,*)jgrid
if(lrlam2(1,j) .le. lrlam2(2,j))goto 10101
goto 10081
10101 continue
i=lrlam2(1,j)
goto 10113
10111 i=i+1
10113 if((i).gt.(lrlam2(2,j)))goto 10112
k=1
goto 10123
10121 k=k+1
10123 if((k).gt.(p))goto 10122
y(k)=sgrid(k,i,j)
goto 10121
10122 continue
write(14,10020)lam1(i),lam2(j),(y(ij),ij=1,p)
goto 10111
10112 continue
goto 10081
10082 continue
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine dumps(diriv)
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
integer neffob(20)
common/neff/neffob
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
character*60 dsname
common/dnm/dsname,ilen
real diriv
external diriv
character*90 plotn
real sep,part
write(6,10010)
10010 format( 36h want to dump a plot of this model? )
read(5,10020)ans
if(ans .ne. no)goto 10041
return
10041 continue
10020 format(a1)
nitav=numsum/nitout
r2=(devnll-devian)*100.0/devnll
if(p .ne. 1 .or. mod .eq. mult)goto 10061
call nplot(plotn)
oldx=-9999
oldy=-9999
i=1
goto 10073
10071 i=i+1
10073 if((i).gt.(neffob(1)))goto 10072
ii=tag(i,1)
if((x(ii,1) .eq. oldx) .and. (y(ii) .eq. oldy))goto 10091
write(2,10100)x(ii,1),y(ii),muhat(ii)
oldx=x(ii,1)
oldy=y(ii)
10091 continue
goto 10071
10072 continue
close(2)
goto 10111
10061 continue
write(6,10120)
10120 format( 32h want to dump plots separately? )
read(5,10020)sep
write(6,10130)
10130 format( 25h want partial residuals? )
read(5,10020)part
j=1
goto 10143
10141 j=j+1
10143 if((j).gt.(p))goto 10142
write(6,10150)j,(labels(k,index(j)),k=1,10)
10150 format( 10h Variable ,i3, 2h: ,10a4)
if(sep .ne. yes)goto 10171
write(6,10180)
10180 format( 18h do you want it...)
read(5,10020)ans
if(ans .ne. no)goto 10201
goto 10141
10201 continue
10171 continue
call nplot(plotn)
totsm=0.0
i=1
goto 10213
10211 i=i+1
10213 if((i).gt.(neffob(j)))goto 10212
totsm=totsm+s(tag(i,j),j)
goto 10211
10212 continue
totsm=totsm/neffob(j)
if(part .ne. yes)goto 10231
i=1
goto 10243
10241 i=i+1
10243 if((i).gt.(neffob(j)))goto 10242
ii=tag(i,j)
sput=s(ii,j)
z=s(ii,j)+(y(ii)-ni(ii)*muhat(ii))*diriv(muhat(ii))
if(mod .ne. mult)goto 10261
z=-z
sput=-sput
10261 continue
write(2,10100)x(ii,j),z,sput,ii,muhat(ii),y(ii)
goto 10241
10242 continue
close(2)
goto 10271
10231 continue
oldx=-9999
i=1
goto 10283
10281 i=i+1
10283 if((i).gt.(neffob(j)))goto 10282
ii=tag(i,j)
if(x(ii,j) .eq. oldx)goto 10301
sput=s(ii,j)
if(mod .ne. mult)goto 10321
sput=-sput
10321 continue
write(2,10100)x(ii,j),sput
oldx=x(ii,j)
10301 continue
goto 10281
10282 continue
close(2)
10271 continue
10221 continue
write(6,10330)totsm
10330 format( 24h average of function is ,f10.3)
goto 10141
10142 continue
10111 continue
10051 continue
if(.not.(vlink))goto 10351
write(6,10360)
10360 format( 14h variable link)
call nplot(plotn)
iii=0
i=1
goto 10373
10371 i=i+1
10373 if((i).gt.(n))goto 10372
write(2,10100)argf(i),slink(i)
goto 10371
10372 continue
close(2)
10351 continue
10380 format( 24h data a; input x y yhat;)
10390 format( 7h data a,i1, 15h; input x yhat;)
10400 format( 39h title .c=blue estimated mean function; /, 37hproc
* gplot; plot y*x yhat*x /overlay;,/, 54hsymbol1 v=star c=red;s
*ymbol2 i=join v=diamond c=green;)
10410 format( 33h label x = eta yhat=f(eta);cards;)
10420 format( 30h title .c=blue estimated link; /, 26hproc gplot; p
*lot yhat*x ;,/, 32hsymbol i=join c=green v=diamond;)
10430 format( 34h title .c=blue estimated function; /, 26hproc gplo
*t; plot yhat*x ;,/, 32hsymbol i=join c=green v=diamond;)
10440 format( 56h title .c=blue estimated function and partial residua
*ls; /, 37hproc gplot; plot y*x yhat*x /overlay;,/, 45hsymbo
*l1 v=star c=red;symbol2 i=join c=green;)
10450 format( 8hlabel x=,10a4, 17hyhat=muhat;cards;)
10460 format( 23hfootnote2 .c=blue span=,f4.2, 5h dof=,f6.2, 8h
*n-nmiss=,i3, 13h #iterations=, i3, 1h(,i3, 2h);,/, 27
*hfootnote3 .c=blue deviance=, f8.2, 4h r2=,f8.2, 2h%;)
10470 format( 8hlabel x=,10a4, 13hyhat=s;cards;)
10480 format(3f10.4)
10100 format(3f15.4,i5,2f10.4)
return
end
subroutine nplot(plotn)
character*60 dsname
common/dnm/dsname,ilen
character*90 plotn
character*30 plf
do 10491 i=1,90
plotn(i:i)= 1h
10491 continue
10492 continue
do 10501 i=1,30
plf(i:i)= 1h
10501 continue
10502 continue
close(2)
write(6,10510)
10510 format( 30h Give a name for the plot file)
read(5,10520)plf
10520 format(a30)
plotn=dsname
plotn(ilen:ilen)= 1h/
plotn(ilen+1:ilen+30)=plf(1:30)
open(2,file=plotn)
rewind 2
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
block data
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real match
integer indms,kms(1000)
common/mtchco/match,indms,kms
data yes/ 1hy/, no/ 1hn/, binom/ 1hb/,gauss/ 1hg/,cox/
* 1hc/
data mult/ 1hm/
data match/ 1hr/
data user/ 1hu/
data n/1000/, p/20/
data verbos/.false./, quiet/.false./
data vlink/.false./
data score/ 1hs/,local/ 1hl/
data thresh/0.001/,thbak/0.001/
data dumlab/ 4h , 4hd2: , 4hd3: , 4hd4: , 4hd5:
* , 4hd6: , 4hd7: , 4hd8: , 4hd9: , 4hd10:, 4hd11
*:, 4hd12:, 4hd13:, 4hd14:, 4hd15:, 4hd16:, 4hd
*17:, 4hd18:, 4hd19:, 4hd20:, 4hcv: , 4hlin:, 4hm
*i: ,7* 4h /
end
external devb,devg,linkg,linkb,dirivg,dirivb,weigtg,weigtb, devu,d
*evc,linku,dirivu,weigtu,linvg,linvb,linvu,linkmc,linvmc,devmc
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
integer neffob(20)
common/neff/neffob
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real z(1000),wz(1000)
common /wvc/z,wz
character*60 dsname
common/dnm/dsname,ilen
real match
integer indms,kms(1000)
common/mtchco/match,indms,kms
real adterm,determ,moterm,runmod,specia,quit
data adterm/ 1ha/, determ/ 1hd/, moterm/ 1hm/,specia/
*1hs/, quit/ 1hq/
data help/ 1hh/
data runmod/ 1hr/
call title
call readd
call name7
answer=yes
10011 if(answer .ne. yes) goto 10012
call model
10021 if(answer .eq. quit) goto 10022
10031 continue
prevan=answer
write(6,10040)
10040 format( 50h add, delete, modify, run, specials, help or quit)
read(5,10050)answer
if(answer .ne. adterm)goto 10071
call addtm
goto10061
10071 if(answer .ne. determ)goto 10081
call deltm
goto10061
10081 if(answer .ne. help)goto 10091
call helpg
goto10061
10091 if(answer .ne. moterm)goto 10101
call modtm
goto10061
10101 if(answer .ne. specia)goto 10111
if(prevan .eq. runmod)goto 10131
write(6,10140)
10140 format( 20h run the model first)
goto 10151
10131 continue
if(mod .ne. gauss)goto 10171
call hat(weigtg)
goto10161
10171 if(mod .ne. binom)goto 10181
call hat(weigtb)
goto10161
10181 if(mod .ne. cox)goto 10191
call hat(weigtg)
goto 10201
10191 continue
write(6,10210)mod
10210 format( 27h sorry, not yet implemented, 16hfor model type ,
*a1)
10201 continue
10161 continue
10151 continue
10121 continue
goto10061
10111 if((answer .ne. quit) .and. (answer .ne. runmod))goto 10221
goto 10032
10221 continue
10061 continue
goto 10031
10032 continue
if(answer .ne. quit)goto 10241
goto 10022
10241 continue
if(mod .ne. gauss)goto 10261
call gam(linkg,dirivg,weigtg,devg,linvg)
call info
call dumps(dirivg)
goto10251
10261 if(mod .ne. binom)goto 10271
call gam(linkb,dirivb,weigtb,devb,linvb)
call info
call dumps(dirivb)
goto10251
10271 if(mod .ne. cox)goto 10281
nq=n
call gam(linkg,dirivg,weigtg,devc,linvg)
call info
call dumps(dirivg)
goto10251
10281 if(mod .ne. user)goto 10291
call gam(linku,dirivu,weigtu,devu,linvu)
call info
call dumps(dirivu)
goto10251
10291 if(mod .ne. mult)goto 10301
call multam
call minfo
call dumps(dirivb)
write(6,10310)
10310 format( 38h want fitted cumulative probabilities )
read(5,10050)ans
if(ans .ne. yes)goto 10331
call dumpm
10331 continue
goto10251
10301 if(mod .ne. match)goto 10341
call gam(linkmc,dirivb,weigtb,devmc,linvmc)
call info
call dumps(dirivb)
goto 10351
10341 continue
write(6,10360)mod
10360 format( 32h sorry, we dont have model type ,a1)
10351 continue
10251 continue
goto 10021
10022 continue
write(6,10370)
10370 format( 15h another model?)
read(5,10050)answer
10050 format(a1)
goto 10011
10012 continue
stop
end
subroutine model
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real match
integer indms,kms(1000)
common/mtchco/match,indms,kms
real v(20),more,quit
logical ok
data quit/ 1hq/
ok=.false.
vlink=.false.
10381 if(ok) goto 10382
write(6,10390)
10390 format( 48h which model...binomial(b), gaussian(g), cox(c),,/,
* 17h multinomial(m), , 44h matched case control(r), user(u)
*or quit(q))
read(5,10400)mod
10400 format(a1)
if((mod .ne. gauss) .and. ((mod .ne. binom) .and. ((mod .ne. cox)
*.and. ((mod .ne. user) .and. ((mod .ne. mult) .and. (mod .ne. matc
*h))))))goto 10421
ok=.true.
goto10411
10421 if(mod .ne. quit)goto 10431
stop
goto 10441
10431 continue
write(6,10450)mod
10450 format( 32h sorry, we dont have model type ,a1)
10441 continue
10411 continue
goto 10381
10382 continue
esttyp=score
if(mod .ne. mult)goto 10471
call multin
return
10471 continue
j=1
goto 10483
10481 j=j+1
10483 if((j).gt.(maxp))goto 10482
indat(j)=0
write(6,10490)j,(labels(k,j),k=1,10)
10490 format(i3, 5h ....,10a4)
goto 10481
10482 continue
ir=maxp+1
10501 if((ir .le. maxp) .and. (ir .gt. 0)) goto 10502
write(6,10510)
10510 format( 19h response variable?)
read(5,*)ir
goto 10501
10502 continue
indat(ir)=2
indy=ir
do 10521 i=1,n
y(i)=datam(i,ir)
10521 continue
10522 continue
ir=maxp+1
10531 if((ir .le. maxp) .and. (ir .ne. indy)) goto 10532
write(6,10540)
10540 format( 55h give index of weight variable 0 is no weight variab
*le)
read(5,*)ir
goto 10531
10532 continue
indw=ir
if(ir .gt. 0)goto 10561
do 10571 i=1,n
w(i)=1.0
10571 continue
10572 continue
goto 10581
10561 continue
totw=0.0
do 10591 i=1,n
w(i)=abs(datam(i,ir))
totw=totw+w(i)
10591 continue
10592 continue
do 10601 i=1,n
w(i)=w(i)*n/totw
10601 continue
10602 continue
indat(ir)=3
10581 continue
10551 continue
nmiss=0
do 10611 i=1,n
if(y(i) .ne. 99999)goto 10631
w(i)=0.0
y(i)=0.0
nmiss=nmiss+1
10631 continue
10611 continue
10612 continue
if(nmiss .le. 0)goto 10651
ratmis=(100*nmiss)/n
write(6,10660)nmiss,ratmis
10660 format (i4, 34h responses missing (=> 99999), or ,f4.1, 14h% o
*f the data;,/, 31h their weights will be set to 0)
10651 continue
p=0
do 10671 i=1,n
ni(i)=1.0
10671 continue
10672 continue
if(mod .ne. binom)goto 10691
write(6,10700)
10700 format( 44h y out of n ... give index of n 0 makes n=1)
ir=maxp+1
10711 if((ir .le. maxp) .and. (ir .ne. indy)) goto 10712
read(5,*)ir
goto 10711
10712 continue
if(ir .gt. 0)goto 10731
indni=0
goto 10741
10731 continue
indni=ir
do 10751 i=1,n
ni(i)=datam(i,ir)
y(i)=y(i)/ni(i)
10751 continue
10752 continue
indat(ir)=4
10741 continue
10721 continue
10691 continue
if(mod .ne. match)goto 10771
write(6,10780)
10780 format( 36h give index of matched set indicator)
ir=maxp+1
10791 if((ir .le. maxp) .and. (ir .ne. indy)) goto 10792
read(5,*)ir
goto 10791
10792 continue
if(ir .le. 0)goto 10811
indms=ir
do 10821 i=1,n
kms(i)=datam(i,ir)
10821 continue
10822 continue
indat(ir)=6
10811 continue
10771 continue
if(mod .ne. cox)goto 10841
write(6,10850)
10850 format( 33h give index of censoring variable)
read(5,*) ic
indat(ic)=5
do 10861 i=1,n
icensq(i)=datam(i,ic)
10861 continue
10862 continue
write(6,10870)
10870 format( 55h there should be no censored obs before the 1st failu
*re)
do 10881 i=1,n
tq(i)=y(i)
10881 continue
10882 continue
call mysort(tq,tagy,1,n)
do 10891 i=1,n
tq(i)=y(tagy(i))
aq(i)=-1.0*icensq(tagy(i))
10891 continue
10892 continue
j=1
10900 continue
jj=j
10911 if(tq(jj+1) .ne. tq(j)) goto 10912
jj=jj+1
goto 10911
10912 continue
if(jj.gt.j) call psort(aq,tagy,j,jj)
j=jj+1
if(j.lt.n) go to 10900
do 10921 i=1,n
tq(i)=y(i)
aq(i)=icensq(i)
bq(i)=w(i)
do 10931 j=1,maxp
xt(i,j)=datam(i,j)
10931 continue
10932 continue
10921 continue
10922 continue
do 10941 i=1,n
y(i)=tq(tagy(i))
icensq(i)=aq(tagy(i))
w(i)=bq(tagy(i))
do 10951 j=1,maxp
datam(i,j)=xt(tagy(i),j)
10951 continue
10952 continue
10941 continue
10942 continue
10841 continue
return
end
subroutine modtm
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
10400 format(a1)
logical ok
vlink=.false.
write(6,10960)
10960 format( 45h do you want to change the thresh-holds (y/n))
read(5,10400)answer
if(answer .ne. yes)goto 10981
write(6,10990)thresh
10990 format( 23h change gam threshold (,f10.8, 6h y/n ))
read(5,10400)answer
if(answer .ne. yes)goto 11011
write(6,11020)
11020 format( 26h give new gam thresh-hold )
read(5,*)thresh
11011 continue
write(6,11030)thbak
11030 format( 28h change back-fit threshold (,f10.8, 6h y/n ))
read(5,10400)answer
if(answer .ne. yes)goto 11051
write(6,11060)
11060 format( 31h give new back-fit thresh-hold )
read(5,*)thbak
11051 continue
10981 continue
ok=.false.
11071 continue
write(6,11080)
11080 format( 40h the following variables may be modified)
j=1
goto 11093
11091 j=j+1
11093 if((j).gt.(p))goto 11092
if(.not.(cross(j)))goto 11111
spans(j)=0.0
11111 continue
if(indmod(j) .gt. 3)goto 11131
write(6,11140)index(j),spans(j),modlab(j),(labels(k,index(j)),k=1,
*10)
11140 format (i3, 10h ....span=,f5.2, 6h .... ,11a4)
11131 continue
goto 11091
11092 continue
write(6,11150)
11150 format( 44h index of covariate to be modified (0= quit))
read(5,*)ir
if(ir .ne. 0)goto 11171
return
11171 continue
if(indat(ir) .eq. 1)goto 11191
write(6,11200)
11200 format( 21h naughty, naughty....)
goto 11211
11191 continue
do 11221 j=1,p
if(index(j) .ne. ir)goto 11241
goto 11222
11241 continue
11221 continue
11222 continue
write(6,11250)ir,spans(j)
11250 format ( 10h variable i3, 10h has span=,f5.2,/, 62h give n
*ew span in interval (-?,1) (1 => linear, - => spline))
modlab(j)=dumlab(1)
indmod(j)=1
read(5,*)spans(j)
cross(j)=.false.
if(spans(j) .ne. 1.0)goto 11271
indmod(j)=2
modlab(j)=dumlab(22)
11271 continue
if((spans(j) .ne. 0.0) .and. (spans(j) .ne. -999))goto 11291
cross(j)=.true.
indmod(j)=3
modlab(j)=dumlab(21)
11291 continue
11211 continue
11181 continue
goto 11071
11072 continue
return
end
subroutine name7
character*60 dsname
common/dnm/dsname,ilen
character*90 devn
character*6 devn2
close(7)
devn=dsname
devn(ilen:ilen)= 1h/
devn(ilen+1:ilen+6)= 6hanodev
open(7,file=devn)
rewind 7
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine gam(link,diriv,weight,dev,linv)
external link,diriv,weight,dev,linv
real linv,diriv,weight,dev,scrat1(1000),scrat2(1000)
real scrat3(1000),scrat4(1000,3)
integer q
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real z(1000),wz(1000)
common /wvc/z,wz
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
itss=30
call init(linv)
if(mod .ne. cox)goto 10021
call initcx
10021 continue
devnew=dev(n,muhat,y,w,ni)
devnll=devnew
devold=10*devnew+10
if(mod .ne. cox)goto 10041
call stpar
10041 continue
nit=0
numsm=0
if(.not.(.not.vlink))goto 10061
10071 if(abs((devold-devnew)/devold) .le. thresh .or. nit .ge. maxit) go
*to 10072
nit=nit+1
if(mod .eq. cox)goto 10091
i=1
goto 10103
10101 i=i+1
10103 if((i).gt.(n))goto 10102
z(i)=etahat(i)+(y(i)-muhat(i))*diriv(muhat(i))
wz(i)=weight(w(i),muhat(i),ni(i))
goto 10101
10102 continue
goto 10111
10091 continue
call calcz
10111 continue
10081 continue
call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,etahat, s,s0,r
*ss,thbak,nitb,verbos,quiet,itss)
numsm=numsm+nitb
call link(n,etahat,muhat)
devold=devnew
devnew=dev(n,muhat,y,w,ni)
if(.not.(.not.quiet))goto 10131
write(6,10140)nit,devnew
10140 format(/10x, 22houter loop (gam) nit= , i3, 5h dev=,f16.5,//
*)
10131 continue
goto 10071
10072 continue
10061 continue
if(.not.(vlink))goto 10161
10171 if(abs((devold-devnew)/devold) .le. thresh .or. nit .ge. 5) goto 1
*0172
devold=devnew
devnei=devold
devoli=devnei*10
nitin=0
10181 if(abs((devoli-devnei)/devoli) .le. thresh .or. nitin .ge. maxit)
*goto 10182
nitin=nitin+1
i=1
goto 10193
10191 i=i+1
10193 if((i).gt.(n))goto 10192
fdirin=fpinv(etahat(i))
z(i)=etahat(i)+ (y(i)-muhat(i))*diriv(muhat(i))*fdirin
wz(i)=weight(w(i),muhat(i),ni(i))/(fdirin*fdirin)
goto 10191
10192 continue
call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,etahat, s,s0,r
*ss,thbak,nitb,verbos,quiet,itss)
numsm=numsm+nitb
i=1
goto 10203
10201 i=i+1
10203 if((i).gt.(n))goto 10202
sleta(i)=flink(etahat(i))
goto 10201
10202 continue
call link(n,etahat,muhat)
devoli=devnei
devnei=dev(n,muhat,y,w,ni)
if(.not.(.not.quiet))goto 10221
write(6,10230)nitin,devnei
10230 format(/10x, 27hinner loop eta (gam) nit= , i3, 5h dev=,f16
*.5,//)
10221 continue
goto 10181
10182 continue
if(.not.(.not.vflag))goto 10251
slinc=devnei
10251 continue
devnes=devnei
devols=devnes*100
nitsin=0
vflag=.true.
call mysort(etahat,sltag,1,n)
do 10261 i=1,n
argf(i)=etahat(sltag(i))
10261 continue
10262 continue
10271 if(abs((devols-devnes)/devols) .le. thresh .or. nitsin .ge. 6) got
*o 10272
nitsin=nitsin+1
i=1
goto 10283
10281 i=i+1
10283 if((i).gt.(n))goto 10282
z(i)=sleta(i) +(y(i)-muhat(i))*diriv(muhat(i))
wz(i)=weight(w(i),muhat(i),ni(i))
goto 10281
10282 continue
call smooth(etahat,z,wz,sltag,0.5,sldof,n,n, .false.,sleta,sv,rss
*,quiet,scrat1)
do 10291 i=1,n
slink(i)=sleta(sltag(i))+sv
10291 continue
10292 continue
call montne(slink,n)
call link(n,slink,scrat1)
i=1
goto 10303
10301 i=i+1
10303 if((i).gt.(n))goto 10302
sleta(sltag(i))=slink(i)
muhat(sltag(i))=scrat1(i)
goto 10301
10302 continue
devols=devnes
devnes=dev(n,muhat,y,w,ni)
if(.not.(.not.quiet))goto 10321
write(6,10330)nitsin,devnes
10330 format(/10x, 30hinner loop s(eta) (gam) nit= , i3, 5h dev=,
*f16.5,//)
10321 continue
goto 10271
10272 continue
call der(n,argf,slink,w,0.001,sldir,scrat4)
dif=slink(n)-slink(1)
idif=1
if(dif .ge. 0)goto 10351
idif=-1
10351 continue
do 10361 i=1,n
if(sldir(i)*idif .ge. 0.01)goto 10381
sldir(i)=idif*0.01
10381 continue
10361 continue
10362 continue
devold=devnew
devnew=devnes
nit=nit+1
if(.not.(.not.quiet))goto 10401
write(6,10410)nit,devnew
10410 format(10x,30( 1h_),/10x, 22houter loop (gam) nit= , i3,
*5h dev=,f16.5,/10x,30( 1h-)/)
10401 continue
goto 10171
10172 continue
10161 continue
devian=devnew
nitout=nit
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine hat(weight)
external weight
real weight
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
integer neffob(20)
common/neff/neffob
real z(1000),wz(1000)
common /wvc/z,wz
integer bacnit
common /bacn/bacnit
character*60 dsname
common/dnm/dsname,ilen
logical parm1,parm2
real zero
integer its
real*8 vf(1000,20),vyhat(1000),trudof,temp1,temp2
real vs(1000,20),veta(1000),vs0,vrss,vscale(1000)
real varco(20,20)
integer indco(20),numco
character*10 plotn
ans=0
10011 if(ans .eq. yes .or. ans .eq. no) goto 10012
write(6,10020)
10020 format( 48h this program computes standard deviation curves,/,
* 53h for the estimated functions, as well as a covariance,/,
* 49h matrix for the estimated constants in the model.,/,/, 55h
* the procedure takes n x (time for final gam iteration),/, 47h
* ( as you ll find out), so dont use carelessly.,/, 26h do you
*wish to continue ?)
read(5,10030)ans
if(ans .ne. yes)goto 10051
goto 10012
10051 continue
if(ans .ne. no)goto 10071
return
10071 continue
goto 10011
10012 continue
numco=1
j=1
goto 10083
10081 j=j+1
10083 if((j).gt.(p))goto 10082
if(spans(j) .ne. 1)goto 10101
numco=numco+1
indco(numco)=j
10101 continue
goto 10081
10082 continue
do 10111 j=1,numco
do 10121 jj=1,j
varco(j,jj)=0.0
10121 continue
10122 continue
10111 continue
10112 continue
parm1=.true.
parm2=.false.
its=bacnit
if(bacnit .le. 10)goto 10141
its=10
10141 continue
zero=0.0
trudof=0d0
vtotw=0.0
i=1
goto 10153
10151 i=i+1
10153 if((i).gt.(n))goto 10152
vyhat(i)=0d0
wz(i)=weight(w(i),muhat(i),ni(i))
vtotw=vtotw+wz(i)
j=1
goto 10163
10161 j=j+1
10163 if((j).gt.(p))goto 10162
vf(i,j)=0d0
goto 10161
10162 continue
goto 10151
10152 continue
ii=1
goto 10173
10171 ii=ii+1
10173 if((ii).gt.(n))goto 10172
if(wz(ii) .ne. 0.0)goto 10191
goto 10171
10191 continue
vs0=wz(ii)/vtotw
i=1
goto 10203
10201 i=i+1
10203 if((i).gt.(n))goto 10202
z(i)=0.0
veta(i)=vs0
j=1
goto 10213
10211 j=j+1
10213 if((j).gt.(p))goto 10212
vs(i,j)=0.0
goto 10211
10212 continue
goto 10201
10202 continue
z(ii)=1
call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,veta, vs,vs0,v
*rss,zero,nitb,parm2,parm1,its)
j=1
goto 10223
10221 j=j+1
10223 if((j).gt.(numco))goto 10222
if(j .ne. 1)goto 10241
sloup1=vs0
goto 10251
10241 continue
sloup1=slopes(indco(j))
10251 continue
10231 continue
jj=1
goto 10263
10261 jj=jj+1
10263 if((jj).gt.(j))goto 10262
if(jj .ne. 1)goto 10281
sloup2=vs0
goto 10291
10281 continue
sloup2=slopes(indco(jj))
10291 continue
10271 continue
varco(j,jj)=varco(j,jj)+sloup1*sloup2/wz(ii)
goto 10261
10262 continue
goto 10221
10222 continue
i=1
goto 10303
10301 i=i+1
10303 if((i).gt.(n))goto 10302
j=1
goto 10313
10311 j=j+1
10313 if((j).gt.(p))goto 10312
temp1=vs(i,j)
vf(i,j)=vf(i,j)+temp1*temp1/wz(ii)
goto 10311
10312 continue
temp2=veta(i)
vyhat(i)=vyhat(i)+temp2*temp2/wz(ii)
goto 10301
10302 continue
write(6,10320)ii,n,nitb
10320 format( 11h completed ,i3, 4h of ,i3, 5h with,i3, 11h it
*erations)
goto 10171
10172 continue
i=1
goto 10333
10331 i=i+1
10333 if((i).gt.(n))goto 10332
trudof=trudof+vyhat(i)*wz(i)
goto 10331
10332 continue
scale=1.0
write(6,10340)trudof
10340 format( 13h true dof is ,f10.2)
if(mod .ne. gauss)goto 10361
scale=sqrt(devian/(n-trudof))
10361 continue
j=1
goto 10373
10371 j=j+1
10373 if((j).gt.(numco))goto 10372
jj=1
goto 10383
10381 jj=jj+1
10383 if((jj).gt.(j))goto 10382
varco(j,jj)=varco(j,jj)*scale*scale
varco(jj,j)=varco(j,jj)
goto 10381
10382 continue
goto 10371
10372 continue
write(6,10390)
10390 format(/ 31h covariance matrix of constants,/)
write(6,10400)(varco(1,j),j=1,numco)
10400 format( 25h intercept s0 -----------,10(e9.2,1x))
jj=2
goto 10413
10411 jj=jj+1
10413 if((jj).gt.(numco))goto 10412
write(6,10420)(labels(k,index(indco(jj))),k=1,6),(varco(jj,j),j=1,
*numco)
10420 format (1x,6a4,10(e9.2,1x))
goto 10411
10412 continue
write(6,10430)
10430 format( 33h want to dump plots with +-2sd? )
read(5,10030)ans
if(ans .ne. no)goto 10451
return
10451 continue
10030 format(a1)
r2=(devnll-devian)*100.0/devnll
j=1
goto 10463
10461 j=j+1
10463 if((j).gt.(p))goto 10462
write(6,10470)(labels(k,index(j)),k=1,10)
10470 format ( 16h do you want ...,10a4)
read(5,10030)ans
if(ans .ne. no)goto 10491
goto 10461
10491 continue
call nplot(plotn)
totsm=0.0
do 10501 i=1,neffob(j)
totsm=totsm+s(tag(i,j),j)
10501 continue
10502 continue
totsm=totsm/neffob(j)
oldx=-9999
i=1
goto 10513
10511 i=i+1
10513 if((i).gt.(neffob(j)))goto 10512
ii=tag(i,j)
if(x(ii,j) .eq. oldx)goto 10531
sput=s(ii,j)-totsm
if(vf(ii,j) .le. 0.0)goto 10551
sputd=sqrt(vf(ii,j))
goto 10561
10551 continue
sputd=0.0
10561 continue
10541 continue
sputd=2*sputd*scale
sputl=sput-sputd
sputu=sput+sputd
write(2,*)x(ii,j),sput,sputl,sputu
oldx=x(ii,j)
10531 continue
goto 10511
10512 continue
write(6,10570)spans(j),dofs(j),trudof,devian,r2
goto 10461
10462 continue
10580 format( 7h data a,i1, 27h; input x yhat yhatl yhatu;)
10590 format( 34h title .c=blue estimated function; /, 49hproc gplo
*t; plot yhat*x yhatl*x yhatu*x/overlay;,/, 33hsymbol1 i=join
*c=green v=diamond;,/, 28hsymbol2 i=join c=red v=star;,/, 2
*8hsymbol3 i=join c=red v=star;)
10570 format( 5hspan=,f5.2, 5h dof=,f6.2, 19h true dof (overall
*),f6.2, 9hdeviance=, f8.2, 4h r2=,f8.2, 2h%;)
10600 format( 8hlabel x=,10a4, 13hyhat=s;cards;)
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine helpg
real a(20),paus,halt
data paus/ 4hpaus/, halt/ 4hstop/
data quit/ 1hq/
open(12,file=8hhelpunix)
rewind 12
10011 continue
read(12,10020)a
10020 format(20a4)
if(a(1) .ne. halt)goto 10041
return
goto10031
10041 if(a(1) .ne. paus)goto 10051
write(6,10060)
10060 format( 40h type q for quit, anything else for more)
read(5,10070)ans
10070 format(a1)
if(ans .ne. quit)goto 10091
return
10091 continue
goto 10101
10051 continue
write(6,10110)a
10101 continue
10031 continue
10110 format( 1h ,20a4)
goto 10011
10012 continue
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine info
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real match
integer indms,kms(1000)
common/mtchco/match,indms,kms
totdof=1.0
do 10011 j=1,p
totdof=totdof+dofs(j)
10011 continue
10012 continue
if(mod .ne. cox)goto 10031
totdof=totdof-1.0
10031 continue
errn=n-totdof
scale=1.0
if(mod .ne. gauss)goto 10051
scale=devian/errn
10051 continue
if(mod .ne. binom)goto 10071
write(6,10080)
10080 format(// 35h logistic additive regression model/)
write(7,10090)
10090 format(// 35h logistic additive regression model/)
goto10061
10071 if(mod .ne. gauss)goto 10101
write(6,10110)
10110 format(// 35h gaussian additive regression model/)
write(7,10120)
10120 format(// 35h gaussian additive regression model/)
goto10061
10101 if(mod .ne. match)goto 10131
write(6,10140)
10140 format(// 47h matched case-control additive regression model/)
write(7,10150)
10150 format(// 47h matched case-control additive regression model/)
goto10061
10131 if(mod .ne. cox)goto 10161
write(6,10170)
10170 format(// 50h additive cox proportional hazard regression model/
*)
write(7,10180)
10180 format(// 50h additive cox proportional hazard regression model/
*)
goto 10191
10161 continue
10191 continue
10061 continue
if(esttyp .ne. score)goto 10211
write(6,10220)
10220 format(// 28h estimation by local scoring/)
write(7,10230)
10230 format(// 28h estimation by local scoring/)
goto 10241
10211 continue
write(6,10250)
10250 format(// 31h estimation by local likelihood/)
write(7,10260)
10260 format(// 31h estimation by local likelihood/)
10241 continue
10201 continue
write(6,10270)(labels(k,indy),k=1,10)
10270 format (/ 23h response variable ....,10a4)
write(7,10280)(labels(k,indy),k=1,10)
10280 format (/ 23h response variable ....,10a4)
if(mod .ne. binom)goto 10301
if(indni .eq. 0)goto 10321
write(6,10330)(labels(k,indni),k=1,10)
10330 format ( 23h binomial denominator..,10a4)
write(7,10340)(labels(k,indni),k=1,10)
10340 format ( 23h binomial denominator..,10a4)
goto 10351
10321 continue
write(6,10360)
10360 format( 20h bernoulli response )
write(7,10370)
10370 format( 20h bernoulli response )
10351 continue
10311 continue
10301 continue
write(6,10380)devian,nitout,numsm
10380 format ( 12h deviance = ,f14.3, 14h # iterations=,i3, 20h #
*smooths/variable=,i4)
write(7,10390)devian,nitout,numsm
10390 format ( 12h deviance = ,f14.3, 14h # iterations=,i3, 20h #
*smooths/variable=,i4)
a=devian/n
write(6,10400)a
10400 format( 20h average deviance = , f12.3)
write(7,10410)a
10410 format( 20h average deviance = , f12.3)
write(6,10420)errn,scale
10420 format( 18h dof of deviance ,f6.2, 18h scale estimate ,f1
*0.3)
write(7,10430)errn,scale
10430 format( 18h dof of deviance ,f6.2, 18h scale estimate ,f1
*0.3)
r2=100*(devnll-devian)/devnll
write(6,10440)r2,devnll
10440 format ( 12h r square = ,f5.2, 24h% of a null deviance of ,f14
*.3)
write(7,10450)r2,devnll
10450 format ( 12h r square = ,f5.2, 24h% of a null deviance of ,f14
*.3)
write(6,10460)
10460 format(/, 65h variable span dof slope n
*ame ,/, 65h -------- ---- --- ------ ---
*------------------------)
write(7,10470)
10470 format(/, 65h variable span dof slope n
*ame ,/, 65h -------- ---- --- ------ ---
*------------------------)
if(mod .eq. cox)goto 10491
write(6,10500)s0
10500 format ( 23h 0 -- 1 ,f10.4, 26h s0---the inte
*rcept term)
write(7,10510)s0
10510 format ( 23h 0 -- 1 ,f10.4, 26h s0---the inte
*rcept term)
10491 continue
j=1
goto 10523
10521 j=j+1
10523 if((j).gt.(p))goto 10522
if(spans(j) .lt. 1.0)goto 10541
write(6,10550)index(j),spans(j), dofs(j),slopes(j),modlab(j),(labe
*ls(k,index(j)),k=1,10)
10550 format ( 1h ,i3,7x,f6.2,2x, f5.3, f9.4,3x,11a4)
write(7,10560)index(j),spans(j), dofs(j),slopes(j),modlab(j),(labe
*ls(k,index(j)),k=1,10)
10560 format ( 1h ,i3,7x,f6.2,2x, f5.3, f9.4,3x,11a4)
goto 10571
10541 continue
write(6,10580)index(j),spans(j),dofs(j),modlab(j),(labels(k,index(
*j)),k=1,10)
10580 format ( 1h ,i3,7x,f6.2,2x, f6.3, 8h smooth,3x,11a4)
write(7,10590)index(j),spans(j),dofs(j),modlab(j),(labels(k,index(
*j)),k=1,10)
10590 format ( 1h ,i3,7x,f6.2,2x, f6.3, 8h smooth,3x,11a4)
10571 continue
10531 continue
goto 10521
10522 continue
write(6,10600)totdof
10600 format ( 29h ----- ,/, 19h
* ,f6.3)
write(7,10610)totdof
10610 format ( 29h ----- ,/, 19h
* ,f6.3)
if(.not.(vlink))goto 10631
slinc=slinc-devian
write(6,10640)sldof,slinc
10640 format( 21h variable link (dof= ,f5.2, 8h) caused, / 21h d
*eviance to drop by ,f10.3)
write(7,10650)sldof,slinc
10650 format( 21h variable link (dof= ,f5.2, 8h) caused, / 21h d
*eviance to drop by ,f10.3)
10631 continue
return
end
real function devu(n,fits,y,w,ni)
real fits(1),y(1),w(1),ni(1)
rss=0.0
do 10661 i=1,n
rat=y(i)/abs(fits(i))
if(rat .gt. 0.0)goto 10681
write(6,10690)
10690 format( 6hshriek)
stop
10681 continue
ratl=alog(rat)
rss=rss+w(i)*(rat-ratl -1)*2.0
10661 continue
10662 continue
devu=rss
return
end
real function linku(n,etahat,muhat)
linku=etahat
return
end
real function linvu(muhat)
real muhat
linvu=muhat
return
end
real function dirivu(muhat)
real muhat
dirivu=1.0
return
end
real function weigtu(w,muhat,ni)
real w,muhat,ni
t=muhat/700.00
t=t*t
if(t .gt. 0.0001)goto 10711
t=0.0001
10711 continue
weigtu=1.0/t
return
end
real function fpinv(arg)
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
if(.not.(vflag))goto 10731
fpinv=1.0/yinter(argf,sldir,n,arg)
goto 10741
10731 continue
fpinv=1.0
10741 continue
10721 continue
return
end
real function flink(arg)
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
if(.not.(vflag))goto 10761
flink=yinter(argf,slink,n,arg)
goto 10771
10761 continue
flink=arg
10771 continue
10751 continue
return
end
subroutine montne(x,n)
real x(n)
integer bb,eb,br,er,bl,el
bb=0
eb=bb
10780 continue
10781 if(eb.ge.n) goto 10782
bb=eb+1
eb=bb
10791 if(eb .ge. n) goto 10792
if(x(bb) .ne. x(eb+1))goto 10811
eb=eb+1
goto 10821
10811 continue
goto 10792
10821 continue
10801 continue
goto 10791
10792 continue
10831 continue
if(eb .ge. n)goto 10851
if(x(eb) .le. x(eb+1))goto 10871
br=eb+1
er=br
10881 if(er .ge. n) goto 10882
if(x(er+1) .ne. x(br))goto 10901
er=er+1
goto 10911
10901 continue
goto 10882
10911 continue
10891 continue
goto 10881
10882 continue
pmn=(x(bb)*(eb-bb+1)+x(br)*(er-br+1))/(er-bb+1)
eb=er
do 10921 i=bb,eb
x(i)=pmn
10921 continue
10922 continue
10871 continue
10851 continue
if(bb.le.1) goto 10781
if(x(bb-1).le.x(bb)) goto 10781
bl=bb-1
el=bl
10931 if(bl .le. 1) goto 10932
if(x(bl-1) .ne. x(el))goto 10951
bl=bl-1
goto 10961
10951 continue
goto 10932
10961 continue
10941 continue
goto 10931
10932 continue
pmn=(x(bb)*(eb-bb+1)+x(bl)*(el-bl+1))/(eb-bl+1)
bb=bl
do 10971 i=bb,eb
x(i)=pmn
10971 continue
10972 continue
goto 10831
10832 continue
goto 10781
10782 continue
return
end
subroutine der (n,x,s,w,fdel,d,sc)
real x(n),s(n),w(n),d(n),sc(n,3)
integer bl,el,bc,ec,br,er
if(x(n) .gt. x(1))goto 10991
do 11001 j=1,n
d(j)=0.0
11001 continue
11002 continue
return
10991 continue
i=n/4
j=3*i
scale=x(j)-x(i)
11011 if(scale.gt.0.0) goto 11012
if(j.lt.n) j=j+1
if(i.gt.1) i=i-1
scale=x(j)-x(i)
goto 11011
11012 continue
del=fdel*scale*2.0
do 11021 j=1,n
sc(j,1)=x(j)
sc(j,2)=s(j)
sc(j,3)=w(j)
11021 continue
11022 continue
call pool (n,sc,sc(1,2),sc(1,3),del)
bc=0
br=bc
er=br
11031 continue
br=er+1
er=br
11041 if(er .ge. n) goto 11042
if(sc(br,1) .ne. sc(er+1,1))goto 11061
er=er+1
goto 11071
11061 continue
goto 11042
11071 continue
11051 continue
goto 11041
11042 continue
if(br .ne. 1)goto 11091
bl=br
el=er
goto 11031
11091 continue
if(bc .ne. 0)goto 11111
bc=br
ec=er
do 11121 j=bl,el
d(j)=(sc(bc,2)-sc(bl,2))/(sc(bc,1)-sc(bl,1))
11121 continue
11122 continue
goto 11031
11111 continue
do 11131 j=bc,ec
d(j)=(sc(br,2)-sc(bl,2))/(sc(br,1)-sc(bl,1))
11131 continue
11132 continue
if(er .ne. n)goto 11151
do 11161 j=br,er
d(j)=(sc(br,2)-sc(bc,2))/(sc(br,1)-sc(bc,1))
11161 continue
11162 continue
goto 11032
11151 continue
bl=bc
el=ec
bc=br
ec=er
goto 11031
11032 continue
return
end
subroutine pool (n,x,y,w,del)
real x(n),y(n),w(n)
integer bb,eb,br,er,bl,el
bb=0
eb=bb
11171 if(eb.ge.n) goto 11172
bb=eb+1
eb=bb
11181 if(eb .ge. n) goto 11182
if(x(bb) .ne. x(eb+1))goto 11201
eb=eb+1
goto 11211
11201 continue
goto 11182
11211 continue
11191 continue
goto 11181
11182 continue
if(eb .ge. n)goto 11231
if(x(eb+1)-x(eb) .ge. del)goto 11251
br=eb+1
er=br
11261 if(er .ge. n) goto 11262
if(x(er+1) .ne. x(br))goto 11281
er=er+1
goto 11291
11281 continue
goto 11262
11291 continue
11271 continue
goto 11261
11262 continue
if(x(er+1)-x(er).lt.x(eb+1)-x(eb))goto 11171
eb=er
pw=w(bb)+w(eb)
px=(x(bb)*w(bb)+x(eb)*w(eb))/pw
py=(y(bb)*w(bb)+y(eb)*w(eb))/pw
do 11301 i=bb,eb
x(i)=px
y(i)=py
w(i)=pw
11301 continue
11302 continue
11251 continue
11231 continue
11311 continue
if(bb.le.1)goto 11312
if(x(bb)-x(bb-1).ge.del)goto 11312
bl=bb-1
el=bl
11321 if(bl .le. 1) goto 11322
if(x(bl-1) .ne. x(el))goto 11341
bl=bl-1
goto 11351
11341 continue
goto 11322
11351 continue
11331 continue
goto 11321
11322 continue
bb=bl
pw=w(bb)+w(eb)
px=(x(bb)*w(bb)+x(eb)*w(eb))/pw
py=(y(bb)*w(bb)+y(eb)*w(eb))/pw
do 11361 i=bb,eb
x(i)=px
y(i)=py
w(i)=pw
11361 continue
11362 continue
goto 11311
11312 continue
goto 11171
11172 continue
return
end
real function yinter(z,fz,n,x)
real z(n),fz(n),x
integer n
if(x .gt. z(1))goto 11381
ii=1
jj=2
11391 if(fz(jj) .ne. fz(ii)) goto 11392
jj=jj+1
goto 11391
11392 continue
goto11371
11381 if(x .lt. z(n))goto 11401
jj=n
ii=n-1
11411 if(fz(ii) .ne. fz(jj)) goto 11412
ii=ii-1
goto 11411
11412 continue
goto 11421
11401 continue
ii=1
11431 if((z(ii).le.x).and.(x.le.z(ii+1))) goto 11432
ii=ii+1
goto 11431
11432 continue
jj=ii+1
11421 continue
11371 continue
if(z(ii) .ne. z(jj))goto 11451
yinter=(fz(ii)+fz(jj))/2
goto 11461
11451 continue
slope=(fz(jj)-fz(ii))/(z(jj)-z(ii))
yinter=fz(ii)+slope*(x-z(ii))
11461 continue
11441 continue
return
end
subroutine mysort(input,tag,from,till)
real input(1)
integer tag(1),from,till
real scrat(1000)
i=from
goto 11473
11471 i=i+1
11473 if((i).gt.(till))goto 11472
tag(i)=i
scrat(i)=input(i)
goto 11471
11472 continue
call sort(scrat,tag,from,till)
return
end
subroutine psort(input,tag,from,till)
real input(1)
integer tag(1),from,till
real scrat(1000)
i=from
goto 11483
11481 i=i+1
11483 if((i).gt.(till))goto 11482
scrat(i)=input(i)
goto 11481
11482 continue
call sort(scrat,tag,from,till)
return
end
subroutine sort (v,a,ii,jj)
c
c puts into a the permutation vector which sorts v into
c increasing order. only elements from ii to jj are considered.
c arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements
c v is returned sorted
c this is a modification of cacm algorithm #347 by r. c. singleton,
c which is a modified hoare quicksort.
c
dimension a(jj),v(1),iu(20),il(20)
integer t,tt
integer a
real v
m=1
i=ii
j=jj
10 if (i.ge.j) go to 80
20 k=i
ij=(j+i)/2
t=a(ij)
vt=v(ij)
if (v(i).le.vt) go to 30
a(ij)=a(i)
a(i)=t
t=a(ij)
v(ij)=v(i)
v(i)=vt
vt=v(ij)
30 l=j
if (v(j).ge.vt) go to 50
a(ij)=a(j)
a(j)=t
t=a(ij)
v(ij)=v(j)
v(j)=vt
vt=v(ij)
if (v(i).le.vt) go to 50
a(ij)=a(i)
a(i)=t
t=a(ij)
v(ij)=v(i)
v(i)=vt
vt=v(ij)
go to 50
40 a(l)=a(k)
a(k)=tt
v(l)=v(k)
v(k)=vtt
50 l=l-1
if (v(l).gt.vt) go to 50
tt=a(l)
vtt=v(l)
60 k=k+1
if (v(k).lt.vt) go to 60
if (k.le.l) go to 40
if (l-i.le.j-k) go to 70
il(m)=i
iu(m)=l
i=k
m=m+1
go to 90
70 il(m)=k
iu(m)=j
j=l
m=m+1
go to 90
80 m=m-1
if (m.eq.0) return
i=il(m)
j=iu(m)
90 if (j-i.gt.10) go to 20
if (i.eq.ii) go to 10
i=i-1
100 i=i+1
if (i.eq.j) go to 80
t=a(i+1)
vt=v(i+1)
if (v(i).le.vt) go to 100
k=i
110 a(k+1)=a(k)
v(k+1)=v(k)
k=k-1
if (vt.lt.v(k)) go to 110
a(k+1)=t
v(k+1)=vt
go to 100
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine init(linv)
external linv
real linv
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
real match
integer indms,kms(1000)
common/mtchco/match,indms,kms
theta=1.0
alpha=1.0
maxit=20
if(mod .ne. gauss .or. .not.(.not.vlink))goto 10021
maxit=1
10021 continue
s0=0
sn=0
i=1
goto 10033
10031 i=i+1
10033 if((i).gt.(n))goto 10032
s0=s0+w(i)*y(i)
sn=sn+w(i)
goto 10031
10032 continue
p0=s0/sn
i=1
goto 10043
10041 i=i+1
10043 if((i).gt.(n))goto 10042
muhat(i)=p0
etahat(i)=linv(muhat(i))
do 10051 j=1,p
s(i,j)=0.0
10051 continue
10052 continue
goto 10041
10042 continue
s0=etahat(1)
if(.not.(vlink))goto 10071
do 10081 i=1,n
sldir(i)=1.0
10081 continue
10082 continue
vflag=.false.
10071 continue
return
end
real function devg(n,fits,y,w,ni)
real fits(1),y(1),w(1),ni(1)
rss=0.0
do 10091 i=1,n
rss=rss+w(i)*(y(i)-fits(i))*(y(i)-fits(i))
10091 continue
10092 continue
devg=rss
return
end
function devb(n,fits,y,w,ni)
real fits(n),y(n),ni(n),w(1)
dev=0.0
ent1=0
ent2=0
do 10101 i=1,n
pr=fits(i)
if(pr .ge. 0.0001)goto 10121
pr=0.0001
10121 continue
if(pr .le. 0.9999)goto 10141
pr=0.9999
10141 continue
if(((1.0-y(i))*y(i)) .gt. 0.0)goto 10161
entrop=0.0
goto 10171
10161 continue
entrop=2.0*w(i)*ni(i)*(y(i)*alog(y(i))+ (1.0-y(i))*alog(1.0-y(i))
*)
10171 continue
10151 continue
entadd=2.*w(i)*ni(i)*(y(i)*alog(pr)+(1.0-y(i))*alog(1.0-pr))
dev=dev+entrop-entadd
ent1=ent1+entrop
ent2=ent2+entadd
10101 continue
10102 continue
devb=dev
return
end
function devmc(n,fits,y,w,ni)
real fits(n),y(n),ni(n),w(1)
dev=0.0
do 10181 i=1,n
pr=fits(i)
if(pr .ge. 0.0001)goto 10201
pr=0.0001
10201 continue
if(pr .le. 0.9999)goto 10221
pr=0.9999
10221 continue
entrop=2.*w(i)*ni(i)*y(i)*alog(pr)
dev=dev-entrop
10181 continue
10182 continue
devmc=dev
return
end
subroutine linkg(n,etahat,muhat)
real muhat(1),etahat(1)
do 10231 i=1,n
muhat(i)=etahat(i)
10231 continue
10232 continue
return
end
real function linvg(muhat)
real muhat
linvg=muhat
return
end
subroutine linkmc(n,etahat,muhat)
real muhat(1),etahat(1),work(1000)
real match
integer indms,kms(1000)
common/mtchco/match,indms,kms
do 10241 i=1,n
muhat(i)=exp(etahat(i))
10241 continue
10242 continue
do 10251 i=1,n
work(i)=0.0
10251 continue
10252 continue
do 10261 i=1,n
work(kms(i))=work(kms(i))+muhat(i)
10261 continue
10262 continue
do 10271 i=1,n
muhat(i)=muhat(i)/work(kms(i))
10271 continue
10272 continue
return
end
real function linvmc(muhat)
real muhat
linvmc=alog(muhat)
return
end
subroutine linkb(n,etahat,muhat)
real muhat(1),etahat(1)
do 10281 i=1,n
pr=exp(etahat(i))
muhat(i)=pr/(1+pr)
10281 continue
10282 continue
return
end
real function linvb(muhat)
real muhat
real logit
logit=muhat
d=1.0-logit
if(d .ge. 0.0001)goto 10301
d=0.0001
10301 continue
logit=logit/d
if(logit .ge. 0.0001)goto 10321
linvb=alog(0.0001)
goto10311
10321 if(logit .le. 0 9999.0)goto 10331
linvb=alog(9999.0)
goto 10341
10331 continue
linvb=alog(logit)
10341 continue
10311 continue
return
end
real function dirivg(muhat)
real muhat
dirivg=1.0
return
end
real function dirivb(muhat)
real muhat
pr=muhat
pr=pr*(1.0-pr)
if(pr .gt. 0.01)goto 10361
pr=0.01
10361 continue
dirivb=1.0/pr
return
end
real function weigtg(w,muhat,ni)
real muhat,ni
weigtg=w
return
end
real function weigtb(w,muhat,ni)
real muhat,ni
pr=dirivb(muhat)
weigtb=w*ni/pr
return
end
subroutine title
write(6,10370)
10370 format (// 64h ******* *** ****** *
* ** ,/, 64h ********* ***** ******
* ** ** ,/, 64h ** ** ** **
* ** *** *** ,/, 64h ** **
*** ** **** **** ,/, 64h **
* ** ** ** ** *** ** ,/, 64h **
***** ********* ** ** * ** ,/, 64h
* ** **** ********* ** ** ** ,/,
* 64h ** ** ** ** ** ** **
* ,/, 64h ********* ** ** ****** ** *
** ,/, 64h ******* ** ** ******
*** ** ,/, 62h
* ,/, 51h generalized additive
* interactive models,//, 53h coded by trevor hastie
* ,/, 61h and robert tibshirani
* (SLAC aug 25 1984),/, 61h last update
* (AT&T labs apr 8 1987),/, //)
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine initcx
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
do 10011 i=1,n
etahat(i)=0
muhat(i)=0
10011 continue
10012 continue
i=1
10021 if(icensq(i) .eq. 1) goto 10022
i=i+1
goto 10021
10022 continue
tq(1)=y(i)
iriskq(1)=i
jj=1
i=i+1
10031 if(i .gt. nq) goto 10032
if(y(i) .eq. y(i-1) .or. icensq(i) .ne. 1)goto 10051
jj=jj+1
tq(jj)=y(i)
iriskq(jj)=i
10051 continue
i=i+1
goto 10031
10032 continue
kq=jj
call calcdd
return
end
subroutine calcdd
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
do 10061 i=1,kq
ddq(i)=0
ii=iriskq(i)
if(i .ge. kq)goto 10081
ij=i+1
ik=iriskq(ij)-1
goto 10091
10081 continue
ik=nq
10091 continue
10071 continue
do 10101 j=ii,ik
ddq(i)=ddq(i)+icensq(j)
10101 continue
10102 continue
10061 continue
10062 continue
return
end
real function devc(n,fits,y,w,ni)
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real fits(1),y(1),w(1),ni(1)
call calcsf(fits)
call calcll(fits,value)
devc=value
return
end
subroutine calcll(fits,value)
real fits(1)
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
j=1
sum2=0
10111 if(icensq(j) .ne. 0) goto 10112
j=j+1
goto 10111
10112 continue
do 10121 i=j,nq
temp=fits(i)
termq(i)=exp(temp)
sum2=sum2+termq(i)
10121 continue
10122 continue
sum1=0
do 10131 i=1,kq
if(i .eq. 1)goto 10151
is=i-1
it=iriskq(is)
iu=iriskq(i)-1
do 10161 j=it,iu
sum2=sum2-termq(j)
10161 continue
10162 continue
10151 continue
sum1=sum1+sfitsq(i)-ddq(i)*log(sum2)
10131 continue
10132 continue
value=-2*sum1
return
end
subroutine calcsf(fits)
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real fits(1)
do 10171 i=1,kq
sfitsq(i)=0
ii=iriskq(i)
if(i .ge. kq)goto 10191
ij=i+1
ik=iriskq(ij)-1
goto 10201
10191 continue
ik=nq
10201 continue
10181 continue
do 10211 j=ii,ik
sfitsq(i)=sfitsq(i)+icensq(j)*fits(j)
10211 continue
10212 continue
10171 continue
10172 continue
return
end
subroutine calcz
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real z(1000),wz(1000)
common /wvc/z,wz
real argf(1000),slink(1000),sldir(1000),sleta(1000),sldof,slinc
integer sltag(1000)
logical vlink,vflag
common /fff/argf,slink,sldir,sltag,vlink,vflag,sleta,sldof,slinc
do 10221 i=1,nq
irindq(i)=0
10221 continue
10222 continue
do 10231 i=1,kq
irindq(iriskq(i))=i
10231 continue
10232 continue
aq(1)=0
ii=iriskq(1)
do 10241 i=ii,nq
aq(1)=aq(1)+exp(etahat(i))
10241 continue
10242 continue
do 10251 i=2,kq
ii=i-1
aq(i)=aq(ii)
it=i-1
iv=iriskq(it)
is=iriskq(i)-1
do 10261 j=iv,is
aq(i)=aq(i)-exp(etahat(j))
10261 continue
10262 continue
10251 continue
10252 continue
bq(nq)=0
bbq(nq)=0
do 10271 i=1,kq
bq(nq)=bq(nq)+ddq(i)*(1.0/aq(i))
bbq(nq)=bbq(nq)+ddq(i)*(1.0/aq(i))**2
10271 continue
10272 continue
ii=nq
do 10281 i=2,nq
ii=ii-1
ij=ii+1
bq(ii)=bq(ij)
bbq(ii)=bbq(ij)
if(irindq(ij) .eq. 0)goto 10301
bq(ii)=bq(ii)-ddq(irindq(ij))*(1.0/aq(irindq(ij)))
bbq(ii)=bbq(ii)-ddq(irindq(ij))*(1.0/aq(irindq(ij)))**2
10301 continue
10281 continue
10282 continue
do 10311 i=1,nq
r1q(i)=icensq(i)-exp(etahat(i))*bq(i)
r2q(i)=-exp(etahat(i))*bq(i)+exp(2*etahat(i))*bbq(i)
10311 continue
10312 continue
do 10321 i=1,nq
z(i)=etahat(i)-r1q(i)/r2q(i)
10321 continue
10322 continue
do 10331 i=1,nq
wz(i)=-1.0*r2q(i)
10331 continue
10332 continue
return
end
subroutine stpar
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real termq(1000),r1q(1000),r2q(1000),tq(1000),sfitsq(1000),aq(1000
*),bq(1000), bbq(1000),etaq(1000),sq(1000,20),wq(1000),ddq(1000)
integer icensq(1000),irindq(1000),iriskq(1000),nq,kq
common /cox/termq,r1q,r2q,tq,sfitsq,aq,bq,bbq,etaq,sq,s0q,wq,ddq,
*icensq,irindq,iriskq,nq,kq
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
s0=0
do 10341 j=1,p
slopes(j)=0
10341 continue
10342 continue
do 10351 i=1,n
do 10361 j=1,p
s(i,j)=slopes(j)*x(i,j)
10361 continue
10362 continue
10351 continue
10352 continue
do 10371 j=1,p
temp=0
do 10381 i=1,n
temp=temp+s(i,j)/n
10381 continue
10382 continue
do 10391 i=1,n
s(i,j)=s(i,j)-temp
10391 continue
10392 continue
10371 continue
10372 continue
do 10401 i=1,n
muhat(i)=0
do 10411 j=1,p
muhat(i)=muhat(i)+s(i,j)
10411 continue
10412 continue
10401 continue
10402 continue
do 10421 i=1,n
etahat(i)=muhat(i)
10421 continue
10422 continue
do 10431 i=1,n
etaq(i)=0
10431 continue
10432 continue
s0q=0
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine multam
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real z(1000),wz(1000)
common /wvc/z,wz
real wm(15,15),wm1(1000)
real scrat1(1000),scrat2(1000)
itss=30
call initm
devnew=devm(0)
devnll=devnew
devold=10*devnew+10
nit=0
numsm=0
10011 if(abs((devold-devnew)/devold) .le. thresh .or. nit .ge. maxit) go
*to 10012
nit=nit+1
call getz
i=1
goto 10023
10021 i=i+1
10023 if((i).gt.(n))goto 10022
z(i)=z(i)+etahat(i)
goto 10021
10022 continue
call bakfit(n,p,z,x,wz,tag,spans,dofs,slopes,cross,etahat, s,s0,r
*ss,thbak,nitb,verbos,quiet,itss)
do 10031 i=1,n
etahat(i)=etahat(i)-s0
10031 continue
10032 continue
j=1
goto 10043
10041 j=j+1
10043 if((j).gt.(mcat-1))goto 10042
thetak(j)=thetak(j)+s0
goto 10041
10042 continue
s0=0
numsm=numsm+nitb
call thetaf
devold=devnew
devnew=devm(1)
if(.not.(.not.quiet))goto 10061
write(6,10070)nit,devnew
10070 format(/10x, 22houter loop (gam) nit= , i3, 5h dev=,f16.5,//
*)
10061 continue
goto 10011
10012 continue
devian=devnew
nitout=nit
return
end
subroutine initm
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real tots(15)
real linvb
maxit=20
do 10081 k=1,mcat
tots(k)=0
10081 continue
10082 continue
totn=0.0
i=1
goto 10093
10091 i=i+1
10093 if((i).gt.(n))goto 10092
k=1
goto 10103
10101 k=k+1
10103 if((k).gt.(mcat))goto 10102
tots(k)=tots(k)+w(i)*ni(i)*ym(i,k)
goto 10101
10102 continue
totn=totn+w(i)*ni(i)
goto 10091
10092 continue
tots(1)=tots(1)/totn
do 10111 k=2,mcat
tots(k)=tots(k)/totn+tots(k-1)
10111 continue
10112 continue
if(abs(tots(mcat)-1.0) .le. 00001)goto 10131
write(6,10140)
10140 format( 18h oops, check initm)
10131 continue
k=1
goto 10153
10151 k=k+1
10153 if((k).gt.(mcat-1))goto 10152
thetak(k)=linvb(tots(k))
goto 10151
10152 continue
i=1
goto 10163
10161 i=i+1
10163 if((i).gt.(n))goto 10162
etahat(i)=0.0
do 10171 j=1,p
s(i,j)=0.0
10171 continue
10172 continue
goto 10161
10162 continue
return
end
real function devm(inp)
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real pm(15),pim(15)
real pr
dev=0.0
ent1=0
ent2=0
k=2
goto 10183
10181 k=k+1
10183 if((k).gt.(mcat-1))goto 10182
if(thetak(k) .ge. thetak(k-1))goto 10201
thetak(k)=thetak(k-1)*(1.00001)
write(6,10210)
10210 format( 27h oops, thetas not ascending)
10201 continue
goto 10181
10182 continue
i=1
goto 10223
10221 i=i+1
10223 if((i).gt.(n))goto 10222
call linkm(mcat,thetak,etahat(i),pim,pm)
k=1
goto 10233
10231 k=k+1
10233 if((k).gt.(mcat))goto 10232
pr=pm(k)
yk=ym(i,k)
if(pr .ge. 0.0001)goto 10251
pr=0.0001
10251 continue
if(yk .gt. 0.000001)goto 10271
entrop=0.0
goto 10281
10271 continue
entrop=2.0*w(i)*ni(i)*yk*alog(yk)
10281 continue
10261 continue
entadd=2.*w(i)*ni(i)*yk*alog(pr)
dev=dev+entrop-entadd
goto 10231
10232 continue
goto 10221
10222 continue
devm=dev
return
end
subroutine linkm(mcat,thetak,etahat,pim,pm)
integer mcat
real thetak(mcat),etahat,pim(mcat),pm(mcat)
eta=thetak(1)+etahat
call linkb(1,eta,pm(1))
pim(1)=pm(1)
k=2
goto 10293
10291 k=k+1
10293 if((k).gt.(mcat-1))goto 10292
eta=thetak(k)+etahat
call linkb(1,eta,pim(k))
pm(k)=pim(k)-pim(k-1)
goto 10291
10292 continue
pm(mcat)=1.0-pim(mcat-1)
return
end
subroutine getz
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real z(1000),wz(1000)
common /wvc/z,wz
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real wm1(15,15),ystar(15)
real wm2(15),wystar(15)
real pim(15),pm(15)
do 10301 k=1,mcat
wystar(k)=0.0
wm2(k)=0.0
10301 continue
10302 continue
i=1
goto 10313
10311 i=i+1
10313 if((i).gt.(n))goto 10312
call linkm(mcat,thetak,etahat(i),pim,pm)
call getw(mcat,ni(i),wm1,pm,pim)
ystar(1)=ym(i,1)-pm(1)
k=2
goto 10323
10321 k=k+1
10323 if((k).gt.(mcat-1))goto 10322
ystar(k)=ystar(k-1)+ym(i,k)-pm(k)
goto 10321
10322 continue
k=1
goto 10333
10331 k=k+1
10333 if((k).gt.(mcat-1))goto 10332
c=pim(k)*(1.0-pim(k))
ystar(k)=ystar(k)/c
goto 10331
10332 continue
k=1
goto 10343
10341 k=k+1
10343 if((k).gt.(mcat-1))goto 10342
wm2(k)=0.0
kk=1
goto 10353
10351 kk=kk+1
10353 if((kk).gt.(mcat-1))goto 10352
wm2(k)=wm2(k)+wm1(k,kk)
goto 10351
10352 continue
goto 10341
10342 continue
z(i)=0.0
wz(i)=0.0
k=1
goto 10363
10361 k=k+1
10363 if((k).gt.(mcat-1))goto 10362
z(i)=z(i)+ystar(k)*wm2(k)
wz(i)=wz(i)+wm2(k)
goto 10361
10362 continue
z(i)=z(i)/wz(i)
goto 10311
10312 continue
return
end
subroutine thetaf
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real wm1(15,15),ystar(15)
real wm2(15,15),wystar(15)
real pim(15),pm(15)
do 10371 k=1,mcat
wystar(k)=0.0
do 10381 kk=1,mcat
wm2(k,kk)=0.0
10381 continue
10382 continue
10371 continue
10372 continue
i=1
goto 10393
10391 i=i+1
10393 if((i).gt.(n))goto 10392
call linkm(mcat,thetak,etahat(i),pim,pm)
call getw(mcat,ni(i),wm1,pm,pim)
ystar(1)=ym(i,1)-pm(1)
k=2
goto 10403
10401 k=k+1
10403 if((k).gt.(mcat-1))goto 10402
ystar(k)=ystar(k-1)+ym(i,k)-pm(k)
goto 10401
10402 continue
k=1
goto 10413
10411 k=k+1
10413 if((k).gt.(mcat-1))goto 10412
c=pim(k)*(1.0-pim(k))
ystar(k)=ystar(k)/c + thetak(k)
goto 10411
10412 continue
k=1
goto 10423
10421 k=k+1
10423 if((k).gt.(mcat-1))goto 10422
kk=1
goto 10433
10431 kk=kk+1
10433 if((kk).gt.(mcat-1))goto 10432
wystar(k)=wystar(k)+wm1(k,kk)*ystar(kk)
wm2(k,kk)=wm2(k,kk)+wm1(k,kk)
goto 10431
10432 continue
goto 10421
10422 continue
goto 10391
10392 continue
call tholve(mcat,wm2,wystar,thetak)
return
end
subroutine tholve(mcat,wm,wystar,thetak)
integer mcat
real wm(15,mcat),wystar(mcat),thetak(mcat)
k=1
goto 10443
10441 k=k+1
10443 if((k).gt.(mcat-1))goto 10442
wm(mcat,k)=wystar(k)
wm(k,mcat)=wystar(k)
goto 10441
10442 continue
wm(mcat,mcat)=1.0
k=1
goto 10453
10451 k=k+1
10453 if((k).gt.(mcat-1))goto 10452
call sweep(15,mcat,k,wm)
goto 10451
10452 continue
k=1
goto 10463
10461 k=k+1
10463 if((k).gt.(mcat-1))goto 10462
thetak(k)=-wm(mcat,k)
goto 10461
10462 continue
return
end
subroutine sweep(maxp,p,k,a)
integer p
real a(maxp,p)
d=a(k,k)
do 10471 j=1,p
a(k,j)=a(k,j)/d
10471 continue
10472 continue
i=1
goto 10483
10481 i=i+1
10483 if((i).gt.(p))goto 10482
if(i .ne. k)goto 10501
goto 10481
10501 continue
b=a(i,k)
j=1
goto 10513
10511 j=j+1
10513 if((j).gt.(p))goto 10512
a(i,j)=a(i,j)-b*a(k,j)
goto 10511
10512 continue
a(i,k)=-b/d
goto 10481
10482 continue
a(k,k)=1/d
return
end
subroutine minfo
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real muhat(1000),etahat(1000)
real s(1000,20),thresh,slopes(20),dofs(20),devian,thbak
logical verbos,quiet
real s0,devnll
integer maxit,nitout,numsm
common /fitss/muhat,etahat,s,thresh,slopes,dofs,devian,thbak, s0,d
*evnll, maxit,verbos,quiet,nitout,numsm
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
totdof=mcat-1.0
do 10521 j=1,p
totdof=totdof+dofs(j)
10521 continue
10522 continue
totn=(mcat-1)*n
errn=totn-totdof
scale=1.0
write(6,10530)
10530 format(// 35h logistic additive regression model/, 42h propor
*tional odds model (mc cullagh 1980)/)
write(6,10540)
10540 format(// 28h estimation by local scoring/)
write(6,10550)devian,nitout,numsm
10550 format ( 12h deviance = ,f14.3, 14h # iterations=,i3, 20h #
*smooths/variable=,i4)
a=devian/totn
write(6,10560)a
10560 format( 20h average deviance = , f12.3)
write(6,10570)errn,scale
10570 format( 18h dof of deviance ,f8.2, 18h scale estimate ,f1
*0.3)
r2=100*(devnll-devian)/devnll
write(6,10580)r2,devnll
10580 format ( 12h r square = ,f5.2, 24h% of a null deviance of ,f14
*.3)
write(6,10590)
10590 format(/, 65h variable span dof slope n
*ame ,/, 65h -------- ---- --- ------ ---
*------------------------)
write(6,10600)thetak(1),(labels(k,indrm(1)),k=1,10)
10600 format( 26h theta(1) 1 ,f9.4,4x,10a4)
if(group .ne. yes)goto 10621
k=2
goto 10633
10631 k=k+1
10633 if((k).gt.(mcat-1))goto 10632
write(6,10640)k,thetak(k),(labels(kk,indrm(k)),kk=1,10)
10640 format( 7h theta(,i2, 17h) 1 ,f9.4,4x,10a4)
goto 10631
10632 continue
goto 10651
10621 continue
k=2
goto 10663
10661 k=k+1
10663 if((k).gt.(mcat-1))goto 10662
write(6,10670)k,thetak(k)
10670 format( 7h theta(,i2, 17h) 1 ,f9.4)
goto 10661
10662 continue
10651 continue
10611 continue
j=1
goto 10683
10681 j=j+1
10683 if((j).gt.(p))goto 10682
if(spans(j) .lt. 1.0)goto 10701
write(6,10710)index(j),spans(j), dofs(j),slopes(j),modlab(j),(labe
*ls(k,index(j)),k=1,10)
10710 format ( 1h ,i3,7x,f6.2,2x, f5.3,1x, f9.4,5x,11a4)
goto 10721
10701 continue
write(6,10730)index(j),spans(j),dofs(j),modlab(j),(labels(k,index(
*j)),k=1,10)
10730 format ( 1h ,i3,7x,f6.2,2x, f6.3, 8h smooth,5x,11a4)
10721 continue
10691 continue
goto 10681
10682 continue
write(6,10740)totdof
10740 format ( 29h ----- ,/, 19h
* ,f6.3)
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine multin
real dumlab(30),modlab(20)
integer indat(20),indmod(20)
common/extra/dumlab,modlab,indat,indmod
real mult,group,ym(1000,15)
real thetak(15)
integer mcat,indrm(15)
common/mltcmn/mult,group,thetak,ym,mcat,indrm
real datam(1000,20),labels(10,20)
integer maxp,n
common /datac/datam,labels,maxp,n
real x(1000,20),xt(1000,20),y(1000),ni(1000),w(1000),spans(20)
real answer,yes,no,binom,gauss,cox,mod,user,esttyp,local,score
integer tag(1000,20),tagy(1000),p,index(20),indy,indni,indw,llinc
logical cross(20)
common/modell/x,xt,y,ni,w,spans, answer,yes,no,binom,gauss,cox,mod
*,user,esttyp, tag,tagy,p,index,indy,indni,indw,cross,llinc,local,s
*core
j=1
goto 10013
10011 j=j+1
10013 if((j).gt.(maxp))goto 10012
indat(j)=0
write(6,10020)j,(labels(k,j),k=1,10)
10020 format(i3, 5h ....,10a4)
goto 10011
10012 continue
10031 continue
write(6,10040)
10040 format( 31h how many response categories k)
read(5,*)mcat
if(mcat .gt. 15)goto 10061
goto 10032
10061 continue
write(6,10070)
10070 format( 19h no more than $mcat)
goto 10031
10032 continue
write(6,10080)
10080 format( 47h is the input grouped( ie k response variables),/,
* 39h or ungrouped (ie y=1 or 2 or ... or k),/, 31h yes = grou
*ped no = ungrouped)
read(5,10090)group
10090 format(a1)
j=1
goto 10103
10101 j=j+1
10103 if((j).gt.(maxp))goto 10102
indat(j)=0
write(6,10110)j,(labels(k,j),k=1,10)
10110 format(i3, 5h ....,10a4)
goto 10101
10102 continue
if(group .ne. yes)goto 10131
write(6,10140)
10140 format( 33h give k response variable indices)
read(5,*)(indrm(j),j=1,mcat)
do 10151 j=1,mcat
indat(indrm(j))=2
10151 continue
10152 continue
indni=0
i=1
goto 10163
10161 i=i+1
10163 if((i).gt.(n))goto 10162
ni(i)=0
j=1
goto 10173
10171 j=j+1
10173 if((j).gt.(mcat))goto 10172
ym(i,j)=datam(i,indrm(j))
ni(i)=ni(i)+ym(i,j)
goto 10171
10172 continue
j=1
goto 10183
10181 j=j+1
10183 if((j).gt.(mcat))goto 10182
ym(i,j)=ym(i,j)/ni(i)
goto 10181
10182 continue
w(i)=1.0
goto 10161
10162 continue
goto 10191
10131 continue
write(6,10200)
10200 format( 35h give index of categorical response)
read(5,*)ir
indat(ir)=2
indni=0
ymin=100
ymax=0
i=1
goto 10213
10211 i=i+1
10213 if((i).gt.(n))goto 10212
yy=datam(i,ir)
if(yy .ge. ymin)goto 10231
ymin=yy
10231 continue
if(yy .le. ymax)goto 10251
ymax=yy
10251 continue
goto 10211
10212 continue
icat=int(ymax-ymin+1)
if(icat .eq. mcat)goto 10271
write(6,10280)mcat,icat
10280 format( 11h there are ,i3, 13h categories, , 4hnot ,i3)
10271 continue
if(icat .ge. mcat)goto 10301
mcat=icat
10301 continue
imin=int(ymin)
ito=imin+mcat-1
write(6,10310)mcat,imin,ito
10310 format(i3, 22h categories are used: ,i3, 4h to ,i3)
i=1
goto 10323
10321 i=i+1
10323 if((i).gt.(n))goto 10322
ni(i)=1
w(i)=1
do 10331 j=1,mcat
ym(i,j)=0.0
10331 continue
10332 continue
iy=int(datam(i,ir))-imin+1
if(iy .le. mcat)goto 10351
iy=mcat
10351 continue
ym(i,iy)=1
goto 10321
10322 continue
10191 continue
10121 continue
p=0
return
end
subroutine getw(mcat,ni,wm1,pm,pim)
real wm1(15,mcat),pm(mcat),pim(mcat)
real temp(15)
real ni
k=1
goto 10363
10361 k=k+1
10363 if((k).gt.(mcat-1))goto 10362
temp(k)=pim(k)*(1.0-pim(k))
kk=k
goto 10373
10371 kk=kk+1
10373 if((kk).gt.(mcat-1))goto 10372
wm1(k,kk)=0.0
wm1(kk,k)=0.0
goto 10371
10372 continue
goto 10361
10362 continue
k=1
goto 10383
10381 k=k+1
10383 if((k).gt.(mcat-1))goto 10382
wm1(k,k)=1.0/pm(k) + 1.0/pm(k+1)
wm1(k,k)=ni*wm1(k,k)*temp(k)*temp(k)
wm1(k,k+1)=ni*(-1.0/pm(k+1))*temp(k)*temp(k+1)
wm1(k+1,k)=wm1(k,k+1)
goto 10381
10382 continue
return
end
subroutine saxpy(n,sa,sx,incx,sy,incy)
c
c constant times a vector plus a vector.
c uses unrolled loop for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),sy(1),sa
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (sa .eq. 0.0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
sy(iy) = sy(iy) + sa*sx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
sy(i) = sy(i) + sa*sx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
sy(i) = sy(i) + sa*sx(i)
sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
50 continue
return
end
subroutine sbart(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,
&ispar,lspar,uspar,tol,isetup,xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,
&abd,p1ip,p2ip,ld4,ldnk,ier)
real xs(n),ys(n),ws(n),knot(nk+4),coef(nk),sz(n),lev(n),crit,spar,
&lspar,uspar,tol,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)
integer n,nk,isetup,icrit,ispar,ld4,ldnk,ier
realt1,t2,ratio, a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, fu,fv,fw,
&fx,x,ax,bx
integer i
i=1
23000 if(.not.(i.le.n))goto 23002
if(.not.(ws(i).gt.0))goto 23003
ws(i)=sqrt(ws(i))
23003 continue
i=i+1
goto 23000
23002 continue
if(.not.(isetup.eq.0))goto 23005
call sgram(sg0,sg1,sg2,sg3,knot,nk)
call stxwx(xs,ys,ws,n,knot,nk,xwy,hs0,hs1,hs2,hs3)
t1=0.
t2=0.
do 23007 i=3,nk-3
t1 = t1 + hs0(i)
23007 continue
do 23009 i=3,nk-3
t2 = t2 + sg0(i)
23009 continue
ratio = t1/t2
isetup = 1
23005 continue
if(.not.(ispar.eq.1))goto 23011
call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
&xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
return
23011 continue
ax=lspar
bx=uspar
c = 0.5*(3. - sqrt(5.0))
eps = 1.00
10 eps = eps/2.00
tol1 = 1.0 + eps
if(.not.(tol1 .gt. 1.00))goto 23013
go to 10
23013 continue
eps = sqrt(eps)
a = ax
b = bx
v = a + c*(b - a)
w = v
x = v
e = 0.0
spar = x
call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
&xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
fx = crit
fv = fx
fw = fx
20 xm = 0.5*(a + b)
tol1 = eps*abs(x) + tol/3.0
tol2 = 2.0*tol1
if(.not.(abs(x - xm) .le. (tol2 - 0.5*(b - a))))goto 23015
go to 90
23015 continue
if(.not.(abs(e) .le. tol1))goto 23017
go to 40
23017 continue
r = (x - w)*(fx - fv)
q = (x - v)*(fx - fw)
p = (x - v)*q - (x - w)*r
q = 2.00*(q - r)
if(.not.(q .gt. 0.0))goto 23019
p = -p
23019 continue
q = abs(q)
r = e
e = d
30 if(.not.(abs(p) .ge. abs(0.5*q*r)))goto 23021
go to 40
23021 continue
if(.not.(p .le. q*(a - x)))goto 23023
go to 40
23023 continue
if(.not.(p .ge. q*(b - x)))goto 23025
go to 40
23025 continue
d = p/q
u = x + d
if(.not.((u - a) .lt. tol2))goto 23027
d = sign(tol1, xm - x)
23027 continue
if(.not.((b - u) .lt. tol2))goto 23029
d = sign(tol1, xm - x)
23029 continue
go to 50
40 if(.not.(x .ge. xm))goto 23031
e = a - x
23031 continue
if(.not.(x .lt. xm))goto 23033
e = b - x
23033 continue
d = c*e
50 if(.not.(abs(d) .ge. tol1))goto 23035
u = x + d
23035 continue
if(.not.(abs(d) .lt. tol1))goto 23037
u = x + sign(tol1, d)
23037 continue
spar = u
call sslvrg(xs,ys,ws,n,knot,nk,coef,sz,lev,crit,icrit,spar,ratio,
&xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,ier)
fu = crit
if(.not.(fu .gt. fx))goto 23039
go to 60
23039 continue
if(.not.(u .ge. x))goto 23041
a = x
23041 continue
if(.not.(u .lt. x))goto 23043
b = x
23043 continue
v = w
fv = fw
w = x
fw = fx
x = u
fx = fu
go to 20
60 if(.not.(u .lt. x))goto 23045
a = u
23045 continue
if(.not.(u .ge. x))goto 23047
b = u
23047 continue
if(.not.(fu .le. fw))goto 23049
go to 70
23049 continue
if(.not.(w .eq. x))goto 23051
go to 70
23051 continue
if(.not.(fu .le. fv))goto 23053
go to 80
23053 continue
if(.not.(v .eq. x))goto 23055
go to 80
23055 continue
if(.not.(v .eq. w))goto 23057
go to 80
23057 continue
go to 20
70 v = w
fv = fw
w = u
fw = fu
go to 20
80 v = u
fv = fu
go to 20
90 continue
spar = x
crit = fx
return
23012 continue
return
end
Caveat receptor. [Jack] dongarra@anl-mcs, [Eric Grosse] research!ehg
Compliments of netlib Fri May 2 13:34:23 EDT 1986
real function sdot(n,sx,incx,sy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
real sx(1),sy(1),stemp
integer i,incx,incy,ix,iy,m,mp1,n
c
stemp = 0.0e0
sdot = 0.0e0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
stemp = stemp + sx(ix)*sy(iy)
ix = ix + incx
iy = iy + incy
10 continue
sdot = stemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
stemp = stemp + sx(i)*sy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
stemp = stemp + sx(i)*sy(i) + sx(i + 1)*sy(i + 1) +
* sx(i + 2)*sy(i + 2) + sx(i + 3)*sy(i + 3) + sx(i + 4)*sy(i + 4)
50 continue
60 sdot = stemp
return
end
subroutine sgram(sg0,sg1,sg2,sg3,tb,nb)
real sg0(nb),sg1(nb),sg2(nb),sg3(nb),tb(nb+4),vnikx(4,3),work(16),
&yw1(4),yw2(4),wpt
integer nb,ileft,ilo,mflag,i,ii,jj
do 23000 i=1,nb
sg0(i)=0.
sg1(i)=0.
sg2(i)=0.
sg3(i)=0.
23000 continue
ilo = 1
do 23002 i=1,nb
call interv(tb(1),(nb+1),tb(i),ileft,mflag)
call bsplvd (tb,4,tb(i),ileft,work,vnikx,3)
do 23004 ii=1,4
yw1(ii) = vnikx(ii,3)
23004 continue
call bsplvd (tb,4,tb(i+1),ileft,work,vnikx,3)
do 23006 ii=1,4
yw2(ii) = vnikx(ii,3) - yw1(ii)
23006 continue
wpt = tb(i+1) - tb(i)
if(.not.(ileft.ge.4))goto 23008
do 23010 ii=1,4
jj=ii
sg0(ileft-4+ii) = sg0(ileft-4+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
jj=ii+1
if(.not.(jj.le.4))goto 23012
sg1(ileft+ii-4) = sg1(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23012 continue
jj=ii+2
if(.not.(jj.le.4))goto 23014
sg2(ileft+ii-4) = sg2(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23014 continue
jj=ii+3
if(.not.(jj.le.4))goto 23016
sg3(ileft+ii-4) = sg3(ileft+ii-4) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23016 continue
23010 continue
goto 23009
23008 continue
if(.not.(ileft.eq.3))goto 23018
do 23020 ii=1,3
jj=ii
sg0(ileft-3+ii) = sg0(ileft-3+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
jj=ii+1
if(.not.(jj.le.3))goto 23022
sg1(ileft+ii-3) = sg1(ileft+ii-3) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23022 continue
jj=ii+2
if(.not.(jj.le.3))goto 23024
sg2(ileft+ii-3) = sg2(ileft+ii-3) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23024 continue
23020 continue
goto 23019
23018 continue
if(.not.(ileft.eq.2))goto 23026
do 23028 ii=1,2
jj=ii
sg0(ileft-2+ii) = sg0(ileft-2+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
jj=ii+1
if(.not.(jj.le.2))goto 23030
sg1(ileft+ii-2) = sg1(ileft+ii-2) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23030 continue
23028 continue
goto 23027
23026 continue
if(.not.(ileft.eq.1))goto 23032
do 23034 ii=1,1
jj=ii
sg0(ileft-1+ii) = sg0(ileft-1+ii) +wpt* (yw1(ii)*yw1(jj) + (yw2(
&ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 +yw2(ii)*yw2(jj)*.3330 )
23034 continue
23032 continue
23027 continue
23019 continue
23009 continue
23002 continue
return
end
subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag)
realabd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk),wjm3(3),wjm2(2),wjm1(1)
&,c0,c1,c2,c3
integerflag,ld4,nk,ldnk,i,j,k
wjm3(1)=0.
wjm3(2)=0.
wjm3(1)=0.
wjm2(1)=0.
wjm2(2)=0.
wjm1(1)=0.
do 23000 i=1,nk
j=nk-i+1
c0 = 1./abd(4,j)
if(.not.(j.le.nk-3))goto 23002
c1 = abd(1,j+3)*c0
c2 = abd(2,j+2)*c0
c3 = abd(3,j+1)*c0
goto 23003
23002 continue
if(.not.(j.eq.nk-2))goto 23004
c1 = 0.
c2 = abd(2,j+2)*c0
c3 = abd(3,j+1)*c0
goto 23005
23004 continue
if(.not.(j.eq.nk-1))goto 23006
c1 = 0.
c2 = 0.
c3 = abd(3,j+1)*c0
goto 23007
23006 continue
if(.not.(j.eq.nk))goto 23008
c1 = 0.
c2 = 0.
c3 = 0.
23008 continue
23007 continue
23005 continue
23003 continue
p1ip(1,j) = 0. - (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3))
p1ip(2,j) = 0. - (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2))
p1ip(3,j) = 0. - (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1))
p1ip(4,j) = c0**2 +c1**2*wjm3(1)+2.*c1*c2*wjm3(2)+2.*c1*c3*wjm3(3)
& +c2**2*wjm2(1)+2.*c2*c3*wjm2(2) +c3**2*wjm1(1)
wjm3(1)=wjm2(1)
wjm3(2)=wjm2(2)
wjm3(3)=p1ip(2,j)
wjm2(1)=wjm1(1)
wjm2(2)=p1ip(3,j)
wjm1(1)=p1ip(4,j)
23000 continue
if(.not.(flag.eq.0))goto 23010
return
23010 continue
do 23012 i=1,nk
j=nk-i+1
k=1
23014 if(.not.(k.le.4.and.j+k-1.le.nk))goto 23016
p2ip(j,j+k-1) = p1ip(5-k,j)
k=k+1
goto 23014
23016 continue
23012 continue
do 23017 i=1,nk
j=nk-i+1
k=j-4
23019 if(.not.(k.ge.1))goto 23021
c0 = 1./abd(4,k)
c1 = abd(1,k+3)*c0
c2 = abd(2,k+2)*c0
c3 = abd(3,k+1)*c0
p2ip(k,j) = 0. - ( c1*p2ip(k+3,j) +c2*p2ip(k+2,j) +c3*p2ip(k+1,j)
&)
k=k-1
goto 23019
23021 continue
23017 continue
return
23011 continue
end
subroutine sknotl(x,n,knot,k)
real x(n),knot(n+6),a1,a2,a3,a4
integer n,k,ndk,j
a1 = log(50.)/log(2.)
a2 = log(100.)/log(2.)
a3 = log(140.)/log(2.)
a4 = log(200.)/log(2.)
if(.not.(n.lt.50))goto 23000
ndk = n
goto 23001
23000 continue
if(.not.(n.ge.50 .and. n.lt.200))goto 23002
ndk = 2.**(a1+(a2-a1)*(n-50.)/150.)
goto 23003
23002 continue
if(.not.(n.ge.200 .and. n.lt.800))goto 23004
ndk = 2.**(a2+(a3-a2)*(n-200.)/600.)
goto 23005
23004 continue
if(.not.(n.ge.800 .and. n.lt.3200))goto 23006
ndk = 2.**(a3+(a4-a3)*(n-800.)/2400.)
goto 23007
23006 continue
if(.not.(n.ge.3200))goto 23008
ndk = 200. + (n-3200)**.2
23008 continue
23007 continue
23005 continue
23003 continue
23001 continue
k = ndk + 6
do 23010 j=1,3
knot(j) = x(1)
23010 continue
do 23012 j=1,ndk
knot(j+3) = x( 1 + (j-1)*(n-1)/(ndk-1) )
23012 continue
do 23014 j=1,3
knot(ndk+3+j) = x(n)
23014 continue
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine smooth(x,y,w,itag,span,dof,n,nef, cross,smo,s0,rss,quie
*t,scrat)
dimension x(n),y(n),w(n),smo(n),scrat(n),itag(n)
logical cross,quiet
real cvrss(6),cvspan(6),cvmin
integer idmin
data cvspan/0.3,0.4,0.5,0.6,0.7,1.0/
if(span .ge. 0.0)goto 10021
call splsm(x,y,w,itag,span,cross,dof, nef,n,smo,s0,rss,quiet,scra
*t)
goto 10031
10021 continue
penal=0.01
cvmin=1e15
idmin=1
if(.not.(cross))goto 10051
k=1
goto 10063
10061 k=k+1
10063 if((k).gt.(6))goto 10062
call smth(x,y,w,itag,cvspan(k), dof,n,nef,.true.,smo,s0,cvrss(k),s
*crat)
if(cvrss(k) .gt. cvmin)goto 10081
cvmin=cvrss(k)
idmin=k
10081 continue
goto 10061
10062 continue
span=cvspan(idmin)
if(penal .le. 0.)goto 10101
cvmin=(1.+penal)*cvmin
k=6
goto 10113
10111 k=k+(-1)
10113 if((-1)*((k)-(1)).gt.0)goto 10112
if(cvrss(k) .gt. cvmin)goto 10131
goto 10112
10131 continue
goto 10111
10112 continue
span=cvspan(k)
10101 continue
if(.not.(.not.quiet))goto 10151
do 10161 k=1,6
write(6,10170)k,cvspan(k),cvrss(k)
10170 format (i4, 7h span= ,f4.1, 8h cvrss= , f10.3)
10161 continue
10162 continue
write(6,10180)span
10180 format ( 20h span chosen .......,f4.1)
10151 continue
10051 continue
call smth(x,y,w,itag,span,dof,n,nef,.false.,smo,s0,rss,scrat)
10031 continue
10011 continue
return
end
subroutine smth(x,y,w,itag,span,dof,n,nef,cross,smo,s0,rss,scrat)
dimension x(n),y(n),w(n),smo(n),scrat(n),itag(n)
real scrat2(1000)
real*8 sumw,xbar,ybar,cov,var,wtot,wspan,lspan,rspan
integer left,right
logical cross,line
wtot=0
do 10191 i=1,nef
wtot=wtot+w(itag(i))
10191 continue
10192 continue
line=.true.
if(span .ge. 1.0)goto 10211
line=.false.
10211 continue
istart=1
cov=0.
var=0.
10221 if(w(itag(istart)).gt.0.0) goto 10222
istart=istart+1
goto 10221
10222 continue
xbar=x(itag(istart))
ybar=y(itag(istart))
sumw=w(itag(istart))
if(.not.(line))goto 10241
istart=istart+1
do 10251 i=istart,nef
xin=x(itag(i))
yin=y(itag(i))
win=w(itag(i))
xbar=(sumw*xbar+xin*win)/(sumw+win)
ybar=(sumw*ybar+yin*win)/(sumw+win)
cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
var=var+win*(xin-xbar)**2*(sumw+win)/sumw
sumw=sumw+win
10251 continue
10252 continue
i=1
goto 10263
10261 i=i+1
10263 if((i).gt.(nef))goto 10262
if(.not.(cross))goto 10281
xout=x(itag(i))
yout=y(itag(i))
win=w(itag(i))
cov=cov-win*(xout-xbar)*(yout-ybar)*sumw/(sumw-win)
var=var-win*(xout-xbar)**2*sumw/(sumw-win)
xbar=(sumw*xbar-win*xout)/(sumw-win)
ybar=(sumw*ybar-win*yout)/(sumw-win)
sumw=sumw-win
10281 continue
if(var .le. 0.)goto 10301
smo(itag(i))=cov*x(itag(i))/var
goto 10311
10301 continue
smo(itag(i))=0
10311 continue
10291 continue
if(.not.(cross))goto 10331
xin=x(itag(i))
yin=y(itag(i))
win=w(itag(i))
xbar=(sumw*xbar+xin*win)/(sumw+win)
ybar=(sumw*ybar+yin*win)/(sumw+win)
cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
var=var+win*(xin-xbar)**2*(sumw+win)/sumw
sumw=sumw+win
10331 continue
goto 10261
10262 continue
if(var .le. 0)goto 10351
s0=ybar-cov*xbar/var
goto 10361
10351 continue
s0=ybar
10361 continue
10341 continue
i=nef+1
goto 10373
10371 i=i+1
10373 if((i).gt.(n))goto 10372
if(var .le. 0.)goto 10391
smo(itag(i))=cov*xbar/var
goto 10401
10391 continue
smo(itag(i))=0
10401 continue
10381 continue
goto 10371
10372 continue
if(var .le. 0)goto 10421
scrat(1)=cov/var
goto 10431
10421 continue
scrat(1)=0.0
10431 continue
10411 continue
dof=1.0
goto 10441
10241 continue
do 10451 i=1,nef
scrat(i)=y(i)
scrat2(i)=w(i)
10451 continue
10452 continue
if(.not.(.not.cross))goto 10471
i=0
10481 if(i.ge.nef-1) goto 10482
i=i+1
m0=i
10491 if(x(itag(i+1)).gt.x(itag(i))) goto 10492
i=i+1
if(i .lt. nef)goto 10491
10492 continue
if(i.eq.m0)goto 10481
ntie=i-m0+1
r=0.
wt=0.
do 10501 jj=m0,i
j=itag(jj)
r=r+y(j)*w(j)
wt=wt+w(j)
10501 continue
10502 continue
if(wt .le. 0.0)goto 10521
r=r/wt
10521 continue
wt=wt/(i-m0+1)
do 10531 j=m0,i
y(itag(j))=r
w(itag(j))=wt
10531 continue
10532 continue
goto 10481
10482 continue
10471 continue
left=istart
right=istart
dof=-1.0
wspan=span*wtot/2
lspan=-w(itag(istart))
rspan=w(itag(istart))
do 10541 i=istart,nef
lspan=lspan+w(itag(i))
rspan=rspan-w(itag(i))
10551 if((rspan+w(itag(right+1))) .ge. wspan .or. right .ge. nef) goto 1
*0552
right=right+1
rspan=rspan+w(itag(right))
xin=x(itag(right))
yin=y(itag(right))
win=w(itag(right))
xbar=(sumw*xbar+xin*win)/(sumw+win)
ybar=(sumw*ybar+yin*win)/(sumw+win)
cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
var=var+win*(xin-xbar)**2*(sumw+win)/sumw
sumw=sumw+win
goto 10551
10552 continue
10561 if(lspan .le. wspan .or. left .ge. i) goto 10562
xout=x(itag(left))
yout=y(itag(left))
win=w(itag(left))
cov=cov-win*(xout-xbar)*(yout-ybar)*sumw/(sumw-win)
var=var-win*(xout-xbar)**2*sumw/(sumw-win)
xbar=(sumw*xbar-win*xout)/(sumw-win)
ybar=(sumw*ybar-win*yout)/(sumw-win)
sumw=sumw-win
lspan=lspan-w(itag(left))
left=left+1
goto 10561
10562 continue
if(.not.(cross))goto 10581
xout=x(itag(i))
yout=y(itag(i))
win=w(itag(i))
cov=cov-win*(xout-xbar)*(yout-ybar)*sumw/(sumw-win)
var=var-win*(xout-xbar)**2*sumw/(sumw-win)
xbar=(sumw*xbar-win*xout)/(sumw-win)
ybar=(sumw*ybar-win*yout)/(sumw-win)
sumw=sumw-win
10581 continue
if(var .le. 0.)goto 10601
smo(itag(i))=ybar+cov*(x(itag(i))-xbar)/var
dof=dof+w(itag(i))/sumw+ (w(itag(i))*(x(itag(i))-xbar)**2)/var
goto 10611
10601 continue
smo(itag(i))=ybar
dof=dof+w(itag(i))/sumw
10611 continue
10591 continue
if(.not.(cross))goto 10631
xin=x(itag(i))
yin=y(itag(i))
win=w(itag(i))
xbar=(sumw*xbar+xin*win)/(sumw+win)
ybar=(sumw*ybar+yin*win)/(sumw+win)
cov=cov+win*(xin-xbar)*(yin-ybar)*(sumw+win)/sumw
var=var+win*(xin-xbar)**2*(sumw+win)/sumw
sumw=sumw+win
10631 continue
10541 continue
10542 continue
if(.not.(.not.cross))goto 10651
i=0
10661 if(i.ge.nef-1) goto 10662
i=i+1
m0=i
10671 if(x(itag(i+1)).gt.x(itag(i))) goto 10672
i=i+1
if(i .lt. nef)goto 10671
10672 continue
if(i.eq.m0)goto 10661
ntie=i-m0+1
r=0.
wt=0.
do 10681 jj=m0,i
j=itag(jj)
r=r+smo(j)*w(j)
wt=wt+w(j)
10681 continue
10682 continue
if(wt .le. 0.0)goto 10701
r=r/wt
10701 continue
do 10711 j=m0,i
smo(itag(j))=r
10711 continue
10712 continue
goto 10661
10662 continue
10651 continue
do 10721 i=1,nef
y(i)=scrat(i)
w(i)=scrat2(i)
10721 continue
10722 continue
totsm=0.0
ybar=0.0
sumw=0.0
do 10731 ii=1,nef
i=itag(ii)
totsm=totsm+smo(i)*w(i)
ybar=ybar+w(i)*y(i)
sumw=sumw+w(i)
10731 continue
10732 continue
totsm=totsm/sumw
do 10741 ii=1,nef
i=itag(ii)
smo(i)=smo(i)-totsm
10741 continue
10742 continue
ii=nef+1
goto 10753
10751 ii=ii+1
10753 if((ii).gt.(n))goto 10752
i=itag(ii)
smo(i)=0.0
ybar=ybar+w(i)*y(i)
sumw=sumw+w(i)
goto 10751
10752 continue
ybar=ybar/sumw
s0=ybar
10441 continue
10231 continue
rss=0.0
do 10761 i=1,n
rss=rss+(w(i)/sumw)*(y(i)-s0-smo(i))**2
10761 continue
10762 continue
return
end
subroutine spbfa(abd,lda,n,m,info)
integer lda,n,m,info
real abd(lda,1)
c
c spbfa factors a real symmetric positive definite matrix
c stored in band form.
c
c spbfa is usually called by spbco, but it can be called
c directly with a saving in time if rcond is not needed.
c
c on entry
c
c abd real(lda, n)
c the matrix to be factored. the columns of the upper
c triangle are stored in the columns of abd and the
c diagonals of the upper triangle are stored in the
c rows of abd . see the comments below for details.
c
c lda integer
c the leading dimension of the array abd .
c lda must be .ge. m + 1 .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c 0 .le. m .lt. n .
c
c on return
c
c abd an upper triangular matrix r , stored in band
c form, so that a = trans(r)*r .
c
c info integer
c = 0 for normal return.
c = k if the leading minor of order k is not
c positive definite.
c
c band storage
c
c if a is a symmetric positive definite band matrix,
c the following program segment will set up the input.
c
c m = (band width above diagonal)
c do 20 j = 1, n
c i1 = max0(1, j-m)
c do 10 i = i1, j
c k = i-j+m+1
c abd(k,j) = a(i,j)
c 10 continue
c 20 continue
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas sdot
c fortran max0,sqrt
c
c internal variables
c
real sdot,t
real s
integer ik,j,jk,k,mu
c begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0e0
ik = m + 1
jk = max0(j-m,1)
mu = max0(m+2-j,1)
if (m .lt. mu) go to 20
do 10 k = mu, m
t = abd(k,j) - sdot(k-mu,abd(ik,jk),1,abd(mu,j),1)
t = t/abd(m+1,jk)
abd(k,j) = t
s = s + t*t
ik = ik - 1
jk = jk + 1
10 continue
20 continue
s = abd(m+1,j) - s
c ......exit
if (s .le. 0.0e0) go to 40
abd(m+1,j) = sqrt(s)
30 continue
info = 0
40 continue
return
end
subroutine spbsl(abd,lda,n,m,b)
integer lda,n,m
real abd(lda,1),b(1)
c
c spbsl solves the real symmetric positive definite band
c system a*x = b
c using the factors computed by spbco or spbfa.
c
c on entry
c
c abd real(lda, n)
c the output from spbco or spbfa.
c
c lda integer
c the leading dimension of the array abd .
c
c n integer
c the order of the matrix a .
c
c m integer
c the number of diagonals above the main diagonal.
c
c b real(n)
c the right hand side vector.
c
c on return
c
c b the solution vector x .
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal. technically this indicates
c singularity but it is usually caused by improper subroutine
c arguments. it will not occur if the subroutines are called
c correctly and info .eq. 0 .
c
c to compute inverse(a) * c where c is a matrix
c with p columns
c call spbco(abd,lda,n,rcond,z,info)
c if (rcond is too small .or. info .ne. 0) go to ...
c do 10 j = 1, p
c call spbsl(abd,lda,n,c(1,j))
c 10 continue
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas saxpy,sdot
c fortran min0
c
c internal variables
c
real sdot,t
integer k,kb,la,lb,lm
c
c solve trans(r)*y = b
c
do 10 k = 1, n
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
t = sdot(lm,abd(la,k),1,b(lb),1)
b(k) = (b(k) - t)/abd(m+1,k)
10 continue
c
c solve r*x = y
c
do 20 kb = 1, n
k = n + 1 - kb
lm = min0(k-1,m)
la = m + 1 - lm
lb = k - lm
b(k) = b(k)/abd(m+1,k)
t = -b(k)
if(lm .gt. 0)call saxpy(lm,t,abd(la,k),1,b(lb),1)
20 continue
return
end
C mortran 2.79 (reserved keyword macros of 09/28/81)
subroutine splsm(x,y,w,itag,span,cross,dof,nef,n,smo, s0,rss,quie
*t,scrat)
real nx(1000),ny(1000),nw(1000),ns(1000)
integer point(1000)
logical quiet,cross
real scrat(1)
real x(1),y(1),w(1),smo(1),rss,s0
real range,min,knot(1020),coef(1020),lev(1000)
real xwy(1000),hs0(1000),hs1(1000),hs2(1000),hs3(1000), sg0(1000)
*,sg1(1000),sg2(1000),sg3(1000)
real abd(4,1000),p1ip(4,1000),p2ip(1)
integer itag(1)
i=1
j=1
eps=.1e-9
10011 if(i .gt. nef) goto 10012
sw=0.0
sy=0.0
cx=x(itag(i))
10021 if(x(itag(i))-cx .ge. eps .or. i .gt. nef) goto 10022
point(itag(i))=j
it=itag(i)
sw=sw+w(it)
sy=sy+w(it)*y(it)
i=i+1
goto 10021
10022 continue
if(sw .le. eps)goto 10041
nx(j)=cx
nw(j)=sw
ny(j)=sy/sw
j=j+1
10041 continue
goto 10011
10012 continue
nspl=j-1
range=nx(nspl)-nx(1)
min=nx(1)
i=1
goto 10053
10051 i=i+1
10053 if((i).gt.(nspl))goto 10052
nx(i)=(nx(i)- min)/range
goto 10051
10052 continue
call sknotl(nx,nspl,knot,k)
nk=k-4
spar=-span
isetup = 0
ier = 1
icrit=1
if(.not.(cross))goto 10071
ispar=0
goto 10081
10071 continue
ispar=1
10081 continue
10061 continue
call sbart(nx,ny,nw,nspl,knot,nk, coef,ns,lev,crit,icrit,spar,ispa
*r,0.,1.,.05,isetup, xwy, hs0,hs1,hs2,hs3, sg0,sg1,sg2,sg3, abd,p1i
*p,p2ip,4,nk,ier)
span=-spar
dof=-1.0
do 10091 i=1,nspl
dof=dof+lev(i)
10091 continue
10092 continue
do 10101 ii=1,nef
i=itag(ii)
smo(i)=ns(point(i))
10101 continue
10102 continue
totsm=0.0
ybar=0.0
sumw=0.0
do 10111 ii=1,nef
i=itag(ii)
totsm=totsm+smo(i)*w(i)
ybar=ybar+w(i)*y(i)
sumw=sumw+w(i)
10111 continue
10112 continue
totsm=totsm/sumw
do 10121 ii=1,nef
i=itag(ii)
smo(i)=smo(i)-totsm
10121 continue
10122 continue
ii=nef+1
goto 10133
10131 ii=ii+1
10133 if((ii).gt.(n))goto 10132
i=itag(ii)
smo(i)=0.0
ybar=ybar+w(i)*y(i)
sumw=sumw+w(i)
goto 10131
10132 continue
ybar=ybar/sumw
s0=ybar
rss=0.0
do 10141 i=1,n
rss=rss+(w(i)/sumw)*(y(i)-s0-smo(i))**2
10141 continue
10142 continue
return
end
subroutine sslvrg(x,y,w,n,knot,nk,coef,sz,lev,crit,icrit,spar,
&ratio,xwy,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,p2ip,ld4,ldnk,
&info)
real x(n),y(n),w(n),knot(nk+4),coef(nk),sz(n),lev(n),crit,ratio,
&spar,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),lambda,b0,b1,
&b2,b3,eps,vnikx(4,1),work(16),xv,bvalue,rss,df
integer n,nk,icrit,ld4,ldnk,i,icoef,ileft,ilo,info,j,mflag
ilo = 1
eps = .1e-10
lambda = ratio*16.**(-2. + spar*(6.))
do 23000 i=1,nk
coef(i) = xwy(i)
23000 continue
do 23002 i=1,nk
abd(4,i) = hs0(i)+lambda*sg0(i)
23002 continue
do 23004 i=1,(nk-1)
abd(3,i+1) = hs1(i)+lambda*sg1(i)
23004 continue
do 23006 i=1,(nk-2)
abd(2,i+2) = hs2(i)+lambda*sg2(i)
23006 continue
do 23008 i=1,(nk-3)
abd(1,i+3) = hs3(i)+lambda*sg3(i)
23008 continue
call spbfa(abd,ld4,nk,3,info)
if(.not.(info.ne.0))goto 23010
return
23010 continue
call spbsl(abd,ld4,nk,3,coef)
icoef = 1
do 23012 i=1,n
xv = x(i)
sz(i) = bvalue(knot,coef,nk,4,xv,0)
23012 continue
if(.not.(icrit.eq.0))goto 23014
return
23014 continue
call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0)
do 23016 i=1,n
xv = x(i)
call interv(knot(1),(nk+1),xv,ileft,mflag)
if(.not.(mflag.eq.-1))goto 23018
ileft = 4
xv = knot(4)+eps
23018 continue
if(.not.(mflag.eq.1))goto 23020
ileft = nk
xv = knot(nk+1)-eps
23020 continue
j=ileft-3
call bsplvd(knot,4,xv,ileft,work,vnikx,1)
b0=vnikx(1,1)
b1=vnikx(2,1)
b2=vnikx(3,1)
b3=vnikx(4,1)
lev(i) = (p1ip(4,j)*b0**2 + 2.*p1ip(3,j)*b0*b1 +2.*p1ip(2,j)*b0*
&b2 + 2.*p1ip(1,j)*b0*b3 +p1ip(4,j+1)*b1**2 + 2.*p1ip(3,j+1)*b1*b2
&+2.*p1ip(2,j+1)*b1*b3 +p1ip(4,j+2)*b2**2 + 2.*p1ip(3,j+2)*b2*b3 +
&p1ip(4,j+3)*b3**2 )*w(i)**2
23016 continue
if(.not.(icrit.eq.1))goto 23022
rss = 0.
df = 0.
do 23024 i=1,n
rss = rss + ((y(i)-sz(i))*w(i))**2
23024 continue
do 23026 i=1,n
df = df + 1.-lev(i)
23026 continue
crit = (rss/n)/((df/n)**2)
goto 23023
23022 continue
crit = 0.
do 23028 i=1,n
crit = crit +(((y(i)-sz(i))*w(i))/(1-lev(i)))**2
23028 continue
23023 continue
return
23015 continue
end
subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3)
real z(k),w(k),x(k),xknot(n+4),y(n),hs0(n),hs1(n),hs2(n),hs3(n),
&eps,vnikx(4,1),work(16)
integer k,n,j,i,ilo,ileft,mflag
do 23000 i=1,n
y(i)=0.
hs0(i)=0.
hs1(i)=0.
hs2(i)=0.
hs3(i)=0.
23000 continue
ilo=1
eps = .1e-9
do 23002 i=1,k
call interv(xknot(1),(n+1),x(i),ileft,mflag)
if(.not.(mflag.eq. 1))goto 23004
if(.not.(x(i).le.(xknot(ileft)+eps)))goto 23006
ileft=ileft-1
goto 23007
23006 continue
return
23007 continue
23004 continue
call bsplvd (xknot,4,x(i),ileft,work,vnikx,1)
j= ileft-4+1
y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1)
hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2
hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1)
hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1)
hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1)
j= ileft-4+2
y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1)
hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2
hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1)
hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1)
j= ileft-4+3
y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1)
hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2
hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1)
j= ileft-4+4
y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1)
hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2
23002 continue
return
end