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 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