# to unbundle, sh this file (in an empty directory) mkdir bart echo bart/bart.d 1>&2 sed >bart/bart.d <<'//GO.SYSIN DD bart/bart.d' 's/^-//' -.BG -.FN bart -.TL -bart: A b-spline adaptive regression technique. -.CS -bart(x, y, w, spar, icrit) -.PP -.AG x -values of the explanatory variable. -There should be at least ten distinct x values. -.AG y -response variable- -the same length as x. -.AG w -optional -vector of weights for weighted regression. -Should have the same length as x and y. -If measurements at different x's have different -variances, w should be inversely proportional -to the variance. -.AG spar -smoothing parameter usually in the range (0,1). -If spar is missing a default amount of smoothing -is choosen by ordinary or generalized cross validation. -.AG icrit -indicates whether ordinary or generalized cross validation score -should be computed. -If icrit is missing the generalized cross validation score -is used. -.RT -The smoothed regression estimates -and -corresponding leverage values -along with the cross validation -score. -The B-spline coefficients and knot vector are also returned. -A cubic b-spline is fit -with -care taken to insure that the -algorithm runs -linear in the number of data points. -The returned structure has the -following components: -.RC tx,ty -ordered x values and corresponding smoothed y values -.RC lev -leverage values - -diagonal elements of the -equivalent hat matrix. -The degrees of freedom for signal is estimated by -the trace of the hat matrix. -.RC crit -cross validation score -.RC coef -coeficients of the B-spline. -.RC tknot -knot vector defining the B-spline. -.bp -.SH REFERENCES -.sp 1 -1. C. DeBoor, "A practical guide to B-splines", -Applied Math Sci 27, Springer, N.Y. 1978. -.sp 1 -2. P. Craven and G. Wahba, "Smoothing noisy data with spline functions: -estimating the correct degree of smoothing -by the method of generalized cross validation", -Numer. Math. 1979, vol 31, pages 377-403. -.sp 1 -3. B. W. Silverman, -"Some aspects of the spline smoothing approach to -nonparametric regression curve fitting", -JRSS Series B 1985, vol 47. -.sp 1 -4. F. O'Sullivan, -Contribution to the discussion of Silverman's paper, -JRSS Series B 1985, vol 47. -.sp 1 -.SH AUTHOR -Finbarr O'Sullivan -.sp 0 -Dept. of Statistics, UCB. -.sp 0 -finbarr@bach.BERKELEY -.sp 3 -.EX - -# Bart fit and Approximate 95% "Confidence" Intervals - - fit _ bart(x,y) # ols bart fit - res _ (y - fit$ty)/(1-fit$lev) # jacknife residuals - sigma _ sqrt(var(res)) # estimated sd - -upper _ fit$ty + 2.0*sigma*sqrt(fit$lev) # upper 95% confidence band -lower _ fit$ty - 2.0*sigma*sqrt(fit$lev) # lower 95% confidence band - -# Smooth and 95% "Confidence" bands - -matplot(fit$x,cbind(upper,fit$y,lower),type="plp",pch=".") - -.KW ~smoothing -.KW ~b-spline -.KW ~cross validation -.WR //GO.SYSIN DD bart/bart.d echo bart/bart.i 1>&2 sed >bart/bart.i <<'//GO.SYSIN DD bart/bart.i' 's/^-//' -FUNCTION bart ( x /REAL/ - y /REAL/ - w /REAL,LENGTH(x)/ - spar /REAL,1/ - icrit /INT,1/ - ) - -STATIC( integer i,ier,ispar,isetup,n,nx,nk) - -if(LENGTH(x)!=LENGTH(y)) FATAL(lengths of x and y must match) -if(LENGTH(x)!=LENGTH(w)) FATAL(lengths of x and w must match) -if(MISSING(w)) {do i=1,LENGTH(x) { w[i] = 1.0 } } -if(MISSING(spar)) { ispar = 0 } -else { ispar = 1 } -if(MISSING(icrit)) { icrit = 1 } -n = LENGTH(x) -isetup = 0 -ier = 1 - -STRUCTURE (knot/REAL,n+6/,xw/REAL,n/) -STRUCTURE (min/REAL,1/,range/REAL,1/) - - call setreg(x,y,w,LENGTH(x),xw,nx,min,range,knot,nk) - - -STRUCTURE ( xwy/REAL,nk/, - hs0/REAL,nk/,hs1/LIKE(hs0)/,hs2/LIKE(hs0)/,hs3/LIKE(hs0)/, - sg0/REAL,nk/,sg1/LIKE(hs0)/,sg2/LIKE(hs0)/,sg3/LIKE(hs0)/, - abd/MATRIX,4,nk/,p1ip/MATRIX,4,nk/,p2ip/MATRIX,nk,nk/) - -STRUCTURE ( tx/REAL,nx/ - ty/REAL,nx/ - lev/REAL,nx/ - crit/REAL,1/ - coef/REAL,nk/,tknot/REAL,nk+4/) - - call sbart(x,y,w,nx,knot,nk, - coef,ty,lev, - crit,icrit,spar,ispar,0.,1.,.05, - isetup, - xwy, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,4,nk,ier) - -if(ier>0) FATAL ( smoothing parameter value too small or too large) - -do i=1,nx { tx[i] = min + range*x[i] } -do i=1,nk+4 { tknot[i] = min + range*knot[i] } - -RETURN(tx,ty,lev,crit,coef,tknot) -END //GO.SYSIN DD bart/bart.i echo bart/bden.d 1>&2 sed >bart/bden.d <<'//GO.SYSIN DD bart/bden.d' 's/^-//' -.BG -.FN bden -.TL -bden: A b-spline density estimator -.CS -bden(x, w, spar, icrit) -.PP -.AG x -a random sample -.AG w -optional -vector of weights attached to each data value. -Higher weights give the corresponding data values more -influence on the estimator. -.AG spar -smoothing parameter usually in the range (0,1). -If spar is missing a default amount of smoothing -is choosen by a linearized likelihood (Habemma Hermans and van den Brock) -or quadratic (Rudemo) -cross-validation. -.AG icrit -indicates whether likelihood (icrit=1) or -quadratic (icrit=2) cross-validation score -should be computed. -If icrit is missing the likelihood cross-validation score -is used. -.RT -The smooth density, -the cumulative, -along with -the B-spline coefficients and knot vector -for the logarithm of the -density estimate are returned. -A penalized maximum likelihood estimator for the logarithm -of the density is used. -The computational algorithm is based on a -Newton-Raphson iteration combined with a Marquart-Levenberg line -search. -The method -for computing the cross-validation scores is -linear in the number of data points. -The returned structure has the -following components: -.RC ox -ordered x values -.RC den -smoothed density estimates evaluated at ox -.RC init -initial histogram estimate evaluated at ox -.RC lev -leverage values at ox, these can be used to construct Bayesian Confidence -bands -.RC scdf -corresponding cumulative distribution function estimate -.RC ecdf -empirical cumulative distribution function -.RC crit -value of the cross-validation score -.RC coef -B-spline coefficients for the logarithmic density. -.RC knot -Knot points for the B-spline representation. -.bp -.SH REFERENCES -.sp 1 -1. C. DeBoor, "A practical guide to B-splines", -Applied Math Sci 27, Springer, N.Y. 1978. -.sp 1 -2. F. O'Sullivan, -Penalized likelihood -density and hazard estimation", -Dept. of Statistics, UC-Berkeley, Tech. Rep. January 1986. -.sp 1 -3. B. W. Silverman, -"On the estimation of a probability density function by the -maximum penalized likelihood method", -Annals of Statistics 1982, vol 10.,No 3,pp 795-810. -.sp 1 -.SH AUTHOR -Finbarr O'Sullivan -.sp 0 -Dept. of Statistics, UCB. -.sp 0 -finbarr@bach.BERKELEY -.sp 2 -.EX - -# Generate a sample from a beta(7,11) density - - x _ rbeta(100,7.,11.) - -# Use bden to estimate the density and construct Bayesian Confidence Intervals - - fit _ bden(x) - lo _ log(fit$den) - 2.*sqrt(fit$lev) ; lo _ exp(lo) - up _ log(fit$den) + 2.*sqrt(fit$lev) ; up _ exp(up) - -# Plot the density and the confidence intervals - - matplot(fit$ox,cbind(lo,fit$den,up),type="plp",pch=".") - -# Compare the empirical and smoothed cdf's - -matplot(fit$ox,cbind(fit$ecdf,fit$scdf),type="pl",pch="-") - - - -.KW ~smoothing -.KW ~b-spline -.KW ~cross-validation -.WR //GO.SYSIN DD bart/bden.d echo bart/bden.i 1>&2 sed >bart/bden.i <<'//GO.SYSIN DD bart/bden.i' 's/^-//' -FUNCTION bden( x /REAL/ - w /REAL,LENGTH(x)/ - a /REAL,1/ - b /REAL,1/ - spar /REAL,1/ - icrit /INT,1/ - ) - -STATIC( integer dcrit,i,ier,ispar,icoef,ind,gxi,gxio, - k,n,nx,nk,nbvcd,nbvcg,ngpt,ilo,ileft,mflag - real av,bv,xv,eps,b0,b1,b2,b3,bvalu,exp,gq1) - -if(LENGTH(x)!=LENGTH(w)) FATAL(lengths of x and w must match) -if(MISSING(w)) {do i=1,LENGTH(x) { w[i] = 1.0 } } -if(MISSING(a) .or. MISSING(b)) { ind = 0 } -else { av = a ; bv = b } -if(MISSING(spar)) { ispar = 0 } -else { ispar = 1 } -if(MISSING(icrit)) { icrit = 3 } - -n = LENGTH(x) -nk = LENGTH(x) + 4 -dcrit = 4 -ngpt = 4*(nk-3) -nbvcd = 4*n+nk -nbvcg = 4*ngpt+nk -ier = 1 -eps = .1e-10 - -STRUCTURE( wknot /REAL,nk+4/ - xw /REAL,nk/ - grid /REAL,ngpt/ - qwts /REAL,ngpt/ - bvcd /REAL,nbvcd/ - sd /INT,nk/ - ld /INT,nk/ - nd /INT,nk/ - bvcg /REAL,nbvcg/ - sg /INT,nk/ - lg /INT,nk/ - ng /INT,nk/ - ecum /REAL,n/ - cmin /REAL,1/ - cmax /REAL,1/ - coef0 /REAL,nk/ - coef1 /REAL,nk/ - hcoef /REAL,nk/ - wcoef /REAL,nk/ - smden /REAL,ngpt/ - levg /LIKE(smden)/ - wtsg /LIKE(smden)/ - vnikx /MATRIX,4,1/ - ) - - - call setden(x,w,n,av,bv,ind,xw,nx, - wknot,nk, - grid,qwts,ngpt, - bvcd,nbvcd,sd,ld,nd, - bvcg,nbvcg,sg,lg,ng, - ecum,coef0,coef1,smden) - -STRUCTURE ( hs0/REAL,nk/,hs1/LIKE(hs0)/,hs2/LIKE(hs0)/,hs3/LIKE(hs0)/, - sg0/REAL,nk/,sg1/LIKE(hs0)/,sg2/LIKE(hs0)/,sg3/LIKE(hs0)/, - cvlik/REAL,nx/,wrk1/REAL,nk/,wrk2/REAL,nk/,wrk3/REAL,nk/, - abd/MATRIX,4,nk/, - p1ip/MATRIX,4,nk/,p2ip/MATRIX,nk,nk/) - -STRUCTURE ( ox /REAL,nx/ - ecdf /LIKE(ox)/ - scdf /LIKE(ox)/ - den /LIKE(ox)/ - init /LIKE(ox)/ - lev /LIKE(ox)/ - crit /REAL,dcrit/ - coef /REAL,nk/,knot/REAL,nk+4/) - - call sbden(x,w,nx,wknot,nk, - coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - wcoef,smden,cvlik,levg,wtsg, - crit,dcrit,icrit, - spar,ispar,0.,1.,.05, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,4,nk, - wrk1,wrk2,wrk3,ier) - - - -# Compute Leverages for Bayesian Confidence Intervals - - do i=1,LENGTH(ox) { xv = x[i] - - - - call intrv(wknot[1],nk+1,xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = wknot[4]+eps } - if(mflag==1) { ileft = nk ; xv = wknot[nk+1]-eps } - j=ileft-3 - call bspvd(wknot,4,1,xv,ileft,4,vnikx,wrk1) - 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 ) - - lev[i] = 2.*w[i]*lev[i] } - - -if(ier>0) FATAL ( smoothing parameter value too small or too large) - - # Empirical and Smoothed cdf's - - - -do i=1,LENGTH(ox) { ecdf[i] = ecum[i] } - - gxi = 0 - - do i=1,LENGTH(ox) { gxio = gxi+1 - - for(k=gxio;k<=ngpt;k=k+1) - {if(grid[k]<=x[i]) {gxi=k} # gxi - else {break} } - - - if(gxio&2 sed >bart/bhaz.d <<'//GO.SYSIN DD bart/bhaz.d' 's/^-//' -.BG -.FN bhaz -.TL -bhaz: A b-spline hazard estimator -.CS -bhaz(t, m, c, w, spar) -.PP -.AG t -a sample of observed failure times. -elements of t need not be distinct. -.AG m -optional -vector indicating the number of subjects associated with -a particular t(i) value. -By default m(i) = 1. -.AG c -optional -vector indicating the proportion of -uncensored failure times at t(i) -By default no censoring is assumed. -.AG w -optional -vector of weights attached to each data value. -One could estimate the baseline hazard in a Cox -regression model by letting the -weights be exp( x(i)'beta ) - -x(i) is the vector of covariates on the i'th individual -and beta is the Cox regression parameter estimate. -.AG spar -smoothing parameter usually in the range (0,1). -If spar is missing a default amount of smoothing -is choosen by a linearized likelihood -cross-validation. -.RT -The smooth hazard, -the corresponding survival function, -the empirical survival function along with -the B-spline coefficients and knot vector -for the logarithm -of the hazard estimate. -A penalized maximum likelihood estimator for the logarithm -of the hazard is used. -The algorithm is based on a -Newton-Raphson iteration combined with a Marquart-Levenberg line -search. -The method -for computing the cross-validation score is -quadratic in the number of data points. -The returned structure has the -following components: -.RC ot -ordered distinct t values -.RC den -smoothed hazard estimates evaluated at ot -.RC ss -corresponding survival curve -.RC es -empirical survival curve -.RC crit -value of the cross-validation score -.RC coef -B-spline coefficients for the logarithmic density. -.RC knot -Knot points for the B-spline representation. -.bp -.SH REFERENCES -.sp 1 -1. J. A. Anderson and A. Senthilselvan, -"Smooth estimates for the hazard function", -J. R. Statist. Soc. B. 42,3,pp322-327. -.sp 1 -2. C. DeBoor, "A practical guide to B-splines", -Applied Math Sci 27, Springer, N.Y. 1978. -.sp 1 -3. F. O'Sullivan, -Computation of cross-validated penalized likelihood -density and hazard estimators", -Dept. of Statistics, UC-Berkeley, Tech. Rep. # 1986. -.sp 1 -.SH AUTHOR -Finbarr O'Sullivan -.sp 0 -Dept. of Statistics, UCB. -.sp 0 -finbarr@bach.BERKELEY -.sp 3 -.EX - -# Generate a sample from a Weibull(1,3) distribution -# with exponential(10.) censoring - - u _ runif(100) ; v _ runif(100) - - y _ (-log(u))**(1./3.) # Failure time - z _ (-log(v))*8.15 # Censoring Variables (leads to about 10% censoring) - - t _ ifelse(y&2 sed >bart/breslw.r <<'//GO.SYSIN DD bart/breslw.r' 's/^-//' -subroutine breslw(t,c,m,w,nt,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - -real t(nt),c(nt),m(nt),w(nt), - cumhaz(nt), - wrk0(nt), # sorted t values - wrk1(nt), # sorted t values - wrk2(nt+6), # t-sorted c values and b-spline knot vector - wrk3(nt), # t-sorted m values - wrk4(nt+2) # t-sorted w values and ecum coeficients - -integer nt,i,nd,icoef - - - # Local - -real a,b,deltl,deltu,eps,sm,bvalu - -integer i,nd,icoef,ileft,mflag - - - eps = .1e-7 - - - - - - # Sort t in increasing order { as in sethaz } - - call scopy(nt,t,1,wrk1,1) ; call scopy(nt,c,1,wrk2,1) - call ssort(wrk1,wrk2,nt,2) - call scopy(nt,t,1,wrk1,1) ; call scopy(nt,m,1,wrk3,1) - call ssort(wrk1,wrk3,nt,2) - call scopy(nt,t,1,wrk1,1) ; call scopy(nt,w,1,wrk4,1) - call ssort(wrk1,wrk4,nt,2) - - # Scale wrk1 - - deltl = (wrk1(nt/10)-wrk1(1))/(.1*nt-1.) - deltu = (wrk1(nt)-wrk1((nt/10)*9))/(.1*nt-1.) - a = wrk1(1) - deltl ; b = wrk1(nt) + deltu - if(b-a < eps**2) { write(6,'(" error in breslow estimator ")') - stop } - - do i=1,nt { wrk1(i) = (wrk1(i)-a)/(b-a) } - call scopy(nt,wrk1,1,wrk0,1) - - # Get Distinct Values - - nd = 1 ; wrk0(nd) = wrk0(1) ; wrk4(nd) = wrk4(1)*wrk3(1) - wrk2(nd) = wrk2(1) ; wrk3(nd) = wrk3(1) - - for(i=2;i<=nt;i=i+1) - { if(wrk1(i)>wrk0(nd)+eps) - { nd=nd+1 ; wrk0(nd)=wrk0(i) ; wrk3(nd)=wrk3(i) - wrk4(nd)=wrk4(i)*wrk3(i) ; wrk2(nd)=wrk2(i) } - else - { wrk3(nd)=wrk3(nd)+wrk3(i) ; wrk4(nd)=wrk4(nd)+wrk4(i)*wrk3(i) - wrk2(nd)=wrk2(nd)+wrk2(i) } - } - - - # Sum m(i) and scale w(i) and c(i) - - sm = 0. - do i=1,nd { sm = sm + wrk3(i) } - do i=1,nd { wrk4(i) = wrk4(i)/sm } - do i=1,nd { wrk2(i) = wrk2(i)/sm } - - - # Sorted t-vector is in wrk1 the corresponding - # c,m,w vectors are in wrk2,3,4 - # there are nd distinct values - - - - # Cumulative Hazard -> { wrk4 (ecum of ihaz) - # use wrk3 like yw1 } - - - wrk3(1) = wrk4(1) - do i=2,nd { wrk3(1) = wrk3(1) + wrk4(i) } - do i=2,nd { wrk3(i) = wrk3(i-1) - wrk4(i-1) } - - # Cumulative Hazard b-spline coefficients -> wrk4 - - wrk4(1) = wrk2(1)/wrk3(1) - do i=2,nd { wrk4(i) = wrk4(i-1) + wrk2(i)/wrk3(i) } - - - do i=1,nd{ wrk2(i) = a + (b-a)*(wrk0(i)) } - - icoef = 1 - do i=1,nt {call intrv(wrk2,nd,t(i),icoef,ileft,mflag) - cumhaz(i) = wrk4(ileft) } - - - # Knots into wrk2 - - #do i=1,3 { wrk2(i) = a + (b-a)*(wrk0(1)-eps) } - #do i=3,nd-1{ wrk2(i+1) = a + (b-a)*(wrk0(i-1)+wrk0(i))*.5 } - #do i=1,3 { wrk2(nd+i) = a + (b-a)*(wrk0(nd)+eps) } - - - # Evaluate estimate at t - - #icoef = 1 - #do i=1,nt { cumhaz(i) = bvalu(wrk2,wrk4,nd,3,0,t(i),icoef,wrk1) } - - -return -end //GO.SYSIN DD bart/breslw.r echo bart/cbnm.r 1>&2 sed >bart/cbnm.r <<'//GO.SYSIN DD bart/cbnm.r' 's/^-//' -subroutine cbnm(x0,x1,x2,x3,y,k,z) - -real x0(k+3),x1(k+3),x2(k+3),x3(k+3),y(k+3),z(k+3) -integer j,k - - - - j=1 - z(j) = x0(1)*y(1)+ - x1(1)*y(2)+x2(1)*y(3)+x3(1)*y(4) - - j=2 - z(j) = x1(1)*y(1)+ - x0(2)*y(2)+ - x1(2)*y(3)+x2(2)*y(4)+x3(2)*y(5) - - j=3 - z(j) = x2(1)*y(1)+x1(2)*y(2)+ - x0(3)*y(3)+ - x1(3)*y(4)+x2(3)*y(5)+x3(3)*y(6) - - do j=4,k{ - z(j) = x3(j-3)*y(j-3)+x2(j-2)*y(j-2)+x1(j-1)*y(j-1)+ - x0(j)*y(j)+ - x1(j)*y(j+1)+x2(j)*y(j+2)+x3(j)*y(j+3) } - - j=k+1 - z(j) = x3(j-3)*y(j-3)+x2(j-2)*y(j-2)+x1(j-1)*y(j-1)+ - x0(j)*y(j)+ - x1(j)*y(j+1)+x2(j)*y(j+2) - - j=k+2 - z(j) = x3(j-3)*y(j-3)+x2(j-2)*y(j-2)+x1(j-1)*y(j-1)+ - x0(j)*y(j)+ - x1(j)*y(j+1) - - j=k+3 - z(j) = x3(j-3)*y(j-3)+x2(j-2)*y(j-2)+x1(j-1)*y(j-1)+ - x0(j)*y(j) - -return -end //GO.SYSIN DD bart/cbnm.r echo bart/cxghs.r 1>&2 sed >bart/cxghs.r <<'//GO.SYSIN DD bart/cxghs.r' 's/^-//' - -subroutine cxghs(t,x,c,m,sm,w,n,iflag, - knot,grad,hs0,hs1,hs2,hs3,hess,ldnk,nk, - lambda,coef0,sg0,sg1,sg2,sg3, - sgrad,shs0,shs1,shs2,shs3) - - -real t(n),x(n),c(n),m(n),sm,w(n), - knot(nk+4), - grad(nk),hs0(nk),hs1(nk),hs2(nk),hs3(nk),hess(ldnk,nk), - lambda,coef0(nk),sg0(nk),sg1(nk),sg2(nk),sg3(nk), - sgrad(nk),shs0(nk),shs1(nk),shs2(nk),shs3(nk) - -integer n,iflag,nk,ldnk - - - - # Local Variables - -real xv,b0,b1,b2,b3,eps,vnikx(4,1),work(15),rito -integer i,j,k,l,ilo,ileft,mflag - - - - # Computes the gradient and Hessian for the Cox Model - # only if iflag==1 is the full Hessian computed - - - - - # Initializations - - ilo = 0 ; eps = .1e-10 - - rito = 0.; do j=1,nk { grad(j) = 0. } ; do j=1,nk { sgrad(j) = 0. } - do j=1,nk { hs0(j) = 0. } ; do j=1,nk { shs0(j) = 0. } - do j=1,nk { hs1(j) = 0. } ; do j=1,nk { shs1(j) = 0. } - do j=1,nk { hs2(j) = 0. } ; do j=1,nk { shs2(j) = 0. } - do j=1,nk { hs3(j) = 0. } ; do j=1,nk { shs3(j) = 0. } - - if(iflag==1) { do j=1,nk { do k=j,nk { hess(k,j) = 0. }}} - - - - - # The Algorithm is straightforward - - i=1 - - for(l=i;l<=n;l=l+1) - { rito = rito + m(l)*w(l) - - # Get the relevant b-spline coefficients - xv = x(l) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - - # Do the appropriate adjustments - - sgrad(j) = sgrad(j) + m(l)*w(l)*b0 - sgrad(j+1) = sgrad(j+1) + m(l)*w(l)*b1 - sgrad(j+2) = sgrad(j+2) + m(l)*w(l)*b2 - sgrad(j+3) = sgrad(j+3) + m(l)*w(l)*b3 - - shs0(j) = shs0(j) + m(l)*w(l)*b0**2. - shs0(j+1) = shs0(j+1) + m(l)*w(l)*b1**2. - shs0(j+2) = shs0(j+2) + m(l)*w(l)*b2**2. - shs0(j+3) = shs0(j+3) + m(l)*w(l)*b3**2. - - shs1(j) = shs1(j) + m(l)*w(l)*b0*b1 - shs1(j+1) = shs1(j+1) + m(l)*w(l)*b1*b2 - shs1(j+2) = shs1(j+2) + m(l)*w(l)*b2*b3 - - shs2(j) = shs2(j) + m(l)*w(l)*b0*b2 - shs2(j+1) = shs2(j+1) + m(l)*w(l)*b1*b3 - - shs3(j) = shs3(j) + m(l)*w(l)*b0*b3 } - - # If c(i) > 0 modify the output vectors - - if(c(i)>0.) { - do j=1,nk { grad(j) = grad(j) + c(i)*sgrad(j)/(rito*sm) } - - do j=1,nk { hs0(j) = hs0(j) + c(i)*shs0(j)/(rito*sm) } - do j=1,nk { hs1(j) = hs1(j) + c(i)*shs1(j)/(rito*sm) } - do j=1,nk { hs2(j) = hs2(j) + c(i)*shs2(j)/(rito*sm) } - do j=1,nk { hs3(j) = hs3(j) + c(i)*shs3(j)/(rito*sm) } - - # Extra Part of the Gradient - xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - - grad(j) = grad(j) - c(i)*b0/sm - grad(j+1) = grad(j+1) - c(i)*b1/sm - grad(j+2) = grad(j+2) - c(i)*b2/sm - grad(j+3) = grad(j+3) - c(i)*b3/sm - - # if iflag == 1 update Full Hessian also - - if(iflag==1) { - do j=1,nk { do k=j,nk { hess(k,j) = hess(k,j) - - c(i)*(sgrad(j)/rito)*(sgrad(k)/rito)/sm }}} - - } - - - for(i=2;i<=n;i=i+1){ - - if(t(i)>t(i-1)) - { rito = rito - m(i-1)*w(i-1) - - # Get the relevant b-spline coefficients - xv = x(i-1) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - - # Do the appropriate adjustments - - sgrad(j) = sgrad(j) - m(i-1)*w(i-1)*b0 - sgrad(j+1) = sgrad(j+1) - m(i-1)*w(i-1)*b1 - sgrad(j+2) = sgrad(j+2) - m(i-1)*w(i-1)*b2 - sgrad(j+3) = sgrad(j+3) - m(i-1)*w(i-1)*b3 - - shs0(j) = shs0(j) - m(i-1)*w(i-1)*b0**2. - shs0(j+1) = shs0(j+1) - m(i-1)*w(i-1)*b1**2. - shs0(j+2) = shs0(j+2) - m(i-1)*w(i-1)*b2**2. - shs0(j+3) = shs0(j+3) - m(i-1)*w(i-1)*b3**2. - - shs1(j) = shs1(j) - m(i-1)*w(i-1)*b0*b1 - shs1(j+1) = shs1(j+1) - m(i-1)*w(i-1)*b1*b2 - shs1(j+2) = shs1(j+2) - m(i-1)*w(i-1)*b2*b3 - - shs2(j) = shs2(j) - m(i-1)*w(i-1)*b0*b2 - shs2(j+1) = shs2(j+1) - m(i-1)*w(i-1)*b1*b3 - - shs3(j) = shs3(j) - m(i-1)*w(i-1)*b0*b3 } - - # If c(i) > 0 modify the output vectors - - if(c(i)>0.) { - do j=1,nk { grad(j) = grad(j) + c(i)*sgrad(j)/(rito*sm) } - - do j=1,nk { hs0(j) = hs0(j) + c(i)*shs0(j)/(rito*sm) } - do j=1,nk { hs1(j) = hs1(j) + c(i)*shs1(j)/(rito*sm) } - do j=1,nk { hs2(j) = hs2(j) + c(i)*shs2(j)/(rito*sm) } - do j=1,nk { hs3(j) = hs3(j) + c(i)*shs3(j)/(rito*sm) } - - # Extra Part of the Gradient - xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - - grad(j) = grad(j) - c(i)*b0/sm - grad(j+1) = grad(j+1) - c(i)*b1/sm - grad(j+2) = grad(j+2) - c(i)*b2/sm - grad(j+3) = grad(j+3) - c(i)*b3/sm - - # if iflag == 1 update Full Hessian also - - if(iflag==1) { - do j=1,nk { do k=j,nk { hess(k,j) = hess(k,j) - - c(i)*(sgrad(j)/rito)*(sgrad(k)/rito)/sm }}} - - } - - - } # i loop - - - # Contribution to the gradient form the prior - - - call cbnm(sg0,sg1,sg2,sg3,coef0,(nk-3),sgrad) - do j=1,nk { grad(j) = grad(j) + 2.*lambda*sgrad(j) } - -return -end //GO.SYSIN DD bart/cxghs.r echo bart/elabsrcx.r 1>&2 sed >bart/elabsrcx.r <<'//GO.SYSIN DD bart/elabsrcx.r' 's/^-//' -Programs developed for the Cox Model but not used. - -a) A search routine - -b) A conjugate gradient minimization algorithm - -c) An olde version of sslvcx routine - -d) gxcv gxcvo attempts at cross validation for the cox model - - - -subroutine srchcx(x,t,c,m,sm,w,n, - sx,sc,nsx,isx, - cmin,cmax,coef,coef0,dir,nk, - abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - smreg,knot, - lambda,pl,plo,ierr,wrk0) - - -real x(n),t(n),c(n),m(n),sm,w(n), - sx(nsx),sc(nsx), - cmin,cmax,coef(nk),coef0(nk),dir(nk), - abd(ld4,nk), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - smreg(n),knot(nk+4), - lambda,pl,plo, - wrk0(nk) - - -integer n,isx(n),nsx,nk,ld4,ierr - - - # Local Variables - -real work(15),bvalu,plm, - exp,diag,alpham,tinit,tdes,tau,snrm2, - sqrt,nco,abs, - x0,x1,x2,x3, - y0,y1,y2,y3, - y0b -integer i,icoef,j,k,cnt,cnt1 - - - - - # - # A trust region implementation of a - # Levenberg-Marquart Line Search for - # a new estimate is used - # - # coef = coef0 - [abd+alpha*I] (-1) dir - # - # where the values of alpha are choosen - # so that the step sizes - # - # norm[coef-coef0]/(norm[coef0]+.01) - # - # are 4**(-k) , k = 1,2, ... 5. - # - # The method of false position (with Illinois Modification) - # is used to find the correct alpha. - # - - - - - - call evalpl(x,t,c,m,sm,w,n, - cmin,cmax,coef,coef0,nk, - smreg,knot, - lambda,sg0,sg1,sg2,sg3, - wrk0,pl) - - - # Size of full step - - nco = snrm2(nk,coef0,1)/sqrt(nk*1.) - do i=1,nk { wrk0(i) = coef(i) - coef0(i) } - tinit = snrm2(nk,wrk0,1)/(sqrt(nk*1.)*(nco+1.)) - - - - # Step Length Search - - - diag = 0. ; do j=1,nk { diag = diag+hs0(j)+2.*lambda*sg0(j)} - diag = diag/nk - - alpham = 0. ; plm = pl ; tau = tinit - - x1 = 0. ; y1 = 1. - - - do k = 1, 10 {# Next Step Size - - tdes = tinit*(2.)**(-k) - - # Initialize x0 and x1 - - y1 = 1. ; cnt = 0 - while( y1 > 0. ) { x1 = 4.*(x1+.5) ; cnt = cnt + 1 - if(cnt>15) { pl = plm + 20. ; break } - call step(x1,tau,nco,diag,cmin,cmax,coef0, - coef,dir,nk,x,n,knot,abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3,lambda) - y1 = tau/tdes - 1. } - if(cnt>15) {go to 11} - x0 = x1 ; y0 = y1; cnt = 0 - while( y0 < 0. ) { x0 = x0/4. ; cnt = cnt + 1 - if(cnt>15) { pl = plm + 20. ; break } - call step(x0,tau,nco,diag,cmin,cmax,coef0, - coef,dir,nk,x,n,knot,abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3,lambda) - y0 = tau/tdes - 1. } - if(cnt>15) {go to 11} - - # x0 is smaller than the correct value for alpha ( y0 > 0. ) - # x1 is larger than the correct value for alpha ( y1 < 0. ) - - - - # The Method of False Position Starts Here - - cnt = 0 - while ( abs(y1) > .01 ) { cnt = cnt + 1 - - x2 = (y1/(y1-y0))*x0 + (y0/(y0-y1))*x1 - - call step(x2,tau,nco,diag,cmin,cmax,coef0, - coef,dir,nk,x,n,knot,abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3,lambda) - - y2 = tau/tdes - 1. - - if(cnt>7) { - if(abs(y2)<.1) { y2 = .0} - else if(abs(y1)<.1) { y2 = .0 ; x2 = x1 } - else if(abs(y0)<.1) { y2 = .0 ; x2 = x0 } - else { write(18,'(" cnt catch ",i3)') k - write(18,*) x0,y0,x2,y2,x1,y1 - call flush(18) - y2 = .0 - if(abs(y0) 0) - { y0b = .5*y0b - x3 = (y0b/(y0b-y3))*x3 + - (y3/(y3-y0b))*x0 - call step(x3,tau,nco,diag,cmin,cmax,coef0, - coef,dir,nk,x,n,knot,abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3,lambda) - - y3 = tau/tdes - 1. - - if(cnt1>7) { - if(abs(y3)<.1) { y3 = .00001} - else { write(18,'(" cnt1 catch ",i3)') k - write(18,*) x0,y0b,x2,y2,x3,y3 - call flush(18) - y3 = .0 } } - - } - - x0 = x2 ; y0 = y2 - x1 = x3 ; y1 = y3 } - - } # End of False Position Loop - - - - do j=1,nk { coef(j) = coef0(j) - coef(j) } - - call evalpl(x,t,c,m,sm,w,n, - cmin,cmax,coef,coef0,nk, - smreg,knot, - lambda,sg0,sg1,sg2,sg3, - wrk0,pl) - -11 if(pl < plm) { alpham = x1 ; plm = pl } - - } # Go onto the next step size - - - - - - # - # Choose the Best Step - # - - - pl = plm - - if(pl > plo ) { ierr = -1 ; pl = plo - do j=1,nk { coef(j) = coef0(j) } - - # Regression Function - - icoef = 1 - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk,4,0,x(i),icoef,work))* - w(i)} - - } - - - else { - - call step(alpham,tau,nco,diag,cmin,cmax,coef0, - coef,dir,nk,x,n,knot,abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3,lambda) - - - do j=1,nk { coef(j) = coef0(j) - coef(j) } - - call evalpl(x,t,c,m,sm,w,n, - cmin,cmax,coef,coef0,nk, - smreg,knot, - lambda,sg0,sg1,sg2,sg3, - wrk0,pl) - - } - - - -return -end - - - - -subroutine cgm(x,t,c,m,sm,w,n, - sx,sc,isx,nsx, - bvcd,nbvcd,sd,nd,ld, - cmin,cmax,coef0,coef,dir, - knot,nk,smreg,cumhaz, - wrk0,wrk1,wrk2,wrk3,wrk4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3,plo,pl) - - - -real x(n),t(n),c(n),m(n),sm,w(n), - sx(nsx),sc(nsx), - bvcd(nbvcd), - cmin,cmax,coef0(nk),coef(nk),dir(nk), - knot(nk+4),smreg(n),cumhaz(n), - wrk0(n+2),wrk1(n),wrk2(n+6),wrk3(n),wrk4(n+2), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - plo,pl - - -integer isx(n),nsx,n,nk,nbvcd, - sd(nk),ld(nk),nd(nk) - - - # Local Variables - -real pl,plo,exp,bvalu,cr,ng0,ng1,tol - -integer i,icoef,j,k - - - - # Initialization - - ng0 = snrm2(n,dir,1)**2. ; tol = .1e-6 - - do k=1,10 { - - # Find best estimate in the direction of dir - - call srchcx(x,t,c,m,sm,w,n,cmin,cmax,coef,coef0,dir,nk, - smreg,knot,sg0,sg1,sg2,sg3,lambda,pl,plo,wrk0) - - # Find the gradient at coef - - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - do i=1,nsx { wrk4(i) = 0. } - do i=1,n { wrk4(isx(i)) = wrk4(isx(i)) + cumhaz(i)*m(i)*w(i) } - - icoef = 1 ; do i=1,nsx { - smreg(i) = exp(bvalu(knot,coef,nk,4,0,sx(i),icoef,wrk1))} - - call sdircx(sc,wrk4,sm,nsx,nk, - bvcd,sd,ld,nd,nbvcd, - smreg,lambda,sg0,sg1,sg2,sg3, - coef,wrk0) - - ng1 = snrm2(n,wrk0,1)**2. ; cr = ng1/ng0 ; ng0 = ng1 - - if(crhcoef(j)) { cmin=hcoef(j) } } - - sm = 0. ; do i=1,n { sm = sm + m(i) } - - # Purpose : Iteratively computes a Cox Regression - # The initial guess is in coef - - lambda = ratio*16.**(-2. + spar*(6.)) - - - # Put the regression estimate into smreg - - icoef = 1 - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk, - 4,0,x(i),icoef,wrk1))*w(i)} - call evalle(t,c,m,sm,smreg,n,le) - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk0) - - plo = le + lambda*sdot(nk,coef,1,wrk0,1) - - - # Initial Cumulative Hazard Estimate - - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - -# -# Newton's Method with a -# Trust Region Implementation of a Levenberg-Marquard Line Search -# -# coef = coef0 - [ H + 2 lambda OMEGA + tau I ] (-1) dir -# -# coef0 -> wrk5 ; dir -> wrk6 -# - - while(iter tol) { - - - iter = iter + 1 - call scopy(nk,coef,1,wrk5,1) - - # Find Newton direction for new estimate - - # Put the cumulative hazard into wrk4 - do i=1,nsx { wrk4(i) = 0. } - do i=1,n { wrk4(isx(i)) = wrk4(isx(i)) + cumhaz(i)*m(i)*w(i) } - - # Put the regression estimate into smreg - icoef = 1 - do i=1,nsx { - smreg(i) = exp(bvalu(knot,coef,nk,4,0,sx(i),icoef,wrk1))} - - call shessc(bvcd,sd,ld,nd,nbvcd, - nk,smreg,wrk4,sm,nsx,hs0,hs1,hs2,hs3) - call sdircx(sc,wrk4,sm,nsx,nk, - bvcd,sd,ld,nd,nbvcd, - smreg, - lambda,sg0,sg1,sg2,sg3, - wrk5,wrk6) - - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - - call spbfa(abd,ld4,nk,3,ierr) - - if(ierr.ne.0) { return } - do j=1,nk { coef(j) = wrk6(j) } - call spbsl(abd,ld4,nk,3,coef) - do j=1,nk { coef(j) = wrk5(j) - coef(j) } - - - # Find a new estimate in the direction of coef - # On return smreg contains the estimated proportionality factor - - call srchcx(x,t,c,m,sm,w,n, - sx,sc,nsx,isx, - cmin,cmax,coef,wrk5,wrk6,nk, - abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - smreg,knot, - lambda,pl,plo,ierr,wrk0) - - if(ierr==-1) {go to 10} - - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - cchk = abs(pl-plo)/pl - plo = pl - - - } - -10 continue - - - if(ierr > 0 ) {crit=10. ; return} - - - # Compute the criterion function if requested - - if(icrit==0) { return} - - else { # Compute Linearized CV score - - do i=1,nsx { wrk4(i) = 0. } - do i=1,n { wrk4(isx(i)) = wrk4(isx(i)) + cumhaz(i)*w(i)*m(i) } - do i=1,nsx { - smreg(i) = exp(bvalu(knot,coef,nk,4,0,sx(i),icoef,wrk1))} - - call shessc(bvcd,sd,ld,nd,nbvcd, - nk,smreg,wrk4,sm,nsx,hs0,hs1,hs2,hs3) - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) - - if(ierr > 0) {return} - - # Get Leverages - - call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) - - do i=1,nsx { xv = sx(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,wrk0) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - - cvlik(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 } - - do i=1,nsx { lev(i) = cvlik(i)*wrk4(i)*smreg(i)/sm } - - # One-step Newton-Raphson Updates - - do i=1,nsx {cvlik(i) = log(smreg(i)) + - (cvlik(i)*(smreg(i)*wrk4(i)-sc(i))/sm)/(1.-lev(i)) } - - - # Kulbach-Liebner Cross-Validation Score - - crit = 0. - do i=1,nsx{crit = crit + exp(cvlik(i))*wrk4(i)-sc(i)*cvlik(i) } - crit = crit/sm - - } - - - -return -end - -subroutine gxcv(x,t,c,m,w,sm,n,sx,sc,isx,nsx,coef,knot,nk,smreg,cumhaz, - bvcd,sd,ld,nd,nbvcd, - gchi,lev,se,lambda,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,hess,ld4,ldnk,wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw) - -real x(n),t(n),c(n),m(n),w(n),sm, - sx(nsx),sc(nsx), - coef(nk),knot(nk+4),smreg(n),cumhaz(n), - bvcd(nbvcd), - gchi,lev(n),se(n),lambda, - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - abd(ld4,nk),hess(ldnk,nk), - wrk0(nw),wrk1(nw),wrk2(nw),wrk3(nw),wrk4(nw), - wrk5(nw),wrk6(nw) - -integer n,isx(n),nsx,nk,nbvcd, - sd(nk),ld(nk),nd(nk),ld4,ldnk,nw - - - - # Local Variables - -real bvalu,sdot,sqrt,b0,b1,b2,b3,vnikx(4,1),work(15), - eps,xv -integer i,icoef,ierr,id,idum,ileft,mflag,j,k,ilo - - - # Setup Preconditioner - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) ; id=1 - - - # Compute pointwise standard errors - - do j=1,nk { wrk6(j) = 0. } ; ilo = 0. ; eps = .1e-9 - - do i=1,nsx { - xv = sx(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - j=ileft-3 - wrk6(j ) = b0 ; wrk6(j+1) = b1 - wrk6(j+2) = b2 ; wrk6(j+3) = b3 - - call pcg(hess,abd,wrk6,wrk5,ldnk,ld4,nk,.000001, - wrk0,wrk1,wrk2,wrk3,wrk4,id) - - se(i) = sdot(nk,wrk5,1,wrk6,1)/sm - write(17,*) sx(i),se(i) - - j=ileft-3 - wrk6(j ) = 0. ; wrk6(j+1) = 0. - wrk6(j+2) = 0. ; wrk6(j+3) = 0. } - - - # Global Chi-squared - - call sdircx(sc,wrk4,sm,nsx,nk,bvcd,sd,ld,nd,nbvcd, - smreg,lambda,sg0,sg1,sg2,sg3,coef,wrk6) - - call pcg(hess,abd,wrk6,wrk5,ldnk,ld4,nk,.000001, - wrk0,wrk1,wrk2,wrk3,wrk4,id) - - write(17,'("Chi-squared")') - gchi = sdot(nk,wrk5,1,wrk6,1)*sm - write(17,*) gchi - - - - # Individual Leverages - - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk, - 4,0,x(i),icoef,wrk1))*w(i)} - - do i=1,n { lev(i) = m(i)*cumhaz(i)*smreg(i)*se(isx(i))**2. } - write(17,'("Leverages")') - write(17,*) (x(i),lev(i),i=1,n) - - -return -end - -subroutine gxcv(x,t,c,m,w,sm,n,sx,sc,isx,nsx,coef,knot,nk,smreg,cumhaz, - bvcd,sd,ld,nd,nbvcd, - gchi,lev,se,lambda,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,hess,ld4,ldnk,wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw) - -real x(n),t(n),c(n),m(n),w(n),sm, - sx(nsx),sc(nsx), - coef(nk),knot(nk+4),smreg(n),cumhaz(n), - bvcd(nbvcd), - gchi,lev(n),se(n),lambda, - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - abd(ld4,nk),hess(ldnk,nk), - wrk0(nw),wrk1(nw),wrk2(nw),wrk3(nw),wrk4(nw), - wrk5(nw),wrk6(nw) - -integer n,isx(n),nsx,nk,nbvcd, - sd(nk),ld(nk),nd(nk),ld4,ldnk,nw - - - - # Local Variables - -real bvalu,sdot,sqrt,b0,b1,b2,b3,vnikx(4,1),work(15), - eps,xv -integer i,icoef,ierr,id,idum,ileft,mflag,j,k,ilo - -write(17,'("n nk sm ")') -write(17,*) n,nk,sm -write(17,'(" ")') - # Setup Preconditioner - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) ; id=1 - - - # Compute pointwise standard errors - - do j=1,nk { wrk6(j) = 0. } ; ilo = 0. ; eps = .1e-9 - - do i=1,nsx { - xv = sx(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - j=ileft-3 - wrk6(j ) = b0 ; wrk6(j+1) = b1 - wrk6(j+2) = b2 ; wrk6(j+3) = b3 - - call pcg(hess,abd,wrk6,wrk5,ldnk,ld4,nk,.00000001, - wrk0,wrk1,wrk2,wrk3,wrk4,id) - - se(i) = sdot(nk,wrk5,1,wrk6,1)/sm - write(17,*) sx(i),se(i) - - j=ileft-3 - wrk6(j ) = 0. ; wrk6(j+1) = 0. - wrk6(j+2) = 0. ; wrk6(j+3) = 0. } - -if(n==n) stop - - # Global Chi-squared - - call sdircx(sc,wrk4,sm,nsx,nk,bvcd,sd,ld,nd,nbvcd, - smreg,lambda,sg0,sg1,sg2,sg3,coef,wrk6) - - call pcg(hess,abd,wrk6,wrk5,ldnk,ld4,nk,.000001, - wrk0,wrk1,wrk2,wrk3,wrk4,id) - - write(17,'("Chi-squared")') - gchi = sdot(nk,wrk5,1,wrk6,1)*sm - write(17,*) gchi - - - - # Individual Leverages - - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk, - 4,0,x(i),icoef,wrk1))*w(i)} - - do i=1,n { lev(i) = m(i)*cumhaz(i)*smreg(i)*se(isx(i))**2. } - write(17,'("Leverages")') - write(17,*) (x(i),lev(i),i=1,n) - - -return -end //GO.SYSIN DD bart/elabsrcx.r echo bart/evalle.r 1>&2 sed >bart/evalle.r <<'//GO.SYSIN DD bart/evalle.r' 's/^-//' -subroutine evalle(t,c,m,sm,w,n,le) - - -integer n - -real t(n),c(n),m(n),w(n),sm, - le - - # Local Variables - -real rit -integer i,k - - - # Computes the -ve log Partial Likelihood - - # Variables are assumed to be sorted in increasing order relative to t - # violation is not checked and is not fatal - - - i = 1 ; le = 0. - rit = 0. - for(k=i;k<=n;k=k+1) - {rit = rit + m(k)*w(k) } - - if(c(i)>0.) { le = le + c(i)*log(rit) - c(i)*log(w(i)) } - - for(i=2;i<=n;i=i+1){ - - if(c(i)>0.){ - if(t(i)>t(i-1)) - { rit = 0. - for(k=i;k<=n;k=k+1) - {rit = rit + m(k)*w(k) }} - - le = le + c(i)*log(rit) - c(i)*log(w(i)) } } - - le = le/sm - -return -end //GO.SYSIN DD bart/evalle.r echo bart/evalpl.r 1>&2 sed >bart/evalpl.r <<'//GO.SYSIN DD bart/evalpl.r' 's/^-//' -subroutine evalpl(x,t,c,m,sm,w,n, - cmin,cmax,coef,coef0,nk, - smreg,knot, - lambda,sg0,sg1,sg2,sg3, - wrk0,pl) - - -real x(n),t(n),c(n),m(n),sm,w(n), - cmin,cmax,coef(nk),coef0(nk), - smreg(n),knot(nk+4), - lambda,sg0(nk),sg1(nk),sg2(nk),sg3(nk), - wrk0(nk),pl - - -integer n,nk - - - - - # Local Variables - - -real work(15),bvalu,corr, - exp,le,sdot,cmin,cmax -integer i,icoef,j - - - -# Evaluate the estimate in coef and compute the penalized partial likelihood -# -# t(n) is assumed to be nondecreasing -# - - do j=1,nk { - if(coef(j)cmax+2.9) { coef(j) = cmax + 2.9 }} - - # Adjust Coefficients aren't violated - - corr = 0. ; icoef = 1 - do i=1,n{ corr = corr + bvalu(knot,coef,nk,4,0,x(i),icoef,work) } - corr = corr/n - - do j=1,nk { coef(j) = coef(j) - corr - if(coef(j)cmax+2.9) { coef(j) = cmax + 2.9 }} - - corr = 0. - do i=1,n{ corr = corr + bvalu(knot,coef,nk,4,0,x(i),icoef,work) } - corr = corr/n - do j=1,nk { coef(j) = coef(j) - corr } ; icoef = 1 - - - - # Proportionality Function - - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk,4,0,x(i),icoef,work))*w(i)} - - - # Value of pl - - call evalle(t,c,m,sm,smreg,n,le) - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk0) - - pl = le + lambda*sdot(nk,coef,1,wrk0,1) - - -return -end //GO.SYSIN DD bart/evalpl.r echo bart/evlkdn.r 1>&2 sed >bart/evlkdn.r <<'//GO.SYSIN DD bart/evlkdn.r' 's/^-//' - -subroutine evlkdn(w,n,coef,nk,smden,knot,grid,gwts,ngpt, - bvcd,sd,nd,ld,nbvcd,le) - -integer n,nk,ngpt, - sd(nk),nd(nk),ld(nk),nbvcd, - - i,icoef,j - - -real w(n),coef(nk),smden(ngpt), - knot(nk+4),grid(ngpt),gwts(ngpt),bvcd(nbvcd),le, - - gq1,work(15),bvalu,xv,exp,log,lcon,con - - - icoef = 1 ; do i=1,ngpt { xv = grid(i) - smden(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,work(1))) } - - con = gq1(smden,gwts,ngpt) ; lcon = log(con) - - do j=1,nk { coef(j) = coef(j) - lcon } - do i=1,ngpt { smden(i) = smden(i)/con } - - - - # Value of le - - le = 1. - do j=1,nk { if(nd(j)>0) { - do i=1,nd(j){ le = le-coef(j)*w(ld(j)+i-1)*bvcd(sd(j)+i-1)}}} - - - - -return -end //GO.SYSIN DD bart/evlkdn.r echo bart/evlkhz.r 1>&2 sed >bart/evlkhz.r <<'//GO.SYSIN DD bart/evlkhz.r' 's/^-//' -subroutine evlkhz(rn,c,n,coef,nk,smhaz,knot,grid,qwts,ngpt, - bvcd,sd,nd,ld,nbvcd,le) - -integer n,nk,ngpt, - sd(nk),nd(nk),ld(nk),nbvcd, - - i,icoef,j - - -real rn(ngpt),c(n),coef(nk),smhaz(ngpt), - knot(nk+4),grid(ngpt),qwts(ngpt),bvcd(nbvcd),le, - - gq2,work(15),bvalu,xv,exp - - - icoef = 1 ; do i=1,ngpt { xv = grid(i) - smhaz(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,work(1))) } - - le = gq2(smhaz,rn,qwts,ngpt) - - do j=1,nk { if(nd(j)>0) { - do i=1,nd(j){ le=le-coef(j)*c(ld(j)+i-1)*bvcd(sd(j)+i-1)}}} - - - -return -end //GO.SYSIN DD bart/evlkhz.r echo bart/gden.r 1>&2 sed >bart/gden.r <<'//GO.SYSIN DD bart/gden.r' 's/^-//' -subroutine gden(den,cum,grid,ngpt,x,n,cseed,tseed) - -real den(ngpt),cum(ngpt),grid(ngpt),x(n) -integer n,ngpt,cseed,tseed - - - # Local - -double precision rans -real y,v,uni -integer i,k,indx - - - cum(1) = 0. - do k=2,ngpt { cum(k) = cum(k-1) + - (den(k)+den(k-1))*(grid(k)-grid(k-1))/2. } - - - do i=1,n { - - # CDF method - uni = rans() ; v = uni - - # Find y s.t. cum(y) = v - - if(v <= cum(1)) {y = grid(1) } - else if(v>=cum(ngpt)) {y = grid(ngpt)} - - else if(cum(1)&2 sed >bart/ghaz.r <<'//GO.SYSIN DD bart/ghaz.r' 's/^-//' -subroutine ghaz(haz,lcen,cum,grid,ngpt,x,c,n,cseed,tseed) - -real haz(ngpt),lcen,cum(ngpt),grid(ngpt),x(n),c(n) -integer n,ngpt,cseed,tseed - - - # Local - -double precision rans -real y,z,v,uni -integer i,k,indx - - cum(1) = 0. - do k=2,ngpt { cum(k) = cum(k-1) + - (haz(k)+haz(k-1))*(grid(k)-grid(k-1))/2. } - - - do i=1,n { - # Censoring time - uni = rans() ; z = -lcen*log(uni) - - - - - # Failure time - uni = rans() ; v = -log(uni) - - # Find y s.t. cum(y) = v - - if(v <= cum(1)) {y = grid(1) } - else if(v>=cum(ngpt)) {y = grid(ngpt)} - - else if(cum(1)=grid(ngpt)){c(i) = 0. } - - } - -return -end //GO.SYSIN DD bart/ghaz.r echo bart/grhss.r 1>&2 sed >bart/grhss.r <<'//GO.SYSIN DD bart/grhss.r' 's/^-//' -subroutine grhss(t,c,m,sm,w,cumhaz,n,x,ldx,grad,hess,ldp,p) - -real t(n),c(n),m(n),sm,w(n),cumhaz(n), - x(ldx,p), - grad(p),hess(ldp,p) - -integer n,p,ldp,ldx - - - - - # Computes the Gradient and Hessian for the Finite Dimensional Model - # Variables are assumed to be sorted in increasing order relative to t - - - - - # Local Variables - -real ritj,ritk,rit,addjk -integer i,j,k,l - - -# Gradient and Hessian Computation just sum over distinct times - - - - - # Gradient and easy part of the Hessian - - do j=1,p { grad(j) = 0. - - for(i=1;i<=n;i=i+1){ - grad(j) = grad(j) + x(i,j)*(m(i)*cumhaz(i)*w(i)-c(i)) } - - grad(j) = grad(j)/sm - - do k=j,p { hess(j,k) = 0. - - for(i=1;i<=n;i=i+1){ - hess(j,k) = hess(j,k) + x(i,j)*x(i,k)*m(i)*cumhaz(i)*w(i) } - - hess(j,k) = hess(j,k)/sm }} - - - - # Finite Sample Contribution to the Hessian - # imposes a constraint on the parameters - - - do j=1,p { - do k=j,p { addjk = 0. - - i = 1 - rit = 0. ; ritj = 0. ; ritk = 0. - for(l=i;l<=n;l=l+1) - { rit = rit + m(l)*w(l) - ritj = ritj + w(l)*m(l)*x(l,j) - ritk = ritk + w(l)*m(l)*x(l,k) } - - if(c(i)>0) { addjk = addjk + c(i)*(ritj/rit)*(ritk/rit) } - - - for(i=2;i<=n;i=i+1){ - - if(c(i)>0){ - if(t(i)>t(i-1)) - { rit = 0. ; ritj = 0. ; ritk = 0. - for(l=i;l<=n;l=l+1) - { rit = rit + w(l)*m(l) - ritj = ritj + w(l)*m(l)*x(l,j) - ritk = ritk + w(l)*m(l)*x(l,k) }} - - addjk = addjk + c(i)*(ritj/rit)*(ritk/rit) } } - - hess(j,k) = hess(j,k) - addjk/sm }} - -return -end //GO.SYSIN DD bart/grhss.r echo bart/gxv.r 1>&2 sed >bart/gxv.r <<'//GO.SYSIN DD bart/gxv.r' 's/^-//' -subroutine gxv(x,t,c,m,w,sm,n,coef,knot,nk,smreg,cumhaz, - crit,lev,se,lambda,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,wrk1,wrk2,wrk3,wrk4,nw) - -real x(n),t(n),c(n),m(n),w(n),sm, - coef(nk),knot(nk+4),smreg(n),cumhaz(n), - crit,lev(n),se(n),lambda, - 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), - wrk0(nw),wrk1(nw),wrk2(nw),wrk3(nw),wrk4(nw) - -integer n,nk,ld4,ldnk,nw - - - # Local Variables - -integer i,icoef,ileft,mflag,j,ilo -real bvalu,sqrt,b0,b1,b2,b3,vnikx(4,1),eps,xv,sqrt,yv,wv,df - - - - - - # Cumulative Hazard - - icoef = 1 - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk, - 4,0,x(i),icoef,wrk0))*w(i)} - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - do i=1,n {smreg(i) = smreg(i)/w(i) } - - - # Individual Leverages - - call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) - - ilo=1 ; eps =.1e-9 - do i=1,n { xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,wrk0) - b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) - - se(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 } - - do i=1,n { lev(i) = se(i)*cumhaz(i)*smreg(i)*w(i)*m(i)/sm } - df = 0. ; do i=1,n { df = df + lev(i) } - if(df<2.) { do i=1,n { lev(i) = (2./df)*lev(i) } ; df = 2.} - do i=1,n { se(i) = sqrt(se(i)/sm) } - - - # Goodness of Fit Score - - crit = 0. - - do i=1,n{ - if(cumhaz(i)>.1e-9){ - wv = sqrt( cumhaz(i)*w(i)*m(i)*smreg(i)+.00001 ) - yv = cumhaz(i)*w(i)*m(i)*smreg(i) - crit = crit + ((c(i)-yv)/wv)**2. - - } } - - crit = crit/sm - - -return -end //GO.SYSIN DD bart/gxv.r echo bart/icreg.r 1>&2 sed >bart/icreg.r <<'//GO.SYSIN DD bart/icreg.r' 's/^-//' -subroutine icreg(x,t,w,c,n,sx,nx,tmed,coef,knot,nk, - smreg,tw,cw,ww,ecum) - - -integer n,nx,nk - -real x(n),t(n),w(n),c(n),sx(nx),tmed,coef(nk),knot(nk+4),smreg(n), - tw(n),cw(n),ww(n),ecum(n) # work vectors - - # Local - -real icoef,bvalu,height,bquad,xv,delt,abs -integer i,indk,nv,iv, - j,js,k,klow,khigh,nxt,nh,max0 - - -# Histogram Cox Regression Estimate - - # For each x-bin a the cumulative hazard is estimated - # the log of the conditional cumulative hazard at x is - # proportional to the Cox regression - - - - - # nh = # of bins = (1/sd)*nx**(1/5) - - nh = max0((1./(.74*(sx((nx/4)*3)-sx(nx/4))))*nx**(.2),1) - - - klow = 1 - # Find indk s.t. knot(indk+3)>=sx(klow) - do j=1,nk-2{ - if(knot(j+3)>=sx(klow)) {indk=j;break} } - nxt = indk - 1 ; js = indk - - do k=2,nh { khigh = 1 + (k-1)*(nx-1)/(nh-1) - - # Pick out data in this x-bin - # and compute Cumulative Hazard - - nv = 0 - - do i=1,n { if(sx(klow)<=x(i) & x(i) <= sx(khigh) ) - { nv = nv + 1 ; tw(nv) = t(i) - cw(nv) = c(i) ; ww(nv) = w(i) }} - - call scopy(nv,tw,1,ecum,1) - call ssort(tw,cw,nv,2) ; call ssort(ecum,ww,nv,2) - - # Find the nearest value to tmed - iv = 1 ; delt = abs(tw(1)-tmed) - for(i=2;i<=nv;i=i+1) {if(abs(tw(i)-tmed)&2 sed >bart/icreg1.r <<'//GO.SYSIN DD bart/icreg1.r' 's/^-//' -subroutine icreg1(x,t,w,c,m,n,coef,knot,nk, - smreg,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4,mx) - -integer n,nk -integer i,icoef,j,ierr - - - -real x(n),t(n),w(n),c(n),m(n),coef(nk),knot(nk+4),smreg(n),cumhaz(n), - wrk0(n+2),wrk1(n),wrk2(n+6),wrk3(n),wrk4(n+2),mx(n)#work vectors - -real beta,hess(1,1),svdhp(1,1),d,grad,beta0,bvalu,bquad,xv - -# Linear Cox Regression Estimate - - typ = 0. ; do i=1,n { typ = typ + x(i) } ; typ = typ/n - do i=1,n { mx(i) = x(i) } - call slvcxp(mx,typ,n,n,1,t,c,m,w,n,beta,smreg,cumhaz, - wrk0,wrk1,wrk2,wrk3,wrk4,beta0,grad,hess,svdhp,d,1,ierr) - - do j=1,nk { coef(j) = beta*knot(j+2) } - - - - # Adjust to have mean zero - - bquad = 0. ; icoef = 1 - do i=1,n { bquad = bquad + bvalu(knot,coef,nk,4,0,x(i),icoef,wrk1) } - bquad = bquad/n - - do j=1,nk { coef(j) = coef(j) - bquad } - icoef = 1 ; do i=1,n { xv = x(i) - smreg(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,wrk1)) } - - -return -end //GO.SYSIN DD bart/icreg1.r echo bart/iden.r 1>&2 sed >bart/iden.r <<'//GO.SYSIN DD bart/iden.r' 's/^-//' -subroutine iden(x,w,range,n,knot,nk,grid,qwts,ngpt,coef,smden,ecum) - - -integer n,nk,ngpt, - - i,icoef,indk,j,js,k,klow,khigh,nxt,nh,max0 - -real x(n),w(n),range,knot(nk+4),grid(ngpt),qwts(ngpt), - coef(nk),smden(ngpt),ecum(n), - - lcon,con,height,gq1,bvalu,xv,work(15) - - - # nh = # of bins = (1/sd)*n**(1/5) - - nh = max0((1./(.74*(x((n/4)*3)-x(n/4))))*n**(.20),1) - - # Cumulative - - ecum(1) = w(1) - do i=2,n { ecum(i) = ecum(i-1) + w(i) } - con = ecum(n) - do i=1,n { ecum(i) = ecum(i)/con } - - # Histogram Density estimate - - klow = 1 - # Find indk s.t. knot(indk+3)>=x(klow) - do j=1,nk-2{ - if(knot(j+3)>=x(klow)) {indk=j;break} } - nxt = indk - 1 ; js = indk - - - do k=2,nh { khigh = 1 + (k-1)*(n-1)/(nh-1) - - height = log( (ecum(khigh)-ecum(klow))/(x(khigh)-x(klow)) ) - - for(j=indk;j<=nk-2;j=j+1) - - { if(x(klow)<= knot(j+3) & knot(j+3) < x(khigh) ) - - { coef(j+1) = height ; nxt = nxt + 1 } - - else { break }} - - indk = nxt + 1 ; klow = khigh } - - # Fix up the begining and ends - - for(j=1;j<=js-1;j=j+1) { coef(j+1) = coef(js+1) } - for(j=nxt;j<=nk-2;j=j+1) { coef(j+1) = coef(nxt) } - - coef(1) = coef(2) - log(2.) - coef(nk) = coef(nk-1) - log(2.) - - - - # Evaluate and Normalize the density estimate - - icoef = 1 ; do i=1,ngpt { xv = grid(i) - smden(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,work(1))) } - con = gq1(smden,qwts,ngpt) ; lcon = log(con) - - do i=1,ngpt { smden(i) = smden(i)/con } - do i=1,nk { coef(i) = coef(i) - lcon } - - -return -end //GO.SYSIN DD bart/iden.r echo bart/iden1.r 1>&2 sed >bart/iden1.r <<'//GO.SYSIN DD bart/iden1.r' 's/^-//' -subroutine iden1(x,w,n,knot,nk,grid,qwts,ngpt,coef,smden) - - -integer n,nk,ngpt, - - i,icoef,j - -real x(n),w(n),knot(nk+4),grid(ngpt),qwts(ngpt), - coef(nk),smden(ngpt),mw,mxw,sdxw,sqrt, - - lcon,con,gq1,bvalu,xv,work(15) - - - # Gaussian Density Estimate - - mw = 0. ; do i=1,n { mw = mw + w(i) } - mxw = 0. ; do i=1,n { mxw = mxw + x(i)*w(i) } ; mxw = mxw/mw - sdxw = 0. ; do i=1,n { sdxw = sdxw + w(i)*(x(i)-mxw)**2 } - sdxw = sqrt(sdxw/mw) - - - do j=1,nk { coef(j) = 0.-.5*((knot(j+2)-mxw)/sdxw)**2. } - - - - # Evaluate and Normalize the density estimate - - icoef = 1 ; do i=1,ngpt { xv = grid(i) - smden(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,work(1))) } - con = gq1(smden,qwts,ngpt) ; lcon = log(con) - - - do i=1,ngpt { smden(i) = smden(i)/con } - do i=1,nk { coef(i) = coef(i) - lcon } - - -return -end //GO.SYSIN DD bart/iden1.r echo bart/ihaz.r 1>&2 sed >bart/ihaz.r <<'//GO.SYSIN DD bart/ihaz.r' 's/^-//' -subroutine ihaz(x,w,c,range,n,knot,nk,grid,qwts,ngpt, - rn,coef,smhaz,norm,ecum,yw1) - -integer n,nk,ngpt, - - i,indk,icoef, - ix,j,js,k,klow,khigh,nxt,nh,max0 - -real x(n),w(n),c(n),range,knot(nk+4),grid(ngpt),qwts(ngpt), - rn(ngpt),coef(nk),smhaz(ngpt), - norm,ecum(n),yw1(n), - - gq1,bvalu,con,lcon,height,work(15),xv - - - - - - # nh = # of bins = (1/sd)*n**(1/5) - - nh = max0((1./(.74*(x((n/4)*3)-x(n/4))))*n**(.20),1) - - - # Evaluate rn on the grid - - ix = 1 - rn(1) = 0. ; do i=1,n { rn(1) = rn(1) + w(i) } - do i=2,ngpt { - - if(grid(i) > x(ix)) - { rn(i) = amax1(rn(i-1)-w(ix),0) - - for(k=ix+1;k<=n;k=k+1) - { if(grid(i)>x(k)) { rn(i)=amax1(rn(i)-w(k),0) } - else { ix = k ; break } } - - } - - else { rn(i) = rn(i-1) } - - } - - - - - # Empirical Cumulative Hazard Breslow-Johansen NPMLE - - yw1(1) = w(1) - do j=2,n{ yw1(1) = yw1(1) + w(j) } - do j=2,n{ yw1(j) = yw1(j-1) - w(j-1) } - - ecum(1) = c(1)/yw1(1) - do i=2,n { ecum(i) = ecum(i-1) + c(i)/yw1(i) } - - - - # Histogram Type Log-Hazard estimate - - klow = 1 - # Find indk s.t. knot(indk+3)>=x(klow) - do j=1,nk-2{ - if(knot(j+3)>=x(klow)) {indk=j;break} } - nxt = indk - 1 ; js = indk - - do k=2,nh { khigh = 1 + (k-1)*(n-1)/(nh-1) - - height = log( .1e-3+(ecum(khigh)-ecum(klow))/(x(khigh)-x(klow)) ) - - - for(j=indk;j<=nk-2;j=j+1) - { - if(x(klow)<= knot(j+3) & knot(j+3) < x(khigh) ) - { coef(j+1) = height ; nxt = nxt + 1 } - - else {break} - - } - - indk = nxt + 1 ; klow = khigh - - } - - - # Fix up the begining and ends - - for(j=1;j<=js-1;j=j+1) { coef(j+1) = coef(js+1) } - for(j=nxt;j<=nk-2;j=j+1) { coef(j+1) = coef(nxt) } - - coef(1) = coef(2) - log(2.) - coef(nk) = coef(nk-1) - log(2.) - - - - - - # Evaluate the hazard estimate - - icoef = 1 ; do i=1,ngpt { xv = grid(i) - smhaz(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,work(1)))} - - - # Normalization for Hazards ecum(n)-ecum(1) ~ scum(n)-scum(1) - - norm = 1./(ecum(n)-ecum(1)) - con = gq1(smhaz,qwts,ngpt)*norm - lcon= log(con) - - do i=1,ngpt { smhaz(i) = smhaz(i)/con } - do i=1,nk { coef(i) = coef(i) - lcon } - - -return -end //GO.SYSIN DD bart/ihaz.r echo bart/ihaz1.r 1>&2 sed >bart/ihaz1.r <<'//GO.SYSIN DD bart/ihaz1.r' 's/^-//' -subroutine ihaz1(x,range,w,c,n,knot,nk,grid,qwts,ngpt, - coef,smhaz,yw2,yw1,wrk0,wrk1,wrk2,wrk3,wrk4) - -integer n,nk,ngpt, - - i,j,icoef - -real x(n),range,w(n),c(n),knot(nk+4), - coef(nk),smhaz(ngpt),grid(ngpt),qwts(ngpt), - yw2(n),yw1(n), - wrk0(n),wrk1(n),wrk2(n+6),wrk3(n),wrk4(n+2), - - bvalu,work(15),xv, - gam0,gam1,alpha,beta,gq1,con,lcon,norm, - sywx,sxwx,swx,swy,sw,amax1 - - - - # Weighted Linear Regression on Cumulative Hazard - - - do i=1,n { yw1(i) = 1. } - call breslw(x,c,yw1,w,n,yw2,wrk0,wrk1,wrk2,wrk3,wrk4) - - sw = 0. - do i=1,n { sw =sw + w(i) } - swx = 0. - do i=1,n { swx=swx + w(i)*x(i) } - swy = 0. - do i=1,n { swy=swy + w(i)*log(yw2(i)+.001) } - sxwx = 0. - do i=1,n { sxwx=sxwx+w(i)*(x(i)-swx/sw)**2 } - sywx = 0. - do i=1,n { sywx=sywx+w(i)*(x(i)-swx/sw)*log(yw2(i)+.001) } - gam1 = sywx/sxwx ; gam0 = swy/sw - gam1*(swx/sw) - - beta = amax1(gam1,.0001) - alpha = gam0 + log(beta) - - - do j=1,nk { coef(j) = alpha + beta*knot(j+2) } - - - - - # Evaluate the hazard estimate - - icoef = 1 ; do i=1,ngpt { xv = grid(i) - smhaz(i) = exp(bvalu(knot,coef,nk,4,0,xv,icoef,work(1)))} - - - # Normalization for Hazards ecum(n)-ecum(1) ~ scum(n)-scum(1) - - norm = 1./(yw2(n)-yw2(1)) - con = gq1(smhaz,qwts,ngpt)*norm - lcon= log(con) ; alpha = alpha - lcon - - - do i=1,ngpt { smhaz(i) = smhaz(i)/con } - do i=1,nk { coef(i) = coef(i) - lcon } - - - -return -end //GO.SYSIN DD bart/ihaz1.r echo bart/lamv.r 1>&2 sed >bart/lamv.r <<'//GO.SYSIN DD bart/lamv.r' 's/^-//' -subroutine lamv(lambda,spar,ratio,dfup,xu,dflo,xl,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - - -integer n,nk,ld4,ierr, - cnt,cnt1 - -real lambda,spar,x(n),w(n),knot(nk+4),lev(n), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - abd(ld4,nk),p1ip(ld4,nk), - ratio,xu,dfup,xl,dflo, - tdes,dof,abs, - x0,x1,x2,x3,y0,y1,y2,y3,y0b,amax1,amin1 - - - - # - # Find the value of lambda for which the degrees of freedom - # is dflo - (1-spar)*(dfup) - # where dfup and dflo are the largest and smallest admissible - # values for the degrees of freedom - - - # Target Value for the degrees of freedom - - tdes = spar*dfup + (1.-spar)*dflo - - - - # Initialize x0 and x1 - - y1 = 1. ; x1 = .1 - while( y1 > 0. ) { x1 = amin1(4.*(x1+.5),xu) ; cnt = cnt + 1 - call dfr(dof,x1,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - y1 = dof/tdes - 1. - if(abs(y1) <.1e-3) { lambda = ratio*x1 ; return}} - - - x0 = x1 ; y0 = y1; cnt = 0 - while( y0 < 0. ) { x0 = amax1(x0/4.,xl) ; cnt = cnt + 1 - call dfr(dof,x0,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - y0 = dof/tdes - 1. - if(abs(y0) <.1e-3) { lambda = ratio*x0 ; return}} - - # x0 is smaller than the correct value for lambda ( y0 > 0. ) - # x1 is larger than the correct value for lambda ( y1 < 0. ) - - - # The Method of False Position Starts Here - - cnt = 0 - while ( abs(y1) > .01 ) { cnt = cnt + 1 - - x2 = (y1/(y1-y0))*x0 + (y0/(y0-y1))*x1 - - call dfr(dof,x2,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - y2 = dof/tdes - 1. - if(abs(y2) <.1e-3) { lambda = ratio*x2 ; return} - - if(cnt>7) { - if(abs(y2)<.1) { y2 = .0} - else if(abs(y1)<.1) { y2 = .0 ; x2 = x1 } - else if(abs(y0)<.1) { y2 = .0 ; x2 = x0 } - else { #write(18,*) x0,y0,x2,y2,x1,y1 - #call flush(18) - y2 = .0 - if(abs(y0) 0 & cnt1 <= 7) - { y0b = .5*y0b ; cnt1 = cnt1+1 - x3 = (y0b/(y0b-y3))*x3 + (y3/(y3-y0b))*x0 - call dfr(dof,x3,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - y3 = dof/tdes - 1. - if(abs(y3) <.1e-3) { lambda = ratio*x3 ; return}} - - if(cnt1>7) { - if(abs(y3)<.1) { y3 = .00001} - else { #write(18,*) x0,y0b,x2,y2,x3,y3 - #call flush(18) - y3 = .0 } } - - x1 = x2 ; y1 = y2 - x0 = x3 ; y0 = y3 - } - - } # End of False Position Loop - - - # Best value corresponds to x1 - - call dfr(dof,x1,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - - lambda = x1*ratio - - -return -end - - -subroutine dfr(dof,xval,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - -integer n,nk,ld4,ierr, - i,ileft,mflag,j,ilo - -real dof,xval,ratio,x(n),w(n), - knot(nk+4),lev(n), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - abd(ld4,nk),p1ip(ld4,nk), - wrk0(15),rcond,t1,t2, - xb,w0,w1, - b0,b1,b2,b3,vnikx(4,1),eps,xv,lam - - - - lam = ratio*xval - - do i=1,nk { abd(4,i) = hs0(i)+lam*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lam*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lam*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lam*sg3(i) } - - - call spbco(abd(1,1),4,nk,3,rcond,p1ip(1,1),ierr) - - if(ierr!=0 .or. rcond<=.1e-5) - - { t1 = 0. ; t2 = 0. - do i=2,nk-1 { t1 = t1 + hs0(i) } - do i=2,nk-1 { t2 = t2 + sg0(i) } - if(t1>t2*lam) { # lam small - interpolate - ierr=-1 ; dof = 1.*min0(n,nk) - do i=1,n { lev(i) = dof/n }} - else { # lam large - regress - ierr=1 ; dof = 2. - xb = 0. ; w0 = 0. ; w1 = 0. - do i=1,n { w0 = w0 + w(i)**2. ; xb = xb + x(i)*w(i)**2. } - xb = xb/w0 - do i=1,n { w1 = w1 + (x(i)-xb)**2. * w(i)**2. } - do i=1,n { lev(i) = w(i)**2/w0 + (x(i)-xb)**2. * w(i)**2./w1}} - - - } - - - if(ierr!=0) {return} - - - - # Individual Leverages - - call sinerp(abd,ld4,nk,p1ip,vnikx,4,0) - - ilo=1 ; eps =.1e-9 - do i=1,n { xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,wrk0) - 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 } - - do i=1,n { lev(i) = lev(i)*w(i)**2 } - - dof = 0. ; do i=1,n { dof = dof + lev(i) } - if(dof<2.) {do i=1,n { lev(i) = (2./dof)*lev(i) } ; dof = 2. } - -return -end - - - -subroutine uplo(xu,dfup,xl,dflo,ratio,x,w,lev,n,knot,nk, - hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,ld4,ierr) - - -integer n,nk,ld4,ierr,i,it - -real xu,dfup,xl,dflo,ratio, - x(n),w(n),lev(n),knot(nk+4), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - abd(ld4,nk),p1ip(ld4,nk), - rcond,lam,t1,t2, - amin1,amax1,x0,x1,x2,y0,y1,y2,a,b - - - # - # Find the smallest and largest allowable degrees of freedom - # - - - t1 = 0. ; t2 = 0. - do i=2,nk-1 {t1 = t1 + hs0(i) ; t2 = t2 + sg0(i) } - ratio = t1/t2 - - # Upper bound - xu = 1. ; it = 0 - while(it<=7) { lam = xu*ratio ; it = it + 1 - do i=1,nk { abd(4,i) = hs0(i)+lam*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lam*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lam*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lam*sg3(i) } - call spbco(abd,4,nk,3,rcond,p1ip(1,1),ierr) - if(ierr!=0 .or. rcond<=.1e-5) { break } - else { x0 = xu ; y0 = rcond/(.1e-5) - 1. - xu = xu*10. } } - - # Fine tune xu a bit - x1 = xu ; y1 = rcond/(.1e-5) - 1. ; it = 0 - a = x0 ; b = x1 - while(it<=9) { x2 = (y1/(y1-y0))*x0 + (y0/(y0-y1))*x1 - x2 = amax1(x2,a) ; x2 = amin1(x2,b) - it = it+1 ; lam = x2*ratio - do i=1,nk { abd(4,i) = hs0(i)+lam*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lam*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lam*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lam*sg3(i) } - call spbco(abd,4,nk,3,rcond,p1ip(1,1),ierr) - y2 = rcond/(.1e-5) - 1. - x0 = x1 ; y0 = y1 ; x1 = x2 ; y1 = y2 } - - - xu = x2 - - call dfr(dfup,xu,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - - - # Lower bound - xl = 1. ; it = 0 - while(it<=6) { lam = xl*ratio ; it = it + 1 - do i=1,nk { abd(4,i) = hs0(i)+lam*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lam*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lam*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lam*sg3(i) } - call spbco(abd,4,nk,3,rcond,p1ip(1,1),ierr) - if(ierr!=0 .or. rcond<=.1e-5) { break } - else { xl = xl*.1 } } - - # Fine tune xl a bit - xl = sqrt(10.)*xl ; lam = xl*ratio ; rcond=0. - do i=1,nk { abd(4,i) = hs0(i)+lam*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lam*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lam*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lam*sg3(i) } - call spbco(abd,4,nk,3,rcond,p1ip(1,1),ierr) - if(ierr!=0 .or. rcond<=.1e-5) { xl = xl*sqrt(10.) } - - call dfr(dflo,xl,ratio,x,w,n,knot,nk, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,ld4,ierr) - - -return -end //GO.SYSIN DD bart/lamv.r echo bart/lamval.r 1>&2 sed >bart/lamval.r <<'//GO.SYSIN DD bart/lamval.r' 's/^-//' -subroutine lamval(lambda,spar,ratio,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - -real lambda,spar,ratio,x(n),t(n),c(n),m(n),w(n),sm, - knot(nk+4),smreg(n),cumhaz(n),lev(n), - 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),wrk0(nw) - -integer n,nk,ld4,ldnk,nw,ierr - - - # Local Variables - -real tdes,dof,abs,x0,x1,x2,x3,y0,y1,y2,y3,y0b -integer cnt,cnt1 - - - - - # - # Find the value of lambda for which the degrees of freedom - # is spar*nk - # - - - - # Target Value for the degrees of freedom - - tdes = 2.01 + spar*(nk-4.01) - - # Initialize x0 and x1 - - y1 = 1. ; x1 = .1 - while( y1 > 0. ) { x1 = 4.*(x1+.5) ; cnt = cnt + 1 - call df(dof,x1,ratio,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - if(ierr!=0) { return } - y1 = dof/tdes - 1.} - - - x0 = x1 ; y0 = y1; cnt = 0 - while( y0 < 0. ) { x0 = x0/4. ; cnt = cnt + 1 - call df(dof,x0,ratio,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - if(ierr!=0) { return } - y0 = dof/tdes - 1.} - - # x0 is smaller than the correct value for lambda ( y0 > 0. ) - # x1 is larger than the correct value for lambda ( y1 < 0. ) - - - - # The Method of False Position Starts Here - - cnt = 0 - while ( abs(y1) > .01 ) { cnt = cnt + 1 - - x2 = (y1/(y1-y0))*x0 + (y0/(y0-y1))*x1 - - call df(dof,x2,ratio,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - if(ierr!=0) { return } - y2 = dof/tdes - 1. - - if(cnt>7) { - if(abs(y2)<.1) { y2 = .0} - else if(abs(y1)<.1) { y2 = .0 ; x2 = x1 } - else if(abs(y0)<.1) { y2 = .0 ; x2 = x0 } - else { write(18,*) x0,y0,x2,y2,x1,y1 - call flush(18) - y2 = .0 - if(abs(y0) 0) - { y0b = .5*y0b ; cnt1 = cnt1+1 - x3 = (y0b/(y0b-y3))*x3 + (y3/(y3-y0b))*x0 - call df(dof,x3,ratio,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - if(ierr!=0) { return } - y3 = dof/tdes - 1. } - - if(cnt1>7) { - if(abs(y3)<.1) { y3 = .00001} - else { write(18,*) x0,y0b,x2,y2,x3,y3 - call flush(18) - y3 = .0 } } - - } - - x1 = x2 ; y1 = y2 - x0 = x3 ; y0 = y3 - } # End of False Position Loop - -lambda = x1*ratio - -return -end - - -subroutine df(dof,xval,ratio,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - -real dof,xval,ratio,x(n),t(n),c(n),m(n),w(n),sm, - knot(nk+4),smreg(n),cumhaz(n),lev(n), - 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), - wrk0(nw) - -integer n,nk,ld4,ldnk,nw,ierr - - - # Local Variables - -real b0,b1,b2,b3,vnikx(4,1),eps,xv,lam -integer i,ileft,mflag,j,ilo - - - - lam = ratio*xval - - do i=1,nk { abd(4,i) = hs0(i)+2.*lam*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lam*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lam*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lam*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) - if(ierr!=0) { - write(4,'("error in Cholesky decomposition in df routine ")') - write(4,*) ierr - return } - - - # Individual Leverages - - call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) - - ilo=1 ; eps =.1e-9 - do i=1,n { xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,wrk0) - 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 } - - do i=1,n { lev(i) = lev(i)*cumhaz(i)*smreg(i)*w(i)*m(i)/sm } - - dof = 0. ; do i=1,n { dof = dof + lev(i) } - if(dof<2.) { do i=1,n { lev(i) = (2./dof)*lev(i) } ; dof = 2. } - -return -end //GO.SYSIN DD bart/lamval.r echo bart/lsamp1.r 1>&2 sed >bart/lsamp1.r <<'//GO.SYSIN DD bart/lsamp1.r' 's/^-//' -subroutine lsamp1(f,n,sbd,ld,ns,lev,wx,we,wf,knot,coef) - -real f(n),sbd(ld,n),lev(n), - we(n),wf(n),wx(n),knot(ns+6),coef(ns+2) - -integer n,ld,ns - - - # Local - -real bvalu -integer i,icoef,ifix,ix,j - - - # Set up sampled Leverages - - call scopy(n,f,1,wf,1) - do i=1,n { wx(i) = i } - call ssort(wf,wx,n,2) - - do i=1,4 { knot(i) = wf(1) } - for(i=2;i<=ns-1;i=i+1){ - knot(i+3) = wf(1+(i-1)*((n-1)/(ns-1)))} - do i=1,4 { knot(ns+2+i) = wf(n) } - - - do i=1,ns{ ix = ifix(wx(1+(i-1)*((n-1)/(ns-1)))) - do j=1,n { we(j) = 0. } ; we(ix) = wf(ix) - call spbsl(sbd,ld,n,2,we) - coef(i+1) = min(.999,we(ix)) } - coef(1) = coef(2) ; coef(ns+2) = coef(ns+1) - - # Compute Leverages at all values - - icoef = 1 - do i=1,n { lev(i) = bvalu(knot,coef,ns+2,4,0,f(i),icoef,wx(1)) } - - - - -return -end //GO.SYSIN DD bart/lsamp1.r echo bart/pcg.r 1>&2 sed >bart/pcg.r <<'//GO.SYSIN DD bart/pcg.r' 's/^-//' -subroutine pcg(a,m,b,x,lda,ldm,n,tol,r0,r1,z0,z1,p,ierr) - -real a(lda,n),m(ldm,n),b(n),x(n),tol, - r0(n),r1(n),z0(n),z1(n),p(n) - -integer lda,ldm,n,ierr - - - # Local - -real sdot,snrm2,alpha,beta -integer j,k - - - - - # - # A Pre-conditioned Conjugate Gradient Algorithm for solving - # - # A x = b - # - # A and M are symetric. - # M is a 3-banded pre-conditioner the cholesky decomposition of - # M is computed. - # - # The algorithm is described in Golub and Van Loan page 374. - # (Ap is stored in z0) - - - call spbfa(m,ldm,n,3,ierr) - if(ierr!=0) { - write(4,'("error in Cholesky decomposition in pcg routine ")') - return } - - - # Initialization - - do j=1,n { x(j) = 0. } - call scopy(n,b,1,r1,1) - - - - - do k=1,n { - if(snrm2(n,r1,1)1) { beta = sdot(n,z1,1,r1,1)/sdot(n,z0,1,r0,1) - do j=1,n { p(j) = z1(j) + beta*p(j) }} - - else { beta = 0. ; call scopy(n,z1,1,p,1) } - - - #Use z0 to store ap - do j=1,n { z0(j) = sdot(n,a(1,j),1,p,1) } - alpha = sdot(n,z1,1,r1,1)/sdot(n,p,1,z0,1) - - do j=1,n { x(j) = x(j) + alpha*p(j) } - - call scopy(n,r1,1,r0,1) - call saxpy(n,-alpha,z0,1,r1,1) - call scopy(n,z1,1,z0,1) - } - - - } - - -return -end //GO.SYSIN DD bart/pcg.r echo bart/sbart.r 1>&2 sed >bart/sbart.r <<'//GO.SYSIN DD bart/sbart.r' 's/^-//' -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) - - - - - # A Cubic B-spline Smoothing routine. - -# -# The algorithm minimises: -# -# (1/n) * sum ws(i)**2 * (ys(i)-sz(i))**2 + lambda* int ( sz"(xs) )**2 dxs -# -# lambda is a function of the spar which is assumed to be between -# 0 and 1 - - - # Input - -# n number of data points -# ys(n) vector of length n containing the observations -# ws(n) vector containing the weights given to each data point -# xs(n) vector containing the ordinates of the observations - - -# nk number of b-spline coefficients to be estimated -# nk <= n+2 -# knot(nk+4) vector of knot points defining the cubic b-spline basis. -# To obtain full cubic smoothing splines one might -# have (provided the xs-values are strictly increasing) - - -# spar penalised likelihood smoothing parameter -# ispar indicator saying if spar is supplied or to be estimated -# lspar, uspar lower and upper values for spar 0.,1. are good values -# tol used in Golden Search routine - -# isetup setup indicator - -# icrit indicator saying which cross validation score -# is to be computed - -# ld4 the leading dimension of abd (ie ld4=4) -# ldnk the leading dimension of p2ip (not referenced) - - - # Output - -# coef(nk) vector of spline coefficients -# sz(n) vector of smoothed z-values -# lev(n) vector of leverages -# crit either ordinary of generalized CV score -# ier error indicator -# ier = 0 ___ everything fine -# ier = 1 ___ spar too small or too big -# problem in cholesky decomposition - - - - # Working arrays/matrix -# xwy X'Wy -# hs0,hs1,hs2,hs3 the diagonals of the X'WX matrix -# sg0,sg1,sg2,sg3 the diagonals of the Gram matrix -# abd(ld4,nk) [ X'WX+lambda*SIGMA] in diagonal form -# p1ip(ld4,nk) inner products between columns of L inverse -# p2ip(ldnk,nk) all inner products between columns of L inverse -# L'L = [X'WX+lambdaSIGMA] NOT REFERENCED - - -integer n,nk,isetup,icrit,ispar,ld4,ldnk,ier, - i - - -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), - - - lamr,dof,t1,t2, - a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, - fu,fv,fw,fx,x, - ax,bx - - - - - - - - - - # Compute SIGMA, X' W**2 X, X' W**2 z, and lamr - - # SIGMA -> sg0,sg1,sg2,sg3 - # X' W**2 X -> hs0,hs1,hs2,hs3 - # X' W**2 Z -> xwy - - if(isetup==0){ - 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 i=2,nk-1 {t1 = t1 + hs0(i) ; t2 = t2 + sg0(i) } - lamr = t1/t2 - call dfr(dof,1.,lamr,xs,ws,ngpt,knot,nk,lev, - hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,ld4,ier) - - isetup = 1 } - - - - - # Compute estimate - - - if(ispar==1) { # Value of spar supplied - - call sslvrg(xs,ys,ws,n,knot,nk, - coef,sz,lev,crit,icrit, - spar,lamr,dof, - xwy, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,ier) - - return } - - else { - - # Use Forsythe Malcom and Moler routine to minimise criterion - - ax=lspar ; bx=uspar # f denotes the value of the criterion - -# -# an approximation x to the point where f attains a minimum on -# the interval (ax,bx) is determined. -# -# -# input.. -# -# ax left endpoint of initial interval -# bx right endpoint of initial interval -# f function subprogram which evaluates f(x) for any x -# in the interval (ax,bx) -# tol desired length of the interval of uncertainty of the final -# result ( .ge. 0.0) -# -# -# output.. -# -# fmin abcissa approximating the point where f attains a minimum -# -# -# the method used is a combination of golden section search and -# successive parabolic interpolation. convergence is never much slower -# than that for a fibonacci search. if f has a continuous second -# derivative which is positive at the minimum (which is not at ax or -# bx), then convergence is superlinear, and usually of the order of -# about 1.324.... -# the function f is never evaluated at two points closer together -# than eps*abs(fmin) + (tol/3), where eps is approximately the square -# root of the relative machine precision. if f is a unimodal -# function and the computed values of f are always unimodal when -# separated by at least eps*abs(x) + (tol/3), then fmin approximates -# the abcissa of the global minimum of f on the interval ax,bx with -# an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, -# then fmin may approximate a local, but perhaps non-global, minimum to -# the same accuracy. -# this function subprogram is a slightly modified version of the -# algol 60 procedure localmin given in richard brent, algorithms for -# minimization without derivatives, prentice - hall, inc. (1973). -# -# -# real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w -# real fu,fv,fw,fx,x -# -# c is the squared inverse of the golden ratio -# - c = 0.5*(3. - sqrt(5.0)) -# -# eps is approximately the square root of the relative machine -# precision. -# - eps = 1.00 - 10 eps = eps/2.00 - tol1 = 1.0 + eps - if (tol1 .gt. 1.00) go to 10 - eps = sqrt(eps) -# -# initialization -# - 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,lamr,dof, - xwy, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,ier) - - fx = crit - fv = fx - fw = fx -# -# main loop starts here -# - 20 xm = 0.5*(a + b) - tol1 = eps*abs(x) + tol/3.0 - tol2 = 2.0*tol1 -# -# check stopping criterion -# - if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 -# -# is golden-section necessary -# - if (abs(e) .le. tol1) go to 40 -# -# fit parabola -# - r = (x - w)*(fx - fv) - q = (x - v)*(fx - fw) - p = (x - v)*q - (x - w)*r - q = 2.00*(q - r) - if (q .gt. 0.0) p = -p - q = abs(q) - r = e - e = d -# -# is parabola acceptable -# - 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 - if (p .le. q*(a - x)) go to 40 - if (p .ge. q*(b - x)) go to 40 -# -# a parabolic interpolation step -# - d = p/q - u = x + d -# -# f must not be evaluated too close to ax or bx -# - if ((u - a) .lt. tol2) d = sign(tol1, xm - x) - if ((b - u) .lt. tol2) d = sign(tol1, xm - x) - go to 50 -# -# a golden-section step -# - 40 if (x .ge. xm) e = a - x - if (x .lt. xm) e = b - x - d = c*e -# -# f must not be evaluated too close to x -# - 50 if (abs(d) .ge. tol1) u = x + d - if (abs(d) .lt. tol1) u = x + sign(tol1, d) - - spar = u - call sslvrg(xs,ys,ws,n,knot,nk, - coef,sz,lev,crit,icrit, - spar,lamr,dof, - xwy, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,ier) - - fu = crit -# -# update a, b, v, w, and x -# - if (fu .gt. fx) go to 60 - if (u .ge. x) a = x - if (u .lt. x) b = x - v = w - fv = fw - w = x - fw = fx - x = u - fx = fu - go to 20 - 60 if (u .lt. x) a = u - if (u .ge. x) b = u - if (fu .le. fw) go to 70 - if (w .eq. x) go to 70 - if (fu .le. fv) go to 80 - if (v .eq. x) go to 80 - if (v .eq. w) go to 80 - go to 20 - 70 v = w - fv = fw - w = u - fw = fu - go to 20 - 80 v = u - fv = fu - go to 20 -# -# end of main loop -# - 90 continue ; spar = x ; crit = fx - return } - - - - - -return -end //GO.SYSIN DD bart/sbart.r echo bart/sbcox.r 1>&2 sed >bart/sbcox.r <<'//GO.SYSIN DD bart/sbcox.r' 's/^-//' -subroutine sbcox(xs,ts,cs,ms,ws,n, - sx,sc,isx,nsx, - bvcd,nbvcd,sd,nd,ld, - coef0,coef1,hcoef,coef,knot,nk, - smreg,cumhaz,cvlik,lev, - crit,icrit, - spar,ispar,lspar,uspar,tol, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw,ier) - - - - - # A Cubic B-spline Cox-Regression Smoothing routine. - -# -# The program minimises: -# -# Cox Partial Likelihood for sz(xs) + lambda* int ( sz"(xs) )**2 dxs -# -# lambda is a function of the spar which is assumed to be between -# 0 and 1 - - - # Input - -# n number of data points -# xs,ts,cs,ms,ws x-values, times , censoring indicators -# multiplicities and prior weights for the data -# sorted by time -# sx,sc,isx sorted x-values and corresponding cs's isx(n) -# gives the sx order for xs,ts,cs,ms,ws. -# hcoef(nk) coefficients of initial estimate - - -# nk number of b-spline coefficients to be estimated -# nk <= n+2 -# knot(nk+4) vector of knot points defining the cubic b-spline basis. -# -# bvcd(ngpt) b-splines evaluated at the sorted data points sx -# sd,nd,ld auxillary information used in processing bvcd - - -# spar penalised likelihood smoothing parameter -# ispar indicator saying if spar is supplied or to be estimated -# lspar, uspar lower and upper values for spar 0.,1. are good values -# tol used in Golden Search routine - -# icrit indicator saying whether information or mse criterion -# is requested - -# ld4 the leading dimension of abd (ie ld4=4) -# ldnk the leading dimension of p2ip (not referenced) - - - # Output - -# coef(nk) spline coefficients -# smreg(n) smoothed hazard estimates at xs -# cumhaz(n) Breslow Cumulative Basline Hazard Estimates at xs -# on input cumhaz has an initial estimate -# cvlik(nsx) cross validated Cox regression estimates at sx -# lev(nsx) Leverage Values at sx -# -# crit Information or MSE score -# ier error indicator -# ier = -1 ___ Newton method failed to decrease pl -# ier = 0 ___ everything fine -# ier = 1 ___ spar too small or too big -# problem in cholesky decomposition - - - - # Working arrays/matrix - -# sg0,sg1,sg2,sg3 the diagonals of the Gram matrix (dim=nk) -# hs0,hs1,hs2,hs3 the diagonals of the Hessian (dim=nk) -# abd(ld4,nk) [ X'WX+lambda*SIGMA] in diagonal form -# p1ip(ld4,nk) inner products between columns of L inverse -# p2ip(ldnk,nk) all inner products between columns of L inverse -# L'L = [X'WX+lambdaSIGMA] NOT REFERENCED -# wrk0,...wrk6 other work arrays dimension - - -integer isx(n),nsx,n,nk,nbvcd, - sd(nk),ld(nk),nd(nk), - icrit,ispar,ld4,ldnk,nw,ier - - -real xs(n),ts(n),cs(n),ms(n),ws(n), - sx(nsx),sc(nsx), - bvcd(nbvcd), - coef0(nk),coef1(nk),hcoef(nk),coef(nk), - knot(nk+4), - smreg(n),cumhaz(n),cvlik(n),lev(n), - crit,spar,lspar,uspar,tol, - 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), - wrk0(nw),wrk1(nw),wrk2(nw),wrk3(nw),wrk4(nw), - wrk5(nw),wrk6(nw) - - - - # Local variables - -real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, - fu,fv,fw,fx,x, - ax,bx - - - - - - # Compute SIGMA - - call sgram(sg0,sg1,sg2,sg3,knot,nk) - - - - # Compute estimate - - - if(ispar==1) { # Value of spar supplied - - call sslvcx(xs,ts,cs,ms,ws,n, - sx,sc,isx,nsx, - bvcd,nbvcd,sd,nd,ld, - coef0,coef1,hcoef,coef,knot,nk,icrit,crit,spar, - smreg,cumhaz,cvlik,lev, - wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,ier) - - return } - - else { - - # Use Forsythe Malcom and Moler routine to minimise criterion - - ax=lspar ; bx=uspar # f denotes the value of the criterion - -# -# an approximation x to the point where f attains a minimum on -# the interval (ax,bx) is determined. -# -# -# input.. -# -# ax left endpoint of initial interval -# bx right endpoint of initial interval -# f function subprogram which evaluates f(x) for any x -# in the interval (ax,bx) -# tol desired length of the interval of uncertainty of the final -# result ( .ge. 0.0) -# -# -# output.. -# -# fmin abcissa approximating the point where f attains a minimum -# -# -# the method used is a combination of golden section search and -# successive parabolic interpolation. convergence is never much slower -# than that for a fibonacci search. if f has a continuous second -# derivative which is positive at the minimum (which is not at ax or -# bx), then convergence is superlinear, and usually of the order of -# about 1.324.... -# the function f is never evaluated at two points closer together -# than eps*abs(fmin) + (tol/3), where eps is approximately the square -# root of the relative machine precision. if f is a unimodal -# function and the computed values of f are always unimodal when -# separated by at least eps*abs(x) + (tol/3), then fmin approximates -# the abcissa of the global minimum of f on the interval ax,bx with -# an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, -# then fmin may approximate a local, but perhaps non-global, minimum to -# the same accuracy. -# this function subprogram is a slightly modified version of the -# algol 60 procedure localmin given in richard brent, algorithms for -# minimization without derivatives, prentice - hall, inc. (1973). -# -# -# real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w -# real fu,fv,fw,fx,x -# -# c is the squared inverse of the golden ratio -# - c = 0.5*(3. - sqrt(5.0)) -# -# eps is approximately the square root of the relative machine -# precision. -# - eps = 1.00 - 10 eps = eps/2.00 - tol1 = 1.0 + eps - if (tol1 .gt. 1.00) go to 10 - eps = sqrt(eps) -# -# initialization -# - a = ax - b = bx - v = a + c*(b - a) - w = v - x = v - e = 0.0 - - spar = x - call sslvcx(xs,ts,cs,ms,ws,n, - sx,sc,isx,nsx, - bvcd,nbvcd,sd,nd,ld, - coef0,coef1,hcoef,coef,knot,nk,icrit,crit,spar, - smreg,cumhaz,cvlik,lev, - wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,ier) - - fx = crit - fv = fx - fw = fx -# -# main loop starts here -# - 20 xm = 0.5*(a + b) - tol1 = eps*abs(x) + tol/3.0 - tol2 = 2.0*tol1 -# -# check stopping criterion -# - if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 -# -# is golden-section necessary -# - if (abs(e) .le. tol1) go to 40 -# -# fit parabola -# - r = (x - w)*(fx - fv) - q = (x - v)*(fx - fw) - p = (x - v)*q - (x - w)*r - q = 2.00*(q - r) - if (q .gt. 0.0) p = -p - q = abs(q) - r = e - e = d -# -# is parabola acceptable -# - 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 - if (p .le. q*(a - x)) go to 40 - if (p .ge. q*(b - x)) go to 40 -# -# a parabolic interpolation step -# - d = p/q - u = x + d -# -# f must not be evaluated too close to ax or bx -# - if ((u - a) .lt. tol2) d = sign(tol1, xm - x) - if ((b - u) .lt. tol2) d = sign(tol1, xm - x) - go to 50 -# -# a golden-section step -# - 40 if (x .ge. xm) e = a - x - if (x .lt. xm) e = b - x - d = c*e -# -# f must not be evaluated too close to x -# - 50 if (abs(d) .ge. tol1) u = x + d - if (abs(d) .lt. tol1) u = x + sign(tol1, d) - - spar = u - call sslvcx(xs,ts,cs,ms,ws,n, - sx,sc,isx,nsx, - bvcd,nbvcd,sd,nd,ld, - coef0,coef1,hcoef,coef,knot,nk,icrit,crit,spar, - smreg,cumhaz,cvlik,lev, - wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,ier) - - fu = crit -# -# update a, b, v, w, and x -# - if (fu .gt. fx) go to 60 - if (u .ge. x) a = x - if (u .lt. x) b = x - v = w - fv = fw - w = x - fw = fx - x = u - fx = fu - go to 20 - 60 if (u .lt. x) a = u - if (u .ge. x) b = u - if (fu .le. fw) go to 70 - if (w .eq. x) go to 70 - if (fu .le. fv) go to 80 - if (v .eq. x) go to 80 - if (v .eq. w) go to 80 - go to 20 - 70 v = w - fv = fw - w = u - fw = fu - go to 20 - 80 v = u - fv = fu - go to 20 -# -# end of main loop -# - 90 continue ; spar = x ; crit = fx - return } - - - - - -return -end //GO.SYSIN DD bart/sbcox.r echo bart/sbden.r 1>&2 sed >bart/sbden.r <<'//GO.SYSIN DD bart/sbden.r' 's/^-//' -subroutine sbden(xs,ws,n,knot,nk, - coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev,wts, - crit,dcrit,icrit, - spar,ispar,lspar,uspar,tol, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - - - - # A Cubic B-spline Log-Density Smoothing routine. - -# -# The program minimises: -# -# int exp{sz(xs)} dxs - sum ws(i)*sz(xs(i)) / sum ws(i) -# + lambda* int ( sz"(xs) )**2 dxs -# -# lambda is a function of the spar which is assumed to be between -# 0 and 1 - - - # Input - -# n number of data points -# ws(n) vector containing the weights given to each data point -# xs(n) vector containing the ordinates of the observations -# hcoef(nk) coefficients of initial histogram estimate -# coef0,coef1 rough and smooth estimates - - -# nk number of b-spline coefficients to be estimated -# nk <= n+2 -# knot(nk+4) vector of knot points defining the cubic b-spline basis. -# -# grid(ngpt) grid points all computations are done relative to this grid -# qwts(ngpt) Gaussian quadrature weights -# -# bvcg(ngpt) b-splines evaluated at the grid points -# sg,ng,lg auxillary information used in processing bvcg -# -# bvcd(ngpt) b-splines evaluated at the data points xs -# sd,nd,ld auxillary information used in processing bvcd - - -# spar penalised likelihood smoothing parameter -# ispar indicator saying if spar is supplied or to be estimated -# lspar, uspar lower and upper values for spar 0.,1. are good values -# tol used in Golden Search routine - -# icrit indicator saying whether information or mse criterion -# is requested - -# ld4 the leading dimension of abd (ie ld4=4) -# ldnk the leading dimension of p2ip (not referenced) - - - # Output - -# coef(nk) vector of spline coefficients -# smden(ngpt) vector of smoothed density estimates -# lev(ngpt) vector of smoothed leverages -# wts(ngpt) vector of weights used in computing leverages -# cvden(n) vector of cross validated density estimates -# evaluated at the xs -# -# crit Information and MSE Linearized CV-scores -# ier error indicator -# ier = -1 ___ Newton method failed to decrease pl -# ier = 0 ___ everything fine -# ier = 1 ___ spar too small or too big -# problem in cholesky decomposition - - - - # Working arrays/matrix - -# sg0,sg1,sg2,sg3 the diagonals of the Gram matrix (dim=nk) -# hs0,hs1,hs2,hs3 the diagonals of the Hessian (dim=nk) -# abd(ld4,nk) [ X'WX+lambda*SIGMA] in diagonal form -# p1ip(ld4,nk) inner products between columns of L inverse -# p2ip(ldnk,nk) all inner products between columns of L inverse -# L'L = [X'WX+lambdaSIGMA] NOT REFERENCED -# wrk1,wrk2,wrk3 other work arrays dimension nk - - - -integer n,nk,ngpt, - nbvcg,sg(nk),ng(nk),lg(nk), - nbvcd,sd(nk),nd(nk),ld(nk), - dcrit,icrit,ispar,ld4,ldnk,ier - - - -real xs(n),ws(n), - knot(nk+4), - coef0(nk),coef1(nk),hcoef(nk), - grid(ngpt),qwts(ngpt),bvcg(nbvcg),bvcd(nbvcd), - coef(nk),smden(ngpt),cvden(n), - lev(ngpt),wts(ngpt), - crit(dcrit),spar,lspar,uspar,tol, - 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), - wrk1(nk),wrk2(nk),wrk3(nk), - - a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, - fu,fv,fw,fx,x,ax,bx - - - - - - - - # Compute SIGMA - - # SIGMA -> sg0,sg1,sg2,sg3 - - call sgram(sg0,sg1,sg2,sg3,knot,nk) - - - # Compute estimate - - - if(ispar==1) { # Value of spar supplied - - call sslvdn(xs,ws,n,knot,nk, - coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev,wts, - crit,dcrit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - return } - - else { - - # Use Forsythe Malcom and Moler routine to minimise criterion - - ax=lspar ; bx=uspar # f denotes the value of the criterion - -# -# an approximation x to the point where f attains a minimum on -# the interval (ax,bx) is determined. -# -# -# input.. -# -# ax left endpoint of initial interval -# bx right endpoint of initial interval -# f function subprogram which evaluates f(x) for any x -# in the interval (ax,bx) -# tol desired length of the interval of uncertainty of the final -# result ( .ge. 0.0) -# -# -# output.. -# -# fmin abcissa approximating the point where f attains a minimum -# -# -# the method used is a combination of golden section search and -# successive parabolic interpolation. convergence is never much slower -# than that for a fibonacci search. if f has a continuous second -# derivative which is positive at the minimum (which is not at ax or -# bx), then convergence is superlinear, and usually of the order of -# about 1.324.... -# the function f is never evaluated at two points closer together -# than eps*abs(fmin) + (tol/3), where eps is approximately the square -# root of the relative machine precision. if f is a unimodal -# function and the computed values of f are always unimodal when -# separated by at least eps*abs(x) + (tol/3), then fmin approximates -# the abcissa of the global minimum of f on the interval ax,bx with -# an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, -# then fmin may approximate a local, but perhaps non-global, minimum to -# the same accuracy. -# this function subprogram is a slightly modified version of the -# algol 60 procedure localmin given in richard brent, algorithms for -# minimization without derivatives, prentice - hall, inc. (1973). -# -# -# real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w -# real fu,fv,fw,fx,x -# -# c is the squared inverse of the golden ratio -# - c = 0.5*(3. - sqrt(5.0)) -# -# eps is approximately the square root of the relative machine -# precision. -# - eps = 1.00 - 10 eps = eps/2.00 - tol1 = 1.0 + eps - if (tol1 .gt. 1.00) go to 10 - eps = sqrt(eps) -# -# initialization -# - a = ax - b = bx - v = a + c*(b - a) - w = v - x = v - e = 0.0 - - spar = x - call sslvdn(xs,ws,n,knot,nk, - coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev,wts, - crit,dcrit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - fx = crit(icrit) - fv = fx - fw = fx -# -# main loop starts here -# - 20 xm = 0.5*(a + b) - tol1 = eps*abs(x) + tol/3.0 - tol2 = 2.0*tol1 -# -# check stopping criterion -# - if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 -# -# is golden-section necessary -# - if (abs(e) .le. tol1) go to 40 -# -# fit parabola -# - r = (x - w)*(fx - fv) - q = (x - v)*(fx - fw) - p = (x - v)*q - (x - w)*r - q = 2.00*(q - r) - if (q .gt. 0.0) p = -p - q = abs(q) - r = e - e = d -# -# is parabola acceptable -# - 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 - if (p .le. q*(a - x)) go to 40 - if (p .ge. q*(b - x)) go to 40 -# -# a parabolic interpolation step -# - d = p/q - u = x + d -# -# f must not be evaluated too close to ax or bx -# - if ((u - a) .lt. tol2) d = sign(tol1, xm - x) - if ((b - u) .lt. tol2) d = sign(tol1, xm - x) - go to 50 -# -# a golden-section step -# - 40 if (x .ge. xm) e = a - x - if (x .lt. xm) e = b - x - d = c*e -# -# f must not be evaluated too close to x -# - 50 if (abs(d) .ge. tol1) u = x + d - if (abs(d) .lt. tol1) u = x + sign(tol1, d) - - spar = u - call sslvdn(xs,ws,n,knot,nk, - coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev,wts, - crit,dcrit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - fu = crit(icrit) -# -# update a, b, v, w, and x -# - if (fu .gt. fx) go to 60 - if (u .ge. x) a = x - if (u .lt. x) b = x - v = w - fv = fw - w = x - fw = fx - x = u - fx = fu - go to 20 - 60 if (u .lt. x) a = u - if (u .ge. x) b = u - if (fu .le. fw) go to 70 - if (w .eq. x) go to 70 - if (fu .le. fv) go to 80 - if (v .eq. x) go to 80 - if (v .eq. w) go to 80 - go to 20 - 70 v = w - fv = fw - w = u - fw = fu - go to 20 - 80 v = u - fv = fu - go to 20 -# -# end of main loop -# - 90 continue ; spar = x ; crit(icrit) = fx - return } - - - - - -return -end //GO.SYSIN DD bart/sbden.r echo bart/sbhaz.r 1>&2 sed >bart/sbhaz.r <<'//GO.SYSIN DD bart/sbhaz.r' 's/^-//' -subroutine sbhaz(xs,ws,cs,rn,n,knot,nk, - norm,coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smhaz,cvhaz,lev,wts, - crit,icrit, - spar,ispar,lspar,uspar,tol, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - - - - # A Cubic B-spline Log-Hazard Smoothing routine. - -# -# The program minimises: -# -# int rn(s) exp{sz(xs)} dxs - sum cs(i)*sz(xs(i)) -# + lambda* int ( sz"(xs) )**2 dxs -# -# lambda is a function of the spar which is assumed to be between -# 0 and 1 - - - # Input - -# n number of data points -# ws(n) vector containing the weights for the i'th case -# cs(n) vector containing the number of failures at each data point -# xs(n) vector containing the ordinates of the observations -# rn(ngpt) empirical risk at each grid point -# hcoef(nk) coefficients of initial estimate -# coef0,coef1 controls on the size of the size of the estimated hazard -# norm normalization for the estimated hazard - - -# nk number of b-spline coefficients to be estimated -# nk <= n+2 -# knot(nk+4) vector of knot points defining the cubic b-spline basis. -# -# grid(ngpt) grid points all computations are done relative to this grid -# qwts(ngpt) quadrature weights -# -# bvcg(ngpt) b-splines evaluated at the grid points -# sg,ng,lg auxillary information used in processing bvcg -# -# bvcd(ngpt) b-splines evaluated at the data points xs -# sd,nd,ld auxillary information used in processing bvcd - - -# spar penalised likelihood smoothing parameter -# ispar indicator saying if spar is supplied or to be estimated -# lspar, uspar lower and upper values for spar 0.,1. are good values -# tol used in Golden Search routine - -# icrit indicator saying whether information or mse criterion -# is requested - -# ld4 the leading dimension of abd (ie ld4=4) -# ldnk the leading dimension of p2ip (not referenced) - - - # Output - -# coef(nk) vector of spline coefficients -# smhaz(ngpt) vector of smoothed hazard estimates -# lev(ngpt) vector of leverages -# wts(ngpt) vector of weights for computing leverages -# cvhaz(n) vector of cross validated hazard estimates -# evaluated at the xs -# -# crit Discrete and Continuous CV score -# ier error indicator -# ier = -1 ___ Newton method failed to decrease pl -# ier = 0 ___ everything fine -# ier = 1 ___ spar too small or too big -# problem in cholesky decomposition - - - - # Working arrays/matrix - -# sg0,sg1,sg2,sg3 the diagonals of the Gram matrix (dim=nk) -# hs0,hs1,hs2,hs3 the diagonals of the Hessian (dim=nk) -# abd(ld4,nk) [ X'WX+lambda*SIGMA] in diagonal form -# p1ip(ld4,nk) inner products between columns of L inverse -# p2ip(ldnk,nk) all inner products between columns of L inverse -# L'L = [X'WX+lambdaSIGMA] NOT REFERENCED -# wrk1,wrk2,wrk3 other work arrays dimension nk - - - -integer n,nk,ngpt, - nbvcg,sg(nk),ng(nk),lg(nk), - nbvcd,sd(nk),nd(nk),ld(nk), - icrit,ispar,ld4,ldnk,ier - - - -real xs(n),ws(n),cs(n),rn(ngpt), - norm,coef0(nk),coef1(nk),hcoef(nk), - knot(nk+4), - grid(ngpt),qwts(ngpt),bvcg(nbvcg),bvcd(nbvcd), - coef(nk),smhaz(ngpt),cvhaz(n), - lev(ngpt),wts(ngpt), - crit(2),spar,lspar,uspar,tol, - 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), - wrk1(nk),wrk2(nk),wrk3(nk), - - a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, - fu,fv,fw,fx,x, - ax,bx - - - - - - - - # Compute SIGMA and the HESSIAN and initialise ratio - - # SIGMA -> sg0,sg1,sg2,sg3 - - - call sgram(sg0,sg1,sg2,sg3,knot,nk) - - - - # Compute estimate - - - if(ispar==1) { # Value of spar supplied - - call sslvhz(xs,ws,cs,rn,n,knot,nk, - norm,coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smhaz,cvhaz,lev,wts, - crit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - return } - - else { - - # Use Forsythe Malcom and Moler routine to minimise criterion - - ax=lspar ; bx=uspar # f denotes the value of the criterion - -# -# an approximation x to the point where f attains a minimum on -# the interval (ax,bx) is determined. -# -# -# input.. -# -# ax left endpoint of initial interval -# bx right endpoint of initial interval -# f function subprogram which evaluates f(x) for any x -# in the interval (ax,bx) -# tol desired length of the interval of uncertainty of the final -# result ( .ge. 0.0) -# -# -# output.. -# -# fmin abcissa approximating the point where f attains a minimum -# -# -# the method used is a combination of golden section search and -# successive parabolic interpolation. convergence is never much slower -# than that for a fibonacci search. if f has a continuous second -# derivative which is positive at the minimum (which is not at ax or -# bx), then convergence is superlinear, and usually of the order of -# about 1.324.... -# the function f is never evaluated at two points closer together -# than eps*abs(fmin) + (tol/3), where eps is approximately the square -# root of the relative machine precision. if f is a unimodal -# function and the computed values of f are always unimodal when -# separated by at least eps*abs(x) + (tol/3), then fmin approximates -# the abcissa of the global minimum of f on the interval ax,bx with -# an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, -# then fmin may approximate a local, but perhaps non-global, minimum to -# the same accuracy. -# this function subprogram is a slightly modified version of the -# algol 60 procedure localmin given in richard brent, algorithms for -# minimization without derivatives, prentice - hall, inc. (1973). -# -# -# real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w -# real fu,fv,fw,fx,x -# -# c is the squared inverse of the golden ratio -# - c = 0.5*(3. - sqrt(5.0)) -# -# eps is approximately the square root of the relative machine -# precision. -# - eps = 1.00 - 10 eps = eps/2.00 - tol1 = 1.0 + eps - if (tol1 .gt. 1.00) go to 10 - eps = sqrt(eps) -# -# initialization -# - a = ax - b = bx - v = a + c*(b - a) - w = v - x = v - e = 0.0 - - spar = x - call sslvhz(xs,ws,cs,rn,n,knot,nk, - norm,coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smhaz,cvhaz,lev,wts, - crit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - fx = crit(icrit) - fv = fx - fw = fx -# -# main loop starts here -# - 20 xm = 0.5*(a + b) - tol1 = eps*abs(x) + tol/3.0 - tol2 = 2.0*tol1 -# -# check stopping criterion -# - if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 -# -# is golden-section necessary -# - if (abs(e) .le. tol1) go to 40 -# -# fit parabola -# - r = (x - w)*(fx - fv) - q = (x - v)*(fx - fw) - p = (x - v)*q - (x - w)*r - q = 2.00*(q - r) - if (q .gt. 0.0) p = -p - q = abs(q) - r = e - e = d -# -# is parabola acceptable -# - 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 - if (p .le. q*(a - x)) go to 40 - if (p .ge. q*(b - x)) go to 40 -# -# a parabolic interpolation step -# - d = p/q - u = x + d -# -# f must not be evaluated too close to ax or bx -# - if ((u - a) .lt. tol2) d = sign(tol1, xm - x) - if ((b - u) .lt. tol2) d = sign(tol1, xm - x) - go to 50 -# -# a golden-section step -# - 40 if (x .ge. xm) e = a - x - if (x .lt. xm) e = b - x - d = c*e -# -# f must not be evaluated too close to x -# - 50 if (abs(d) .ge. tol1) u = x + d - if (abs(d) .lt. tol1) u = x + sign(tol1, d) - - spar = u - call sslvhz(xs,ws,cs,rn,n,knot,nk, - norm,coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcg,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smhaz,cvhaz,lev,wts, - crit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ier) - - fu = crit(icrit) -# -# update a, b, v, w, and x -# - if (fu .gt. fx) go to 60 - if (u .ge. x) a = x - if (u .lt. x) b = x - v = w - fv = fw - w = x - fw = fx - x = u - fx = fu - go to 20 - 60 if (u .lt. x) a = u - if (u .ge. x) b = u - if (fu .le. fw) go to 70 - if (w .eq. x) go to 70 - if (fu .le. fv) go to 80 - if (v .eq. x) go to 80 - if (v .eq. w) go to 80 - go to 20 - 70 v = w - fv = fw - w = u - fw = fu - go to 20 - 80 v = u - fv = fu - go to 20 -# -# end of main loop -# - 90 continue ; spar = x ; crit(icrit) = fx - return } - - - - - -return -end //GO.SYSIN DD bart/sbhaz.r echo bart/sbvc.r 1>&2 sed >bart/sbvc.r <<'//GO.SYSIN DD bart/sbvc.r' 's/^-//' - -subroutine sbvc (knots,nb,data,m,bvec,s,l,n,work) - -# Purpose : compute the cubic B-splines at the data points -# The data are assumed to be in increasing order -# "Natural" cubic B-splines are used - knots -# of multiplicity 3 at the ends. - -real knots(nb+4), data(m), bvec(4*m+nb), work(nb) - -integer nb,m,s(nb),l(nb),n(nb),avail,istart,i,j - - # nb : number of B-splines in the basis - # s(j) : index of first value of B(j) in bvec - # l(j) : index of first data point where B(j) is computed - # n(j) : number of evaluations of B(j) - # avail: next available position in bvec - # istart: where to start looking for data in the support of B(j) - - - avail = 1 ; istart = 1 - do j=1,nb { work(j) = 0.0 } # used to define the j'th B-spline - - - do j=1,nb { # Find l(j) - l(j) = 0 - do i=istart,m { if(data(i)<=knots(j+4)& - data(i)>=knots(j) ) - {l(j)=i ; goto 5 }} -5 continue - - if(l(j)==0) { n(j) = 0 ; goto 15 } - - # Find n(j) - n(j) = 0 - do i=l(j),m { if(data(i)>knots(j+4)) {goto 10} - else { n(j) = n(j) + 1 } } - -10 continue - - s(j) = avail - work(j) = 1.0 - call fill(bvec(s(j)),data(l(j)),n(j),knots,work,nb) - work(j) = 0.0 - avail = s(j) + n(j) ; istart = l(j) - -15 continue - - if(data(m)&2 sed >bart/scvdn.r <<'//GO.SYSIN DD bart/scvdn.r' 's/^-//' -subroutine scvdn(x,w,n,crit,knot,nk, - grid,qwts,ngpt, - bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev, - abd,p1ip,p2ip,ld4,ldnk) - - -integer n,nk,ngpt,nbvcd, - sd(nk),ld(nk),nd(nk), - ld4,ldnk - -integer i,j,ilo,mflag,ileft,icoef - - -real x(n),w(n),crit(4), - knot(nk+4), - grid(ngpt),qwts(ngpt), - bvcd(nbvcd), - coef(nk),smden(ngpt),cvden(n),lev(ngpt), - abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk) - -real bvalu,vnikx(4,1),eps,work(15),gq1,mdf,le, - xv,b0,b1,b2,b3,sw,sw0 - - - eps = .1e-9 ; ilo = 1 - - - call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) - - - - # Leave-out-one estimate - - do i=1,n { xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - 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 - - lev(i) = lev(i)*w(i) - - } - - icoef = 1 ; do i=1,n { xv = x(i) - cvden(i) = bvalu(knot,coef,nk,4,0,xv,icoef,work(1)) - lev(i) } - - - - # Continuous Leverages - - do i=1,ngpt { xv = grid(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - 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 - - lev(i) = lev(i)*qwts(i)*smden(i) - - } - - - - -# Evaluate CV Scores - - - crit(1) = 0. - do i=1,n { crit(1) = crit(1) - cvden(i)*w(i) } - - mdf = 0. - do i=1,ngpt { mdf = mdf + lev(i) } - if(mdf<2.) { do i=1,ngpt {lev(i) = (2./mdf)*lev(i)} ; mdf = 2.} - - le = 1. - do j=1,nk { if(nd(j)>0) { - do i=1,nd(j){ le = le-coef(j)*w(ld(j)+i-1)*bvcd(sd(j)+i-1)}}} - - crit(3) = le + mdf/n - - # Robustified version of the score - le = 0. ; sw = 0. ; sw0 = 0. - icoef = 1 ; do i=1,n { xv = x(i) ; sw = sw + w(i) - if(.01*n<=i&i<=.99*n) { sw0 = sw0 + w(i) - le = le - w(i)*bvalu(knot,coef,nk,4,0,xv,icoef,work(1)) }} - - crit(3) = 1. + (sw/sw0)*le + mdf/n - - - - le = 0. # For ISE score -do j=1,nk { if(nd(j)>0) { -do i=1,nd(j){ le = le-2.*exp(coef(j)*w(ld(j)+i-1)*bvcd(sd(j)+i-1))*w(i)}}} - - do i=1,ngpt { smden(i) = smden(i)**2 } - do i=1,n { cvden(i) = exp(cvden(i)) } - - crit(2) = gq1(smden,qwts,ngpt) - - crit(4) = crit(2) + le + 2. * mdf/n - - do i=1,n { crit(2) = crit(2) - 2.*cvden(i)*w(i) } - - - do i=1,ngpt { smden(i) = sqrt(smden(i)) } - - -return -end //GO.SYSIN DD bart/scvdn.r echo bart/scvhz.r 1>&2 sed >bart/scvhz.r <<'//GO.SYSIN DD bart/scvhz.r' 's/^-//' -subroutine scvhz(x,w,c,rn,n,crit,icrit,knot,nk, - grid,qwts,ngpt, - bvcd,nbvcd,sd,nd,ld, - bvcgr,nbvcg,sg,ng,lg, - coef,smhaz,cvlik,lev, - abd,p1ip,p2ip,ld4,ldnk, - wrk,wrk1,xi) - -integer n,icrit,nk,ngpt,nbvcd,nbvcg, - sd(nk),ld(nk),nd(nk), - sg(nk),lg(nk),ng(nk), - ld4,ldnk, - - i,inbv,min0, - niu,niuo,nil,nilo,gxi,gxio,k,j,ni, - ilo,ileft,mflag,i1,i2,ifix - - -real x(n),w(n),c(n),rn(ngpt),crit(2), - knot(nk+4), - grid(ngpt),qwts(ngpt), - bvcd(nbvcd),bvcgr(nbvcg), - coef(nk),smhaz(ngpt),cvlik(n),lev(ngpt), - abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk), - wrk(nk),wrk1(nk),xi(nk), - - bvalu,bjxi,work(15), - gq1,gq2,cumhaz,lhazx,le,mdf, - xv,eps,b0,b1,b2,b3,vnikx(4,1),work(15) - - - - -if(icrit==1) { - - # A linearized CV score (leave out one type) - - - call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,1) - - niu = 1 ; nil = 1 ; gxi = 1 ; cumhaz = 0. - - do j=1,nk { wrk1(j) = 0. } - - - - do i=1,n { - niuo = niu ; nilo = nil ; gxio = gxi - - if(i==1) { nil = 1 } # nil - else { for(j=nilo;j<=nk;j=j+1) - {if(grid(lg(j)+ng(j))<=x(i-1)){nil=j+1} - else {break} } } - - for(k=gxio;k<=ngpt;k=k+1) - {if(grid(k)<=x(i)) {gxi=k} # gxi - else {break} } - - for(j=niuo;j<=nk;j=j+1) - {if(grid(lg(j))<=x(i)) {niu=j} # niu - else {break} } - - - # Compute xi , cumhaz and lhazx - - for(j=1 ;j< nil;j=j+1) - - { xi(j) = w(i)*wrk(j) } - - - for(j=nil;j<=niu;j=j+1) - - { ni = min0(gxi-lg(j)+1,ng(j)) - if(ni>1){ - wrk(j)=gq2(smhaz(lg(j)),bvcgr(sg(j)),qwts(lg(j)),ni)} - - - if(gxi<=lg(j)+ng(j)) - { wrk1(j) = 1. ; inbv = 1 - bjxi = bvalu(knot,wrk1,nk,4,0,x(i),inbv,work) - wrk1(j) = 0. - - xi(j) = w(i)*wrk(j) - c(i)*bjxi } - - else - - { xi(j) = w(i)*wrk(j) } - - } - - if(gxio&2 sed >bart/sdircx.r <<'//GO.SYSIN DD bart/sdircx.r' 's/^-//' -subroutine sdircx(c,cumhaz,sm,n,nb, - bvcd,sd,ld,nd,nbvcd, - smreg, - lambda,sg0,sg1,sg2,sg3, - coef0,ans) - -real c(n),cumhaz(n),sm, - bvcd(nbvcd), - smreg(n), - lambda,sg0(nb),sg1(nb),sg2(nb),sg3(nb), - coef0(nb),ans(nb) - -integer n,nb,sd(nb),ld(nb),nd(nb),nbvcd - - - # Local Variables - -real sum -integer i,j - - - call cbnm(sg0,sg1,sg2,sg3,coef0,(nb-3),ans) - - - do j=1,nb{ - # discrete part -> sum - sum = 0. - if(nd(j)==0) { go to 10 } - do i=1,nd(j) { sum = sum + - (smreg(ld(j)+i-1)*cumhaz(ld(j)+i-1)-c(ld(j)+i-1))* - bvcd(sd(j)+i-1) } - -10 continue - - ans(j) = 2.*lambda*ans(j) + sum/sm - - } - -return -end //GO.SYSIN DD bart/sdircx.r echo bart/sdirdn.r 1>&2 sed >bart/sdirdn.r <<'//GO.SYSIN DD bart/sdirdn.r' 's/^-//' -subroutine sdirdn(w,n,bvcg,sg,lg,ng,nbvcg,nb,gwts,ngpt, - bvcd,sd,ld,nd,nbvcd, - f1, - lambda,sg0,sg1,sg2,sg3, - coef0,ans) - - -integer n,sg(nb),lg(nb),ng(nb),nbvcg, - nb,ngpt,sd(nb),ld(nb),nd(nb),nbvcd, - - i,j - -real w(n), - bvcg(nbvcg),gwts(ngpt), - bvcd(nbvcd), - f1(ngpt), - lambda,sg0(nb),sg1(nb),sg2(nb),sg3(nb), - coef0(nb),ans(nb), - - sum,gq2 - - - - - call cbnm(sg0,sg1,sg2,sg3,coef0,(nb-3),ans) - - - - do j=1,nb{ - - # discrete part -> sum - sum = 0. - if(nd(j)==0) { go to 10 } - do i=1,nd(j) { sum = sum + - bvcd(sd(j)-1+i)*w(ld(j)-1+i) } - -10 continue - - ans(j) = 2.*lambda*ans(j) + - gq2(f1(lg(j)),bvcg(sg(j)),gwts(lg(j)),ng(j)) - sum - - } - - -return -end //GO.SYSIN DD bart/sdirdn.r echo bart/sdirhz.r 1>&2 sed >bart/sdirhz.r <<'//GO.SYSIN DD bart/sdirhz.r' 's/^-//' -subroutine sdirhz(c,n,bvcg,sg,lg,ng,nbvcg,nb,qwts,ngpt, - bvcd,sd,ld,nd,nbvcd, - f1,f3, - lambda,sg0,sg1,sg2,sg3, - coef0,ans) - -integer n,sg(nb),lg(nb),ng(nb),nbvcg, - nb,ngpt,sd(nb),ld(nb),nd(nb),nbvcd, - - i,j - -real c(n), - bvcg(nbvcg), qwts(ngpt), - bvcd(nbvcd), - f1(ngpt),f3(ngpt), - lambda,sg0(nb),sg1(nb),sg2(nb),sg3(nb), - coef0(nb),ans(nb), - - sum,gq3 - - - - - - call cbnm(sg0,sg1,sg2,sg3,coef0,(nb-3),ans) - - - do j=1,nb{ - # discrete part -> sum - sum = 0. - if(nd(j)==0) { go to 10 } - do i=1,nd(j) { sum = sum + - bvcd(sd(j)-1+i)*c(ld(j)-1+i) } - -10 continue - - ans(j) = 2.*lambda*ans(j) + - - gq3(f1(lg(j)),f3(lg(j)),bvcg(sg(j)),qwts(lg(j)),ng(j)) - - - sum - - } - -return -end //GO.SYSIN DD bart/sdirhz.r echo bart/search.r 1>&2 sed >bart/search.r <<'//GO.SYSIN DD bart/search.r' 's/^-//' -subroutine search(x,ldx,t,c,m,sm,w,n, - coef,coef0,dir,p, - svdhp,d,wd,ldp, - smreg,le,lo,ierr,wrk0) - - -real x(ldx,p),t(n),c(n),m(n),sm,w(n), - coef(p),coef0(p),dir(p), - svdhp(ldp,p),d(p),wd(p), - smreg(n),le,lo,wrk0(p) - - -integer ldx,n,p,ldp,ierr - - - # Local Variables - -real lm,sdot, - exp,diag,alpham,tinit,tdes,tau,snrm2, - sqrt,nco,abs, - x0,x1,x2,x3, - y0,y1,y2,y3, - y0b -integer i,j,k,cnt,cnt1 - - - - # (THIS IS TOO ELABORATE - HALVING SHOULD BE FINE) - - # - # A trust region implementation of a - # Levenberg-Marquart Line Search for - # a new estimate is used - # - # coef = coef0 - [abd+alpha*I] (-1) dir - # - # where the values of alpha are choosen - # so that the step sizes - # - # norm[coef-coef0]/(norm[coef0]+.01) - # - # are 4**(-k) , k = 1,2, ... 10. - # - # The method of false position (with Illinois Modification) - # is used to find the correct alpha. - # - - - - - - do i=1,n {smreg(i) = exp(sdot(p,coef,1,x(i,1),ldx))*w(i) } - call evalle(t,c,m,sm,smreg,n,le) - - # Size of full step - - nco = snrm2(p,coef0,1)/sqrt(p*1.) - do i=1,p { wd(i) = coef(i) - coef0(i) } - tinit = snrm2(p,wd,1)/(sqrt(p*1.)*(nco+1.)) - - - - # Step Length Search - - - diag = 0. ; do j=1,p { diag = diag+d(j) } - diag = diag/p - - alpham = 0. ; lm = le ; tau = tinit - - x1 = 0. ; y1 = 1. - - - do k = 1, 10 {# Next Step Size - - tdes = tinit*(2.)**(-k) - - # Initialize x0 and x1 - - y1 = 1. - while( y1 > 0. ) { x1 = 4.*(x1+.5) - call stepp(x1,tau,nco,diag,coef,dir,p, - svdhp,d,wd,ldp,wrk0) - y1 = tau/tdes - 1. } - - x0 = x1 ; y0 = y1 - - while( y0 < 0. ) { x0 = x0/4. - call stepp(x0,tau,nco,diag,coef,dir,p, - svdhp,d,wd,ldp,wrk0) - y0 = tau/tdes - 1. } - - # x0 is smaller than the correct value for alpha ( y0 > 0. ) - # x1 is larger than the correct value for alpha ( y1 < 0. ) - - - - # The Method of False Position Starts Here - - cnt = 0 - while ( abs(y1) > .01 ) { cnt = cnt + 1 - - x2 = (y1/(y1-y0))*x0 + (y0/(y0-y1))*x1 - - call stepp(x2,tau,nco,diag,coef,dir,p, - svdhp,d,wd,ldp,wrk0) - - y2 = tau/tdes - 1. - - if(cnt>15) { - if(abs(y2)<.1) { y2 = .0} - else if(abs(y1)<.1) { y2 = .0 ; x2 = x1 } - else if(abs(y0)<.1) { y2 = .0 ; x2 = x0 } - else { write(18,'(" cnt catch ",i3)') k - write(18,*) x0,y0,x2,y2,x1,y1 - call flush(18) - y2 = .0 - if(abs(y0) 0) - { y0b = .5*y0b - x3 = (y0b/(y0b-y3))*x3 + - (y3/(y3-y0b))*x0 - call stepp(x3,tau,nco,diag,coef,dir,p, - svdhp,d,wd,ldp,wrk0) - - y3 = tau/tdes - 1. - - if(cnt1>20) { - if(abs(y3)<.1) { y3 = .00001} - else { write(18,'(" cnt1 catch ",i3)') k - write(18,*) x0,y0b,x2,y2,x3,y3 - call flush(18) - y3 = .0 } } - - } - - x0 = x2 ; y0 = y2 - x1 = x3 ; y1 = y3 } - - - } # End of False Position Loop - - - - do j=1,p { coef(j) = coef0(j) - coef(j) } - - do i=1,n {smreg(i) = exp(sdot(p,coef,1,x(i,1),ldx))*w(i) } - call evalle(t,c,m,sm,smreg,n,le) - - if(le < lm) { alpham = x1 ; lm = le } - - } # Go onto the next step size - - - - - - # - # Choose the Best Step - # - - - le = lm - - if(le > lo ) { ierr = -1 ; le = lo - do j=1,p { coef(j) = coef0(j) } - - # Regression Function - - do i=1,n {smreg(i) = exp(sdot(p,coef,1,x(i,1),ldx))*w(i) } - - } - - - else { - - - call stepp(alpham,tau,nco,diag,coef,dir,p, - svdhp,d,wd,ldp,wrk0) - - if(tau>tinit) { do j=1,p { coef(j) = coef(j)*(tinit/tau) }} - - do j=1,p { coef(j) = coef0(j) - coef(j) } - do i=1,n {smreg(i) = exp(sdot(p,coef,1,x(i,1),ldx))*w(i) } - call evalle(t,c,m,sm,smreg,n,le) - - } - - -return -end //GO.SYSIN DD bart/search.r echo bart/setup.r 1>&2 sed >bart/setup.r <<'//GO.SYSIN DD bart/setup.r' 's/^-//' -subroutine setreg(x,y,w,n,xw,nx,min,range,knot,nk) - -real x(n),y(n),w(n),xw(n), - min,range,knot(n+6) - -integer n,nk,nx - - - # Local - -real eps -integer ycnt,i,k - - - - - - call scopy(n,x,1,xw,1) - call ssort(x,w,n,2) - call ssort(xw,y,n,2) - - range = x(n)-x(1) ; min = x(1) ; eps = .1e-9 - do i=1,n { x(i) = (x(i)-min)/range } - call scopy(n,x,1,xw,1) - - - nx = 1 ; x(nx) = x(1) ; w(nx) = w(1) ; y(nx) = y(1) ; ycnt = 1 - - for(i=2;i<=n;i=i+1) - - { if(xw(i)>x(nx)+eps) - { y(nx) = y(nx)/ycnt - ycnt = w(i) - nx = nx + 1 - x(nx) = x(i) - y(nx) = y(i)*w(i) - w(nx) = w(i) } - else - { y(nx) = y(nx)+y(i)*w(i) - ycnt = ycnt + w(i) - w(nx) = w(nx) + w(i) } - - } - y(nx) = y(nx)/ycnt - - - - - - call sknotl(x,nx,knot,k) ; nk = k-4 - - - - -return -end - - -subroutine setden(x,w,n,a,b,ind,xw,nx, - knot,nk,grid,qwts,ngpt, - bvcd,nbvcd,sd,ld,nd, - bvcg,nbvcg,sg,lg,ng, - ecum,coef0,coef1,smden) - - - -integer ind,n,nk,nx, - ngpt,nbvcd,nbvcg, - sd(nk),ld(nk),nd(nk), - sg(nk),lg(nk),ng(nk), - - i,j,k - -real x(n),w(n),a,b,xw(n+2), - knot(nk+4), - grid(ngpt),qwts(ngpt),bvcd(nbvcd),bvcg(nbvcg), - ecum(n),coef0(nk),coef1(nk),smden(ngpt), - - deltl,deltu,eps,xs,len,sm,av(4),wv(4) - - - - - eps = .1e-9 - call ssort(x,w,n,2) - - -# Collapse to distinct x's - - # Scale x's - deltl = (x(n/10)-x(1))/(.1*n-1.) - deltu = (x(n)-x((n/10)*9))/(.1*n-1.) - if(ind!=1) { a = x(1) - deltl ; b = x(n) + deltu } - do i=1,n { x(i) = (x(i)-a)/(b-a) } - call scopy(n,x,1,xw(2),1) - - - nx = 1 ; x(nx) = x(1) ; w(nx) = w(1) - - for(i=2;i<=n;i=i+1) - { if(xw(i+1)>x(nx)+eps) - { nx=nx+1 ; x(nx)=x(i) ; w(nx)=w(i) } - - else - { w(nx)=w(nx)+w(i)} - } - - -# Sum and scale w(i) - - sm = 0. - do i=1,nx { sm = sm + w(i) } - do i=1,nx { w(i) = w(i)/sm } - - xw(1) = 0. ; call scopy(nx,x,1,xw(2),1) ; xw(nx+2) = 1. - - call sknotl(xw,nx+2,knot,k) ; nk = k-4 - - - - - # Setup the grid (4 Gaussian Quadrature Points between knots ) - # Gaussian quadature points and weights ( standardized interval [-1,1] ) - - av(1) = -.861136 ; av(2) = -.339981 ; av(3) = .339981 ; av(4) = .861136 - wv(1) = .347855 ; wv(2) = .652145 ; wv(3) = .652145 ; wv(4) = .347855 - - ngpt = 0 - - do j=4,nk { xs = knot(j) ; len = (knot(j+1)-knot(j))/2. - do k=1,4 { - ngpt = ngpt + 1 ; grid(ngpt) = xs + len*(1.+av(k)) - qwts(ngpt) = len*wv(k) } - } - - - - # Smooth (ML) estimate - - call iden1(x,w,nx,knot,nk,grid,qwts,ngpt,coef1,smden) - - # Initial histogram estimate (rough estimate) - - call iden(x,w,b-a,nx,knot,nk,grid,qwts,ngpt, - coef0,smden,ecum) - - - - - # Evaluate Basis Vectors on Grid and Data - - nbvcg = 4*ngpt+nk - call sbvc (knot,nk,grid,ngpt,bvcg,sg,lg,ng,xw) - nbvcd = 4*nx+nk - call sbvc (knot,nk,x,nx,bvcd,sd,ld,nd,xw) - - - -return -end - - - - - -subroutine sethaz(x,w,m,c,n,a,b,ind,xw,nx,knot,nk, - grid,qwts,ngpt,rn, - bvcd,nbvcd,sd,ld,nd, - bvcg,nbvcg,sg,lg,ng, - ecum,norm,coef0,coef1,smhaz,yw1, - wrk0,wrk1,wrk2,wrk3,wrk4) - - -integer ind,n,nk,nx, - ngpt,nbvcd,nbvcg, - sd(nk),ld(nk),nd(nk), - sg(nk),lg(nk),ng(nk), - - i,j,k - -real x(n),w(n),m(n),c(n),a,b,xw(n+2), - knot(nk+4), - grid(ngpt),qwts(ngpt),rn(ngpt), - bvcd(nbvcd),bvcg(nbvcg), - ecum(n),coef0(nk),coef1(nk),norm,smhaz(ngpt),yw1(n), - wrk0(n),wrk1(n),wrk2(n+6),wrk3(n),wrk4(n+2), - - deltl,deltu,eps,xs,len,sm,av(4),wv(4) - - - - - eps = .1e-9 - - # Sort data and collapse to distinct x's - - call scopy(n,x,1,xw,1) ; call ssort(xw,w,n,2) - call scopy(n,x,1,xw,1) ; call ssort(xw,m,n,2) - call ssort(x,c,n,2) - - # Scale x's - deltl = (x(n/10)-x(1))/(.1*n-1.) - deltu = (x(n)-x((n/10)*9))/(.1*n-1.) - if(ind!=1) { a = x(1) - deltl ; b = x(n) + deltu } - do i=1,n { x(i) = (x(i)-a)/(b-a) } - call scopy(n,x,1,xw(2),1) - -nx = 1 ; x(nx) = x(1) ; w(nx) = w(1) ; c(nx) = c(1) ; m(nx) = m(1) - - for(i=2;i<=n;i=i+1) - { if(xw(i+1)>x(nx)+eps) - { nx=nx+1 ; x(nx)=x(i) ; m(nx)=m(i) ; w(nx)=w(i) ; c(nx)=c(i) } - else - { m(nx)=m(nx)+m(i) ; w(nx)=w(nx)+w(i) ; c(nx)=c(nx)+c(i) } - } - - - # Sum m(i) and scale w(i) and c(i) - - sm = 0. - do i=1,nx { sm = sm + m(i) } - do i=1,nx { w(i) = w(i)/sm } - do i=1,nx { c(i) = c(i)/sm } - - xw(1) = 0. ; call scopy(nx,x,1,xw(2),1) ; xw(nx+2) = 1. - call sknotl(xw,nx+2,knot,k) ; nk = k-4 - - - # Setup the grid (4 Gaussian Quadrature Points between knots ) - # Gaussian quadature points and weights ( standardized interval [-1,1] ) - - av(1) = -.861136 ; av(2) = -.339981 ; av(3) = .339981 ; av(4) = .861136 - wv(1) = .347855 ; wv(2) = .652145 ; wv(3) = .652145 ; wv(4) = .347855 - - ngpt = 0 - - do j=4,nk { xs = knot(j) ; len = (knot(j+1)-knot(j))/2. - do k=1,4 { - ngpt = ngpt + 1 ; grid(ngpt) = xs + len*(1.+av(k)) - qwts(ngpt) = len*wv(k) } - } - - - # Smooth Weibull Estimate - - call ihaz1(x,b-a,w,c,nx,knot,nk,grid,qwts,ngpt, - coef1,smhaz,ecum,yw1,wrk0,wrk1,wrk2,wrk3,wrk4) - - # Initial histogram estimate - - - call ihaz(x,w,c,b-a,nx,knot,nk,grid,qwts,ngpt, - rn,coef0,smhaz,norm,ecum,yw1) - - - - # Evaluate Basis Vectors on Grid and Data - - nbvcg = 4*ngpt+nk - call sbvc (knot,nk,grid,ngpt,bvcg,sg,lg,ng,xw) - nbvcd = 4*nx+nk - call sbvc (knot,nk,x,nx,bvcd,sd,ld,nd,xw) - - - - -return -end - - - - - -subroutine setcrg(x,t,c,w,m,n,sx,sc,isx,nx,a,b,ind,tmed,coef0,coef1,knot,nk, - cumhaz,smreg,bvcd,nbvcd,sd,ld,nd,xw,wrk0,wrk1,wrk2,wrk3,wrk4) - - -real x(n),t(n),c(n),w(n),m(n), - sx(n),sc(n), - a,b,tmed,coef0(n+2),coef1(n+2),knot(n+6), - cumhaz(n),smreg(n), - bvcd(nbvcd), - xw(n),wrk0(n+2),wrk1(n),wrk2(n+6),wrk3(n),wrk4(n+2) - - -integer ind,isx(n),n,nk,nx, - nbvcd,sd(n),ld(n),nd(n) - - - # Local - -real deltu,deltl,eps -integer ifix,i,k - - - eps = .1e-9 - - - # Get isx - call scopy(n,x,1,wrk1,1) ; do i=1,n { sx(i) = i } - call ssort(wrk1,sx,n,2) - do i=1,n { isx(i) = ifix(sx(i)) } - - # Arrange x in increasing distinct order - # carry along c, scale x-vector - - call scopy(n,x,1,wrk1,1) ; call scopy(n,c,1,sc,1) - call ssort(wrk1,sc,n,2) - call scopy(n,wrk1,1,sx,1) - - deltl = (wrk1(n/10)-wrk1(1))/(.1*n-1.) - deltu = (wrk1(n)-wrk1((n/10)*9))/(.1*n-1.) - if(ind!=1) { a = wrk1(1) - deltl ; b = wrk1(n) + deltu } - do i=1,n { x(i) = (x(i)-a)/(b-a) ; wrk1(i) = (wrk1(i)-a)/(b-a)} - call scopy(n,wrk1,1,sx,1) - - - nx = 1 ; sx(nx) = sx(1) ; sc(nx) = sc(1) - - for(i=2;i<=n;i=i+1) - - { if(wrk1(i)>sx(nx)+eps) - { nx = nx + 1 ; sx(nx) = sx(i) ; sc(nx) = sc(i)} - else { sc(nx) = sc(nx)+sc(i) } - } - - - call sknotl(sx,nx,knot,k) ; nk = k-4 - if(ind==1) { do i=1,3 { knot(i)=a ; knot(k-i+1) = b }} - - - - # Smooth Estimate - - call icreg1(x,t,w,c,m,n,coef1,knot,nk, - smreg,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4,xw) - - # Histogram Estimate - - call icreg(x,t,m,c,n,sx,nx,tmed,coef0,knot,nk, - smreg,wrk1,wrk2,wrk3,wrk4) - - - nbvcd = 4*n+nk - call sbvc(knot,nk,sx,nx,bvcd,sd,ld,nd,wrk1) - - - - -return -end //GO.SYSIN DD bart/setup.r echo bart/sg.r 1>&2 sed >bart/sg.r <<'//GO.SYSIN DD bart/sg.r' 's/^-//' -subroutine sg(sg0,sg1,sg2,sg3,tb,nb,id) - - -integer nb,id,ileft,ilo,mflag, - i,ii,j,jj # indices - -real sg0(nb),sg1(nb),sg2(nb),sg3(nb),tb(nb+4), - fx(4,4),a(4),w(4),vx,vnikx(4,3),work(15), - wpt - - - #PURPOSE - -# -# Computes a Gram matrix from cubic B-splines -# - -# sgm[0-3](nb) Symmetric 7-banded matrix -# whose (i,j)'th element contains the integral of -# B (id) (i,.) B (id) (j,.) , i=1,2 ... nb and j=i,...nb. -# Only the upper four diagonals are computed. -# sg0 - main diagonal -# sg1 - next diagonal band above sg0 -# sg2 - next diagonal band above sg1 -# sg3 - next diagonal band above sg2 - -# -# (id=0,1,2) If id=2 then the Gram matrix corresponding to -# cubic smoothing splines is computed. - - - - -# Uses BSPVD and INTRV in the CMLIB - - - - - # Method involves a 4-point Gauss-Legendre quadrature formula - # to evaluate the integrals (this is exact for 7'th degree polynomials) - - - # Gaussian quadature points and weights ( standardized interval [-1,1] ) - - a(1) = -.861136 ; a(2) = -.339981 ; a(3) = .339981 ; a(4) = .861136 - w(1) = .347855 ; w(2) = .652145 ; w(3) = .652145 ; w(4) = .347855 - - #Initialise the sigma vectors - do i=1,nb{ sg0(i)=0.;sg1(i)=0.;sg2(i)=0.;sg3(i)=0.} - - ilo = 1 - - do i=1,nb { - wpt = (tb(i+1)-tb(i))/2. - - do j=1,4 { - vx = tb(i) + wpt*(1.+a(j)) - call intrv(tb(1),(nb+1),vx,ilo,ileft,mflag) - call bspvd (tb,4,id+1,vx,ileft,4,vnikx,work) - call scopy(4,vnikx(1,id+1),1,fx(1,j),1) } - -# Evaluate Contributions to the simga vectors - -if(ileft>=4){ -do ii=1,4 { - -jj=ii -sg0(ileft-4+ii) = sg0(ileft-4+ii) + - wpt*(w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4)) -jj=ii+1 -if(jj<=4) {sg1(ileft+ii-4) = sg1(ileft+ii-4) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4))} -jj=ii+2 -if(jj<=4) {sg2(ileft+ii-4) = sg2(ileft+ii-4) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4))} -jj=ii+3 -if(jj<=4) {sg3(ileft+ii-4) = sg3(ileft+ii-4) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4))} - } - } - -else if(ileft==3){ -do ii=1,3 { - -jj=ii -sg0(ileft-3+ii) = sg0(ileft-3+ii) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4)) -jj=ii+1 -if(jj<=3) {sg1(ileft+ii-3) = sg1(ileft+ii-3) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4))} -jj=ii+2 -if(jj<=3) {sg2(ileft+ii-3) = sg2(ileft+ii-3) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4))} - } - } - -else if(ileft==2){ -do ii=1,2 { - -jj=ii -sg0(ileft-2+ii) = sg0(ileft-2+ii) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4)) -jj=ii+1 -if(jj<=2) {sg1(ileft+ii-2) = sg1(ileft+ii-2) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4))} - } - } - -else if(ileft==1){ -do ii=1,1 { - -jj=ii -sg0(ileft-1+ii) = sg0(ileft-1+ii) + - wpt* (w(1)*fx(ii,1)*fx(jj,1)+w(2)*fx(ii,2)*fx(jj,2)+ - w(3)*fx(ii,3)*fx(jj,3)+w(4)*fx(ii,4)*fx(jj,4)) - } - }} - - -return -end //GO.SYSIN DD bart/sg.r echo bart/sgram.r 1>&2 sed >bart/sgram.r <<'//GO.SYSIN DD bart/sgram.r' 's/^-//' -subroutine sgram(sg0,sg1,sg2,sg3,tb,nb) - -real sg0(nb),sg1(nb),sg2(nb),sg3(nb),tb(nb+4), - vnikx(4,3),work(15),yw1(4),yw2(4), - wpt - -integer nb,ileft,ilo,mflag, - i,ii,jj # indices - - - #PURPOSE - -# Calculation of the cubic B-spline smoothness prior -# for "usual" interior knot setup. - - -# Uses BSPVD and INTRV in the CMLIB - - - - -# sgm[0-3](nb) Symmetric matrix -# whose (i,j)'th element contains the integral of -# B''(i,.) B''(j,.) , i=1,2 ... nb and j=i,...nb. -# Only the upper four diagonals are computed. - - - #Initialise the sigma vectors - do i=1,nb{ sg0(i)=0.;sg1(i)=0.;sg2(i)=0.;sg3(i)=0.} - - ilo = 1 - - do i=1,nb { - - # Calculate a linear approximation to the - # second derivative of the non-zero B-splines - # over the interval [tb(i),tb(i+1)]. - - call intrv(tb(1),(nb+1),tb(i),ilo,ileft,mflag) - - #Left end second derivatives - call bspvd (tb,4,3,tb(i),ileft,4,vnikx,work) - - # Put values into yw1 - do ii=1,4 { yw1(ii) = vnikx(ii,3) } - - - #Right end second derivatives - call bspvd (tb,4,3,tb(i+1),ileft,4,vnikx,work) - - # Slope*(length of interval) in Linear Approximation to B'' - do ii=1,4 { yw2(ii) = vnikx(ii,3) - yw1(ii) } - wpt = tb(i+1) - tb(i) - - -# Calculate Contributions to the simga vectors - -if(ileft>=4){ -do 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(jj<=4) {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 )} -jj=ii+2 -if(jj<=4) {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 )} -jj=ii+3 -if(jj<=4) {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 )} - - } - } - -else if(ileft==3){ -do 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(jj<=3) {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 )} -jj=ii+2 -if(jj<=3) {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 )} - - } - } - -else if(ileft==2){ -do 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(jj<=2) {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 )} - - } - } - -else if(ileft==1){ -do 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 ) - - } - }} - - -return -end //GO.SYSIN DD bart/sgram.r echo bart/shessc.r 1>&2 sed >bart/shessc.r <<'//GO.SYSIN DD bart/shessc.r' 's/^-//' -subroutine shessc(bvec,sd,ld,nd,nbvc,nb,f,cumhaz,sm,n,hs0,hs1,hs2,hs3) - -real bvec(nbvc),f(n),cumhaz(n),sm, - hs0(nb),hs1(nb),hs2(nb),hs3(nb) - -integer sd(nb),ld(nb),nd(nb),nbvc,nb,n,i,j,k,ne - - - do j=1,nb{ - - k=j ; ne = nd(j) - hs0(j) = 0. - for(i=1;i<=ne;i=i+1){ - hs0(j) = hs0(j) + f(ld(k)+i-1)*cumhaz(ld(k)+i-1)* - bvec(sd(k)+i-1)*bvec(sd(j)+(ld(k)-ld(j))+i-1) } - - hs0(j) = hs0(j)/sm - - k=j+1 ; if((j+1)<=nb) { ne = nd(j)- (ld(k)-ld(j)) - hs1(j) = 0. - for(i=1;i<=ne;i=i+1){ - hs1(j) = hs1(j) + f(ld(k)+i-1)*cumhaz(ld(k)+i-1)* - bvec(sd(k)+i-1)*bvec(sd(j)+(ld(k)-ld(j))+i-1) }} - - hs1(j) = hs1(j)/sm - - k=j+2 ; if((j+2)<=nb) { ne = nd(j)- (ld(k)-ld(j)) - hs2(j) = 0. - for(i=1;i<=ne;i=i+1){ - hs2(j) = hs2(j) + f(ld(k)+i-1)*cumhaz(ld(k)+i-1)* - bvec(sd(k)+i-1)*bvec(sd(j)+(ld(k)-ld(j))+i-1) }} - - hs2(j) = hs2(j)/sm - - k=j+3 ; if((j+3)<=nb) { ne = nd(j)- (ld(k)-ld(j)) - hs3(j) = 0. - for(i=1;i<=ne;i=i+1){ - hs3(j) = hs3(j) + f(ld(k)+i-1)*cumhaz(ld(k)+i-1)* - bvec(sd(k)+i-1)*bvec(sd(j)+(ld(k)-ld(j))+i-1) }} - - hs3(j) = hs3(j)/sm - - } - - -return -end - - - - -subroutine add(t,x,w,c,m,sm,n,knot,hess,ldk,nk,coef,rit) - -real t(n),x(n),w(n),c(n),m(n),sm, - knot(nk+6), - hess(ldk,nk), - coef(nk),rit(nk) - -integer n,ldk,nk - - - # Local - -real rito,bvalu,work(15) - -integer i,j,k,l,icoef - - -# Computes the additional Part of the Hessian in the Cox Model -# this is a slow routine. It might be worth creating a design -# matrix (n X nk) once and for all. - - - do j=1,nk { coef(j) = 0. - do k=j,nk { hess(k,j) = 0. }} - - - - - - i = 1 ; rito = 0. - for(l=i;l<=n;l=l+1){ rito = rito + m(l)*w(l)} - - do j=1,nk { coef(j) = 1. ; rit(j) = 0. ; icoef = 1. - - for(l=i;l<=n;l=l+1) - { rit(j) = rit(j) + m(l)*w(l)* - bvalu(knot,coef,nk,4,0,x(l),icoef,work)} - - coef(j) = 0. } - if(c(i)>0.) { - do j=1,nk{ - do k=j,nk{ hess(k,j) = hess(k,j) - - m(i)*(rit(j)/rito)*(rit(k)/rito)/sm }}} - - - for(i=2;i<=n;i=i+1){ - - if(c(i)>0.) { if(t(i)>t(i-1)) { - rito = 0. - for(l=i;l<=n;l=l+1){ rito = rito + m(l)*w(l)} - - do j=1,nk { coef(j) = 1. ; rit(j) = 0. ; icoef = 1. - - for(l=i;l<=n;l=l+1) - { rit(j) = rit(j) + m(l)*w(l)* - bvalu(knot,coef,nk,4,0,x(l),icoef,work)} - - coef(j) = 0. } - - } - do j=1,nk{ - do k=j,nk{ hess(k,j) = hess(k,j) - - m(i)*(rit(j)/rito)*(rit(k)/rito)/sm }}} - - } - - - -return -end //GO.SYSIN DD bart/shessc.r echo bart/shessd.r 1>&2 sed >bart/shessd.r <<'//GO.SYSIN DD bart/shessd.r' 's/^-//' -subroutine shessd(bvec,s,l,n,nbvc,nb,f,gwts,ngpt,hs0,hs1,hs2,hs3) - -real bvec(nbvc),gwts(ngpt),f(ngpt), - hs0(nb),hs1(nb),hs2(nb),hs3(nb), - gq3 - -integer s(nb),l(nb),n(nb),nbvc,nb,ngpt,j,k,ne - - - do j=1,nb{ - - k=j ; ne = n(j) - hs0(j) = gq3(f(l(k)),bvec(s(j)),bvec(s(k)),gwts(l(k)),ne) - - k=j+1 ; if((j+1)<=nb) { ne = n(j)- (l(k)-l(j)) - hs1(j) = gq3(f(l(k)),bvec(s(j)+(l(k)-l(j))), - bvec(s(k)),gwts(l(k)),ne) } - - k=j+2 ; if((j+2)<=nb) { ne = n(j)- (l(k)-l(j)) - hs2(j) = gq3(f(l(k)),bvec(s(j)+(l(k)-l(j))), - bvec(s(k)),gwts(l(k)),ne) } - - k=j+3 ; if((j+3)<=nb) { ne = n(j)- (l(k)-l(j)) - hs3(j) = gq3(f(l(k)),bvec(s(j)+(l(k)-l(j))), - bvec(s(k)),gwts(l(k)),ne) } - } - - -return -end //GO.SYSIN DD bart/shessd.r echo bart/shessh.r 1>&2 sed >bart/shessh.r <<'//GO.SYSIN DD bart/shessh.r' 's/^-//' -subroutine shessh(bvec,s,l,n,nbvc,nb,f,g,gwts,ngpt,hs0,hs1,hs2,hs3) - -integer s(nb),l(nb),n(nb),nbvc,nb,ngpt,j,k,ne - -real bvec(nbvc),gwts(ngpt),f(ngpt),g(ngpt), - hs0(nb),hs1(nb),hs2(nb),hs3(nb), - gq4 - - - do j=1,nb{ - - k=j ; ne = n(j) - hs0(j) = gq4(f(l(k)),g(l(k)),bvec(s(j)), - bvec(s(k)),gwts(l(k)),ne) - - k=j+1 ; if((j+1)<=nb) { ne = n(j)- (l(k)-l(j)) - hs1(j) = gq4(f(l(k)),g(l(k)),bvec(s(j)+(l(k)-l(j))), - bvec(s(k)),gwts(l(k)),ne) } - - k=j+2 ; if((j+2)<=nb) { ne = n(j)- (l(k)-l(j)) - hs2(j) = gq4(f(l(k)),g(l(k)),bvec(s(j)+(l(k)-l(j))), - bvec(s(k)),gwts(l(k)),ne) } - - k=j+3 ; if((j+3)<=nb) { ne = n(j)- (l(k)-l(j)) - hs3(j) = gq4(f(l(k)),g(l(k)),bvec(s(j)+(l(k)-l(j))), - bvec(s(k)),gwts(l(k)),ne) } - } - - -return -end //GO.SYSIN DD bart/shessh.r echo bart/shist.r 1>&2 sed >bart/shist.r <<'//GO.SYSIN DD bart/shist.r' 's/^-//' -subroutine shist(x,y,w,n,histx,histy,knot,xw,min,range,nx,nk) - -real x(n),y(n),w(n), - min,range, - histx(n+2),histy(n+2), - knot(n+6),xw(n) - - -integer n,nx,nk - - - # Local - -real eps,con,xh,yh,intg1,bvalu,xv,work(15) -integer i,icoef,indk,j,js,k,klow,khigh,nxt,nh,min0,max0,ycnt - - - # Sort x's in increasing order - - call scopy(n,x,1,xw,1) - call ssort(x,w,n,2) - call ssort(xw,y,n,2) - - range = x(n)-x(1) ; min = x(1) ; eps = .1e-9 - do i=1,n { x(i) = (x(i)-min)/range } - call scopy(n,x,1,xw,1) - - - nx = 1 ; x(nx) = x(1) ; w(nx) = w(1) ; y(nx) = y(1) ; ycnt = 1 - - for(i=2;i<=n;i=i+1) - - { if(xw(i)>x(nx)+eps) - { y(nx) = y(nx)/ycnt - ycnt = 1 - nx = nx + 1 - x(nx) = x(i) - y(nx) = y(i) - w(nx) = w(i) } - else - { y(nx) = y(nx)+y(i) - ycnt = ycnt + 1 - w(nx) = w(nx) + w(i) } - } - - - # Put knots at distinct knot points - - call sknotl(x,nx,knot,k) ; nk = k-4 - - - - - - - # Form histx, histy ( A bit rough ) - - # nh = # of bins = (5./sd)*nx**(1/5) - - nh = max0((5.0/(.74*(x((nx/4)*3)-x(nx/4))))*nx**(.2),1) - nh = min0(nh,nx) - - # Cumulative - - xw(1) = w(1) - do i=2,nx { xw(i) = xw(i-1) + w(i) } - con = xw(nx) - do i=1,nx { xw(i) = xw(i)/con } - - # Histogram estimates of x-density and regression - - klow = 1 - # Find indk s.t. knot(indk+3)>=x(klow) - do j=1,nk-2{ - if(knot(j+3)>=x(klow)) {indk=j;break} } - nxt = indk - 1 ; js = indk - - - do k=2,nh { khigh = 1 + (k-1)*(nx-1)/(nh-1) - - xh = (xw(khigh)-xw(klow))/(x(khigh)-x(klow)) - - yh = 0. ; con = 0. - - for(i=klow;i<=khigh;i=i+1) - {yh = yh + y(i)*w(i) ; con = con + w(i) } - - yh = yh/con - - for(j=indk;j<=nk-2;j=j+1) - - { if(x(klow)<= knot(j+3) & knot(j+3) < x(khigh) ) - - { histx(j+1) = xh ; histy(j+1) = yh ; nxt = nxt + 1 } - - else { break }} - - indk = nxt + 1 ; klow = khigh } - - # Fix up the begining and ends - - for(j=1;j<=js-1;j=j+1) { histx(j+1) = histx(js+1) - histy(j+1) = histy(js+1)} - for(j=nxt;j<=nk-2;j=j+1){histx(j+1) = histx(nxt) - histy(j+1) = histy(nxt) } - - histx(1) = histx(2) - histx(nk) = histx(nk-1) - histy(1) = histy(2) - histy(nk) = histy(nk-1) - - - - # Evaluate and Normalize the density estimate - - icoef = 1 ; do i=1,nx { xv = x(i) - xw(i) = bvalu(knot,histx,nk,4,0,xv,icoef,work(1))} - - con = intg1(xw,x,nx) - do i=1,nk { histx(i) = histx(i)/con } - - -return -end //GO.SYSIN DD bart/shist.r echo bart/sinerp.r 1>&2 sed >bart/sinerp.r <<'//GO.SYSIN DD bart/sinerp.r' 's/^-//' -subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag) - -real abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk), - wjm3(3),wjm2(2),wjm1(1),c0,c1,c2,c3 - -integer flag,ld4,nk,ldnk,i,j,k - - - - - # Purpose : Computes Inner Products between columns of L(-1) - # where L = abd is a Banded Matrix with 3 subdiagonals - - # The algorithm works in two passes: - # - # Pass 1 computes (cj,ck) k=j,j-1,j-2,j-3 ,j=nk, .. 1 - # Pass 2 computes (cj,ck) k<=j-4 (If flag == 1 ). - # - # A refinement of Elden's trick is used. - # - - - - - # Pass 1 - - - wjm3(1)=0. ; wjm3(2)=0. ; wjm3(1)=0. - wjm2(1)=0. ; wjm2(2)=0. - wjm1(1)=0. - - do i=1,nk { j=nk-i+1 - - c0 = 1./abd(4,j) - if(j<=nk-3) {c1 = abd(1,j+3)*c0 - c2 = abd(2,j+2)*c0 - c3 = abd(3,j+1)*c0 } - else if(j==nk-2) {c1 = 0. - c2 = abd(2,j+2)*c0 - c3 = abd(3,j+1)*c0 } - else if(j==nk-1) {c1 = 0. - c2 = 0. - c3 = abd(3,j+1)*c0 } - else if(j==nk) {c1 = 0. - c2 = 0. - c3 = 0.} - - - 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) - - } - - - if(flag==0) {return} - - - # Pass 2 - - - else { # Compute p2ip - - do i=1,nk { j=nk-i+1 - for(k=1;k<=4 & j+k-1<=nk;k=k+1) - { p2ip(j,j+k-1) = p1ip(5-k,j) }} - - - do i=1,nk { j=nk-i+1 - for(k=j-4;k>=1;k=k-1){ - - 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) ) } - } - - return - - } - - -end //GO.SYSIN DD bart/sinerp.r echo bart/sknotl.r 1>&2 sed >bart/sknotl.r <<'//GO.SYSIN DD bart/sknotl.r' 's/^-//' -subroutine sknotl(x,n,knot,k) - -real x(n),knot(n+6),a1,a2,a3,a4 -integer n,k,ndk,j - - - # Allocate knots acording to the cutoffs given below - - - # Cutoff constants - - a1 = log(50.)/log(2.) ; a2 = log(100.)/log(2.) - a3 = log(140.)/log(2.) ; a4 = log(200.)/log(2.) - - # Cutoff Criteria - - if(n<50) { ndk = n } - else if (n>=50 & n<200) { ndk = 2.**(a1+(a2-a1)*(n-50.)/150.) } - else if (n>=200 & n<800) { ndk = 2.**(a2+(a3-a2)*(n-200.)/600.) } - else if (n>=800 & n<3200) { ndk = 2.**(a3+(a4-a3)*(n-800.)/2400.) } - else if (n>=3200) { ndk = 200. + (n-3200)**.2 } - - k = ndk + 6 - - - # Allocate Knots ( note no account is taken of any weighting vector ) - - do j=1,3 { knot(j) = x(1) } - do j=1,ndk { knot(j+3) = x( 1 + (j-1)*(n-1)/(ndk-1) ) } - do j=1,3 { knot(ndk+3+j) = x(n) } - -return -end //GO.SYSIN DD bart/sknotl.r echo bart/slap.r 1>&2 sed >bart/slap.r <<'//GO.SYSIN DD bart/slap.r' 's/^-//' -subroutine slap(f,z,us,lev,n, - crit,icrit,spar,ispar,lspar,uspar,tol, - sbd,ld,wf,we,wx,clev,knot,ns,ier) - - - # A Discrete Biharmonic Smoother - -# -# The algorithm solves: -# -# [ B + del*f ] u = del*f *z -# -# on the uniform grid (i-1)/(n-1) , i=1,2, ... ,n+2 -# -# del is a function of the spar which is assumed to be between -# 0 and 1 - - - # Input - -# n+2 number of grid points -# f(n+2) x-density -# z(n+2) z-averages - - -# spar smoothing parameter -# ispar indicator saying if spar is supplied or to be estimated -# lspar, uspar lower and upper values for spar 0.,1. are good values -# tol used in Golden Search routine - - -# icrit indicator saying which cross validation score -# is to be computed - -# ld the leading dimension of sbd (ie ld4=4) - - - # Output - -# us(n+2) smoothed values -# lev(n+2) vector of leverages -# crit either ordinary of generalized CV score -# ier error indicator -# ier = 0 ___ everything fine -# ier = 1 ___ spar too small or too big -# problem in cholesky decomposition - - - - # Working arrays/matrix - -# we,wf,wx, odd jobs -# clev(ns+2) -# knot(ns+6) Used in the computation of the cv-approximation -# sbd(ld,n+2) Standard 6-point approximation to B in banded form - - -real f(n+2),z(n+2), - us(n+2),lev(n+2), - crit,spar,lspar,uspar,tol, - knot(ns+6),clev(ns+2), - we(n+2),wf(n+2),wx(n+2), - sbd(ld,n+2) - - -integer n,ns,icrit,ispar,ld,ier - - - - - # Local variables - -real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, - fu,fv,fw,fx,x, - ax,bx - - - - - - - - # Compute estimate - - - if(ispar==1) { # Value of spar supplied - - call sslvlp(f,z,us,lev,n, - knot,clev,ns, - sbd,ld, - spar,crit,icrit,ier, - we,wf,wx) - - return } - - else { - - # Use Forsythe Malcom and Moler routine to minimise criterion - - ax=lspar ; bx=uspar # f denotes the value of the criterion - -# -# an approximation x to the point where f attains a minimum on -# the interval (ax,bx) is determined. -# -# -# input.. -# -# ax left endpoint of initial interval -# bx right endpoint of initial interval -# f function subprogram which evaluates f(x) for any x -# in the interval (ax,bx) -# tol desired length of the interval of uncertainty of the final -# result ( .ge. 0.0) -# -# -# output.. -# -# fmin abcissa approximating the point where f attains a minimum -# -# -# the method used is a combination of golden section search and -# successive parabolic interpolation. convergence is never much slower -# than that for a fibonacci search. if f has a continuous second -# derivative which is positive at the minimum (which is not at ax or -# bx), then convergence is superlinear, and usually of the order of -# about 1.324.... -# the function f is never evaluated at two points closer together -# than eps*abs(fmin) + (tol/3), where eps is approximately the square -# root of the relative machine precision. if f is a unimodal -# function and the computed values of f are always unimodal when -# separated by at least eps*abs(x) + (tol/3), then fmin approximates -# the abcissa of the global minimum of f on the interval ax,bx with -# an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, -# then fmin may approximate a local, but perhaps non-global, minimum to -# the same accuracy. -# this function subprogram is a slightly modified version of the -# algol 60 procedure localmin given in richard brent, algorithms for -# minimization without derivatives, prentice - hall, inc. (1973). -# -# -# real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w -# real fu,fv,fw,fx,x -# -# c is the squared inverse of the golden ratio -# - c = 0.5*(3. - sqrt(5.0)) -# -# eps is approximately the square root of the relative machine -# precision. -# - eps = 1.00 - 10 eps = eps/2.00 - tol1 = 1.0 + eps - if (tol1 .gt. 1.00) go to 10 - eps = sqrt(eps) -# -# initialization -# - a = ax - b = bx - v = a + c*(b - a) - w = v - x = v - e = 0.0 - - spar = x - call sslvlp(f,z,us,lev,n, - knot,clev,ns, - sbd,ld, - spar,crit,icrit,ier, - we,wf,wx) - - fx = crit - fv = fx - fw = fx -# -# main loop starts here -# - 20 xm = 0.5*(a + b) - tol1 = eps*abs(x) + tol/3.0 - tol2 = 2.0*tol1 -# -# check stopping criterion -# - if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 -# -# is golden-section necessary -# - if (abs(e) .le. tol1) go to 40 -# -# fit parabola -# - r = (x - w)*(fx - fv) - q = (x - v)*(fx - fw) - p = (x - v)*q - (x - w)*r - q = 2.00*(q - r) - if (q .gt. 0.0) p = -p - q = abs(q) - r = e - e = d -# -# is parabola acceptable -# - 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 - if (p .le. q*(a - x)) go to 40 - if (p .ge. q*(b - x)) go to 40 -# -# a parabolic interpolation step -# - d = p/q - u = x + d -# -# f must not be evaluated too close to ax or bx -# - if ((u - a) .lt. tol2) d = sign(tol1, xm - x) - if ((b - u) .lt. tol2) d = sign(tol1, xm - x) - go to 50 -# -# a golden-section step -# - 40 if (x .ge. xm) e = a - x - if (x .lt. xm) e = b - x - d = c*e -# -# f must not be evaluated too close to x -# - 50 if (abs(d) .ge. tol1) u = x + d - if (abs(d) .lt. tol1) u = x + sign(tol1, d) - - spar = u - call sslvlp(f,z,us,lev,n, - knot,clev,ns, - sbd,ld, - spar,crit,icrit,ier, - we,wf,wx) - - fu = crit -# -# update a, b, v, w, and x -# - if (fu .gt. fx) go to 60 - if (u .ge. x) a = x - if (u .lt. x) b = x - v = w - fv = fw - w = x - fw = fx - x = u - fx = fu - go to 20 - 60 if (u .lt. x) a = u - if (u .ge. x) b = u - if (fu .le. fw) go to 70 - if (w .eq. x) go to 70 - if (fu .le. fv) go to 80 - if (v .eq. x) go to 80 - if (v .eq. w) go to 80 - go to 20 - 70 v = w - fv = fw - w = u - fw = fu - go to 20 - 80 v = u - fv = fu - go to 20 -# -# end of main loop -# - 90 continue ; spar = x ; crit = fx - return } - - - - - -return -end //GO.SYSIN DD bart/slap.r echo bart/slvcxp.r 1>&2 sed >bart/slvcxp.r <<'//GO.SYSIN DD bart/slvcxp.r' 's/^-//' -subroutine slvcxp(x,typx,ldx,n,p,t,c,m,w,n,beta,smreg,cumhaz, - wrk0,wrk1,wrk2,wrk3,wrk4,beta0,grad, - hess,svdhp,d,ldp,ierr) - - -real x(ldx,p),typx(p),t(n),c(n),m(n),w(n), - beta(p),smreg(n),cumhaz(n), - wrk0(n+2),wrk1(n),wrk2(n+6),wrk3(n),wrk4(n+2), - beta0(p),grad(p),hess(ldp,p), - svdhp(ldp,p),d(p) - -integer ldx,p,n,ldp,ierr - - - - # Local Variables - -real le,lo,ls,sm,sdot,tol,cchk,abs,sum -integer i,iter,j,max,k - - - - max = 15 ; tol = .1e-6 - cchk = tol+10. ; ierr = 0 ; iter = 0 - - - - # - # Purpose : Iteratively computes a p-dimensional Cox Regression Model - # - - - # Put the proportionality estimate into smreg - - do i=1,n {smreg(i) = exp(sdot(p,beta,1,x(i,1),ldx))*w(i) } - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - - # Starting Value for -ve Log Likelihood - - sm = 0. ; do i=1,n { sm = sm + m(i)} - - call evalle(t,c,m,sm,smreg,n,le) - lo = le ; ls = le - -# -# Newton's Method with a -# Trust Region Implementation of a Levenberg-Marquard Line Search -# -# beta = beta0 - [ H + alpha ] (-1) dir -# -# beta0 -> beta0 ; dir -> grad -# - - - while(iter tol) { - - iter = iter + 1 - do j=1,p { beta0(j) = beta(j) } - - # Find Newton direction for new estimate - - # Put the regression estimate into smreg - do i=1,n {smreg(i) = exp(sdot(p,beta,1,x(i,1),ldx))*w(i) } - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - call grhss(t,c,m,sm,smreg,cumhaz,n,x,ldx,grad,hess,ldp,p) - - - # Solve for the Update - - do j=1,p { do k=j,p { hess(k,j) = hess(j,k) }} - call ssvdc(hess,ldp,p,p,d,wrk1,svdhp,ldp,hess,ldp,wrk0,20,ierr) - if(ierr!=0) {write(6,'("Fucked")') ; return} - do j=1,p { beta(j) = sdot(p,grad,1,svdhp(1,j),1)/d(j) } - - call scopy(p,beta,1,wrk0,1) - do j=1,p { beta(j) = sdot(p,wrk0,1,svdhp(j,1),ldp) } - - sum = 0. - do j=1,p { sum = sum + abs(beta(j)*typx(j)) } - if(sum > 14. ) { do j=1,p { beta(j) = beta(j)*14./sum }} - - do j=1,p { beta(j) = beta0(j) - beta(j) } - - - # Find a new estimate in the direction of beta - - - - call search(x,ldx,t,c,m,sm,w,n, - beta,beta0,grad,p, - svdhp,d,wrk1,ldp, - smreg,le,lo,ierr,wrk0) - - - - cchk = abs(le-lo)/abs(le) - lo = le - - } - - -return -end //GO.SYSIN DD bart/slvcxp.r echo bart/smv5.r 1>&2 sed >bart/smv5.r <<'//GO.SYSIN DD bart/smv5.r' 's/^-//' -subroutine smv5(x,y,n) - - # Five-Point Weighted moving average - # The answer is put into y - -real x(n),y(n) -integer n - - # Local - -real w(5),v(5),val -integer j,jj - - - if(n<=7) {write(6,'("Error in mv5; too few data points")') - stop } - - - -jj = 1 ;do j=1,5 { w(j)=(x(6)-x(j))/(x(6)-x(jj)) - v(j)= y(j) } -y(jj) = val(v,w,5) - -jj = 2 ;do j=3,5 { w(j)=(x(6)-x(j))/(x(6)-x(jj)) } - w(1)= w(2) ; w(2) = 1. -y(jj) = val(v,w,5) - -jj = 3 ;do j=4,5 { w(j)=(x(6)-x(j))/(x(6)-x(jj)) } - w(1)= w(5) ; w(2) = w(4) ; w(3) = 1. -y(jj) = val(v,w,5) - - -do jj=4,n-3 { - do j=1,2 { w(j)=(x(jj)-x(jj-3+j))/(x(jj)-x(jj-3)) - w(j+3)=(x(jj+3)-x(jj+j))/(x(jj+3)-x(jj))} - - v(1)=v(2);v(2)=v(3);v(3)=v(4);v(4)=v(5);v(5)=y(jj+2) - -y(jj) = val(v,w,5) } - - - v(1)=v(2);v(2)=v(3);v(3)=v(4);v(4)=v(5);v(5)=y(n) - -jj = n-2 ;do j=1,2 { w(j)=(x(jj)-x(n-3+j))/(x(jj)-x(n-3)) } - w(4)= w(1) ; w(5) = w(2) -y(jj) = val(v,w,5) - -jj = n-1 ;do j=1,3 { w(j)=(x(jj)-x(n-3+j))/(x(jj)-x(n-3)) } - w(4)= 1. ; w(5) = w(1) -y(jj) = val(v,w,5) - -jj = n ;do j=1,4 { w(j)=(x(jj)-x(n-3+j))/(x(jj)-x(n-3)) } - w(5)= 1. -y(jj) = val(v,w,5) - - -return -end - - -real function val(v,w,n) - -real v(n),w(n),svw,sw,sml,lrg -integer n,i - - - - sml = v(1)*w(1) ; lrg = v(1)*w(1) - smlw = w(1) ; lrgw = w(1) - - for(i=2;i<=n;i=i+1) { if(sml>v(i)*w(i)) {sml=v(i)*w(i);smlw=w(i)} - if(lrg&2 sed >bart/srchcx.r <<'//GO.SYSIN DD bart/srchcx.r' 's/^-//' -subroutine srchcx(x,t,c,m,sm,w,n,cmin,cmax,coef,coef0,dir,nk, - smreg,knot,sg0,sg1,sg2,sg3,lambda,pl,plo,wrk0) - - -real x(n),t(n),c(n),m(n),sm,w(n), - cmin,cmax,coef(nk),coef0(nk),dir(nk), - smreg(n),knot(nk+4), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - lambda,pl,plo,wrk0(nk) - -integer n,nk - - - # Local Variables - -real work(15),bvalu,corr,plm,le, - exp,alpha,alpham,tau,snrm2 - -integer i,icoef,j,k - - - - - # - # A Step Size Halving Line Search Method - # - # coef = coef0 - alpha * (newton dir) - # - # where the values of alpha are choosen - # so that the step sizes - # - # norm[coef-coef0]/(norm[hcoef]+.1) - # - # are 2**(-k) , k = 1,2, ... 15. - # - # - - - - # Adjust Coeficient so that it's not too far from coef0 - - tau = snrm2(nk,dir,1)/(cmax-cmin+.1) - do j=1,nk { dir(j) = (dir(j)/tau)*5. } - - - - # Step Size Halving - - - plm = plo ; alpha = 0. ; alpham = 0. - - - do k=1,15 { alpha = .5**(k-1.) - - do j=1,nk {coef(j) = coef0(j) + dir(j)*alpha } ; icoef = 1 - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk,4,0,x(i),icoef,work))*w(i)} - - call evalle(t,c,m,sm,smreg,n,le) - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk0) - pl = le + lambda*sdot(nk,coef,1,wrk0,1) - - if(pl&2 sed >bart/srchdn.r <<'//GO.SYSIN DD bart/srchdn.r' 's/^-//' -subroutine srchdn(w,n,coef,coef0,cmin,cmax, - knot,dir,nk,sg0,sg1,sg2,sg3, - smden,grid,qwts,ngpt,bvcd,sd,nd,ld,nbvcd, - lambda,pl,plo,wrk) - - -integer n,nk,ngpt, - sd(nk),nd(nk),ld(nk),nbvcd, - - j,k - - -real w(n),coef(nk),coef0(nk),cmin,cmax, - knot(nk+4),dir(nk),sg0(nk),sg1(nk),sg2(nk),sg3(nk), - smden(ngpt),grid(ngpt),qwts(ngpt),bvcd(nbvcd), - lambda,pl,plo,wrk(nk), - - alpha,alpham,le,plm, - sdot,snrm2,tau,amax1 - - - - # - # A Step Size Halving Line Search Method - # - # coef = coef0 - alpha * (newton dir) - # - # where the values of alpha are choosen - # so that the step sizes - # - # norm[coef-coef0]/(norm[hcoef]+.1) - # - # are 3**(-k) , k = 1,2, ... 5. - # - # - - # Adjust Coeficient so that it's not too far from coef0 - - tau = amax1(1.,snrm2(nk,dir,1)/(cmax-cmin+.1)) - - do j=1,nk { dir(j) = dir(j)/tau } - - do j=1,nk { coef(j) = coef0(j) + dir(j) } - - - # Evaluate a Normalized Newton Step and if it's ok return - - call evlkdn(w,n,coef,nk,smden, - knot,grid,qwts,ngpt,bvcd,sd,nd,ld,nbvcd,le) - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk) - - pl = le + lambda*sdot(nk,coef,1,wrk,1) - - if(pl/plo<.95) {return} - - - - # Else do Step Size Halving - - - if(pl&2 sed >bart/srchhz.r <<'//GO.SYSIN DD bart/srchhz.r' 's/^-//' -subroutine srchhz(c,rn,n,cmin,cmax,norm,coef,coef0,dir,nk, - sg0,sg1,sg2,sg3, - smhaz,knot, - grid,qwts,ngpt, - bvcd,sd,nd,ld,nbvcd, - lambda,pl,plo,wrk) - - -integer n,nk,ngpt,sd(nk),nd(nk),ld(nk),nbvcd, - - j,k - - -real c(n),rn(ngpt),cmin,cmax,norm, - coef(nk),coef0(nk),dir(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - smhaz(ngpt),knot(nk+4), - grid(ngpt),qwts(ngpt), - bvcd(nbvcd), - lambda,pl,plo,wrk(nk), - - plm,alpham,alpha,le,sdot,snrm2,tau,amax1 - - - # - # A Step Size Halving Line Search Method - # - # coef = coef0 - alpha * (newton dir) - # - # where the values of alpha are choosen - # so that the step sizes - # - # norm[coef-coef0]/(norm[hcoef]+.1) - # - # are 2**(-k) , k = 1,2, ... 6. - # - # - - # Adjust Coeficient so that it's not too far from coef0 - - tau = amax1(1.,snrm2(nk,dir,1)/(cmax-cmin+.1)) - - do j=1,nk { dir(j) = dir(j)/tau } - - do j=1,nk { coef(j) = coef0(j) + dir(j) } - - - # Evaluate a Normalized Newton Step and if it's ok return - - call evlkhz(rn,c,n,coef,nk,smhaz,knot,grid,qwts,ngpt, - bvcd,sd,nd,ld,nbvcd,le) - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk) - - pl = le + lambda*sdot(nk,coef,1,wrk,1) - - if(pl/plo<.95) {return} - - - # Else do Step Size Halving - - - if(pl&2 sed >bart/sslvcx.r <<'//GO.SYSIN DD bart/sslvcx.r' 's/^-//' -subroutine sslvcx(x,t,c,m,w,n, - sx,sc,isx,nsx, - bvcd,nbvcd,sd,nd,ld, - coef0,coef1,hcoef,coef,knot,nk,icrit,crit,spar, - smreg,cumhaz,cvlik,lev, - wrk0,wrk1,wrk2,wrk3,wrk4,wrk5,wrk6,nw, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,hess,ld4,ldnk,ierr) - - -real x(n),t(n),c(n),m(n),w(n), - sx(nsx),sc(nsx), - bvcd(nbvcd), - coef0(nk),coef1(nk),hcoef(nk),coef(nk),knot(nk+4), - crit,spar, - smreg(n),cumhaz(n),cvlik(n),lev(n), - wrk0(nw),wrk1(nw),wrk2(nw),wrk3(nw),wrk4(nw), - wrk5(nw),wrk6(nw), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - abd(ld4,nk),p1ip(ld4,nk),hess(ldnk,nk) - - -integer isx(n),nsx,n,nk,nbvcd,nw, - sd(nk),ld(nk),nd(nk), - icrit,ld4,ldnk,ierr - - - # Local Variables - -real le,pl,plo,bvalu,sm,dfr,t1,t2,lamr,dof,d0,d1,cmin,cmax, - lambda,sdot,tol,cchk - -integer i,iter,icoef,j,k,max - - - - max = 15; tol = .1e-4 - cchk = tol+10. ; ierr = 0 ; iter = 0 - - # Purpose : Iteratively computes a Cox Regression - - - # Set Initial Estimate and lambda - - # Initial Estimate - do j=1,nk { hcoef(j) = (1.-spar)*coef0(j)+spar*coef1(j) } - cmin=hcoef(1) ; cmax=hcoef(1) - do j=2,nk { if(cmaxhcoef(j)) { cmin=hcoef(j) } } - - icoef = 1 - do i=1,n { smreg(i) = exp(bvalu(knot,hcoef,nk, - 4,0,x(i),icoef,wrk0))*w(i) } - call breslw(t,c,m,smreg,n,cumhaz,wrk0,wrk1,wrk2,wrk3,wrk4) - - sm = 0. ; do i=1,n { sm = sm + m(i) } - call cxghs(t,x,c,m,sm,smreg,n,0, - knot,wrk6,hs0,hs1,hs2,hs3,p2ip,ldnk,nk, - 0.,coef,sg0,sg1,sg2,sg3, - wrk0,wrk1,wrk2,wrk3,wrk4) - - t1 = 0. ; t2 = 0. - do i=2,nk-1 {t1 = t1 + hs0(i) ; t2 = t2 + 2.*sg0(i) } - lamr = t1/t2 - - call df(dof,1.,lamr,x,t,c,m,w,sm,n,knot,nk,smreg,cumhaz, - lev,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,wrk0,nw,ierr) - d0 = .75 ; if(nk>60) {d1 = nk/5. -2. } - else {d1 = nk -2. } - d0 = 1./d0 ; d1 = 1./d1 - lambda = (dof-2.)*(d0*d1)/(spar*d1+(1.-spar)*d0) - - lambda = (lambda)**4 * lamr - - do j=1,nk { coef(j) = hcoef(j) } - call evalle(t,c,m,sm,smreg,n,le) - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk0) - - plo = le + lambda*sdot(nk,coef,1,wrk0,1) - - -# -# Newton's Method with a Step Size Halving Strategy -# -# coef = coef0 - (coef-coef0)*alpha -# -# coef0 -> wrk5 ; dir -> wrk6 -# - - - while(iter tol) { - - - iter = iter + 1 - call scopy(nk,coef,1,wrk5,1) - - # Find Newton direction for new estimate - - # Put the regression estimate into smreg - - icoef = 1 - do i=1,n {smreg(i) = exp(bvalu(knot,coef,nk, - 4,0,x(i),icoef,wrk1))*w(i)} - - call cxghs(t,x,c,m,sm,smreg,n,1, - knot,wrk6,hs0,hs1,hs2,hs3,hess,ldnk,nk, - lambda,wrk5,sg0,sg1,sg2,sg3, - wrk0,wrk1,wrk2,wrk3,wrk4) - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - do j=1,nk { hess(j,j) = abd(4,j) + hess(j,j) } - do j=1,nk-1 { hess(j+1,j) = abd(3,j+1) + hess(j+1,j) } - do j=1,nk-2 { hess(j+2,j) = abd(2,j+2) + hess(j+2,j) } - do j=1,nk-3 { hess(j+3,j) = abd(1,j+3) + hess(j+3,j) } - do j=1,nk { do k=j,nk { hess(j,k) = hess(k,j) }} - - call pcg(hess,abd,wrk6,coef,ldnk,ld4,nk,.000001, - wrk0,wrk1,wrk2,wrk3,wrk4,ierr) - - if(ierr!=0) {go to 10} - - do j=1,nk { wrk6(j) = - coef(j) } - - # Find a new estimate in the direction of wrk6 - # On return smreg contains the estimated proportionality factor - - call srchcx(x,t,c,m,sm,w,n,cmin,cmax,coef,wrk5,wrk6,nk, - smreg,knot,sg0,sg1,sg2,sg3,lambda,pl,plo,wrk0) - - cchk = abs(pl-plo)/pl - plo = pl - - - - } - -10 continue - - - if(ierr > 0 ) {crit=10. ; return} - - - # Compute the criterion function if requested - - if(icrit==0) { return} - - else { # Compute Linearized CV score - - call gxv(x,t,c,m,w,sm,n,coef,knot,nk,smreg,cumhaz, - crit,lev,cvlik,lambda,hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3, - abd,p1ip,hess,ld4,ldnk,wrk0,wrk1,wrk2,wrk3,wrk4,nw) - - dfr = 0. ; do i=1,n { dfr = dfr + 1. - lev(i) } - - crit = crit/(dfr**2) - - } - - - -return -end //GO.SYSIN DD bart/sslvcx.r echo bart/sslvdn.r 1>&2 sed >bart/sslvdn.r <<'//GO.SYSIN DD bart/sslvdn.r' 's/^-//' -subroutine sslvdn(x,w,n,knot,nk, - coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcgr,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev,wts, - crit,dcrit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ierr) - -integer n,nk,ngpt,nbvcg,nbvcd, - sd(nk),ld(nk),nd(nk), - sg(nk),lg(nk),ng(nk), - dcrit,icrit,ld4,ldnk,ierr, - - i,iter,j,max - - -real x(n),w(n), - knot(nk+4), - coef0(nk),coef1(nk),hcoef(nk), - grid(ngpt),qwts(ngpt), - bvcgr(nbvcg),bvcd(nbvcd), - coef(nk),smden(ngpt),cvden(n), - lev(ngpt),wts(ngpt), - crit(dcrit),spar, - 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), - wrk1(nk),wrk2(nk),wrk3(nk), - - le,pl,plo,pls,t1,t2,lamr,dof,d0,d1,cmin,cmax, - lambda,sdot,tol,cchk,sqrt - - - - - - max = 15; tol = .1e-3 ; cchk = tol+10. - - - # Purpose : Solves the smoothing problem and computes - # A linearized CV score - - # Set Initial Estimate and lambda - - # Initial Estimate - do j=1,nk { hcoef(j) = (1.-spar)*coef0(j)+spar*coef1(j) } - cmin=hcoef(1) ; cmax=hcoef(1) - do j=2,nk { if(cmaxhcoef(j)) { cmin=hcoef(j) } } - - call evlkdn(w,n,hcoef,nk,smden,knot,grid,qwts,ngpt, - bvcd,sd,nd,ld,nbvcd,le) - call shessd(bvcgr,sg,lg,ng,nbvcg,nk,smden, - qwts,ngpt,hs0,hs1,hs2,hs3) - do i=1,ngpt { wts(i) = sqrt(qwts(i)*smden(i)) } - - t1 = 0. ; t2 = 0. - do i=2,nk-1 {t1 = t1 + hs0(i) ; t2 = t2 + 2.*sg0(i) } - lamr = t1/t2 - - call dfr(dof,1.,lamr,grid,wts,ngpt,knot,nk,lev, - hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,ld4,ierr) - - d0 = .75 ; if(nk>60) {d1 = nk/5. -2. } - else {d1 = nk -2. } - d0 = 1./d0 ; d1 = 1./d1 - lambda = (dof-2.)*(d0*d1)/(spar*d1+(1.-spar)*d0) - - lambda = (lambda)**4 * lamr - - do j=1,nk { coef(j) = hcoef(j) } - - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk2) - - # Starting Value for Penalized Likelihood - - pls = le + lambda*sdot(nk,coef,1,wrk2,1) - plo = pls - - - - - - # Use a damped Newton's method to compute the estimate - # - # coef = coef0 - tau*[H+2*lambda*OMEGA] (-1) dir - # - # coef0 -> wrk1 ; dir -> wrk2 - # - - - - - iter = 0 - - while(iter tol) { - - iter = iter + 1 - do j=1,nk { wrk1(j) = coef(j) } - - - call shessd(bvcgr,sg,lg,ng,nbvcg, - nk,smden,qwts,ngpt,hs0,hs1,hs2,hs3) - - call sdirdn(w,n,bvcgr,sg,lg,ng,nbvcg,nk,qwts,ngpt, - bvcd,sd,ld,nd,nbvcd, - smden,lambda,sg0,sg1,sg2,sg3, - wrk1,wrk2) - - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) - if(ierr.ne.0) { return } - call spbsl(abd,ld4,nk,3,wrk2) - do j=1,nk { wrk2(j) = 0. - wrk2(j) } - - # Find new estimate in the direction of coef - - - call srchdn(w,n,coef,wrk1,cmin,cmax, - knot,wrk2,nk,sg0,sg1,sg2,sg3, - smden,grid,qwts,ngpt,bvcd,sd,nd,ld,nbvcd, - lambda,pl,plo,wrk3) - if(pl>plo) { plo = pl } - cchk = abs(pl-plo)/abs(pl) ; plo = pl - - } - - - - if(pls<=pl .or. ierr > 0) {return} - - - - - # Compute the criterion function if requested - - if(icrit==0) { return} - - else { # Compute Linearized CV scores - - call shessd(bvcgr,sg,lg,ng,nbvcg, - nk,smden,qwts,ngpt,hs0,hs1,hs2,hs3) - - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) - if(ierr.ne.0) { return } - - call scvdn(x,w,n,crit,knot,nk, - grid,qwts,ngpt,bvcd,nbvcd,sd,nd,ld, - coef,smden,cvden,lev,abd,p1ip,p2ip,ld4,ldnk) - - - } - - - -return -end //GO.SYSIN DD bart/sslvdn.r echo bart/sslvhz.r 1>&2 sed >bart/sslvhz.r <<'//GO.SYSIN DD bart/sslvhz.r' 's/^-//' -subroutine sslvhz(x,w,c,rn,n,knot,nk, - norm,coef0,coef1,hcoef, - grid,qwts,ngpt, - bvcgr,nbvcg,sg,ng,lg, - bvcd,nbvcd,sd,nd,ld, - coef,smhaz,cvlik,lev,wts, - crit,icrit,spar, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3,ierr) - - -integer n,nk,ngpt,nbvcg,nbvcd, - sd(nk),ld(nk),nd(nk), - sg(nk),lg(nk),ng(nk), - icrit,ld4,ldnk,ierr, - - i,iter,j,max - -real x(n),w(n),c(n),rn(ngpt), - knot(nk+4), - norm,coef0(nk),coef1(nk),hcoef(nk), - grid(ngpt),qwts(ngpt), - bvcgr(nbvcg),bvcd(nbvcd), - coef(nk),smhaz(ngpt),cvlik(n), - lev(ngpt),wts(ngpt), - crit(2),spar, - 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), - wrk1(nk),wrk2(nk),wrk3(nk), - - le,pl,plo,pls,t1,t2,lamr,dof,d0,d1,cmin,cmax, - lambda,sqrt, - sdot,tol,cchk - - - - - max = 15; tol = .1e-3 ; cchk = tol+10. ; ierr = 0 - - - - # Purpose : Solves the smoothing problem and computes - # A linearized CV score - - - # Set Initial Estimate and lambda - - # Initial Estimate - do j=1,nk { hcoef(j) = (1.-spar)*coef0(j)+spar*coef1(j) } - cmin=hcoef(1) ; cmax=hcoef(1) - do j=2,nk { if(cmaxhcoef(j)) { cmin=hcoef(j) } } - - call evlkhz(rn,c,n,hcoef,nk,smhaz,knot,grid,qwts,ngpt, - bvcd,sd,nd,ld,nbvcd,le) - do i=1,ngpt { wts(i) = sqrt(qwts(i)*smhaz(i)*rn(i)) } - call shessh(bvcgr,sg,lg,ng,nbvcg,nk, - smhaz,rn,qwts,ngpt,hs0,hs1,hs2,hs3) - t1 = 0. ; t2 = 0. - do i=2,nk-1 {t1 = t1 + hs0(i) ; t2 = t2 + 2.*sg0(i) } - lamr = t1/t2 - - call dfr(dof,1.,lamr,grid,wts,ngpt,knot,nk,lev, - hs0,hs1,hs2,hs3,sg0,sg1,sg2,sg3,abd,p1ip,ld4,ierr) - - d0 = .75 ; if(nk>60) {d1 = nk/5. -2. } - else {d1 = nk -2. } - d0 = 1./d0 ; d1 = 1./d1 - lambda = (dof-2.)*(d0*d1)/(spar*d1+(1.-spar)*d0) - lambda = (lambda)**4 * lamr - - do j=1,nk { coef(j) = hcoef(j) } - call cbnm(sg0,sg1,sg2,sg3,coef,(nk-3),wrk2) - - - # Starting Value for Penalized Likelihood - - pls = le + lambda*sdot(nk,coef,1,wrk2,1) - - plo = pls - - - - - -# Use Newton's Method with Levenberg-Marquard Line Search to compute estimate -# -# coef = coef0 - [ H + 2 lambda OMEGA ] (-1) dir -# -# coef0 -> wrk1 ; dir -> wrk2 -# - - iter = 0 - - - while(iter tol) { - - iter = iter + 1 - do j=1,nk { wrk1(j) = coef(j) } - - - # Find Newton direction for new estimate - - - call shessh(bvcgr,sg,lg,ng,nbvcg, - nk,smhaz,rn,qwts,ngpt,hs0,hs1,hs2,hs3) - call sdirhz(c,n,bvcgr,sg,lg,ng,nbvcg,nk,qwts,ngpt, - bvcd,sd,ld,nd,nbvcd, - smhaz,rn, - lambda,sg0,sg1,sg2,sg3, - wrk1,wrk2) - - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - - call spbfa(abd,ld4,nk,3,ierr) - - if(ierr.ne.0) { return } - call spbsl(abd,ld4,nk,3,wrk2) - do j=1,nk { wrk2(j) = 0. - wrk2(j) } - - # Find new estimate in the direction of coef - - call srchhz(c,rn,n,cmin,cmax,norm,coef,wrk1, - wrk2,nk,sg0,sg1,sg2,sg3, - smhaz,knot, - grid,qwts,ngpt, - bvcd,sd,nd,ld,nbvcd, - lambda,pl,plo,wrk3) - - -#write(6,*) pl , plo - - if(pl>plo) {plo = pl} - cchk = abs(pl-plo)/abs(pl) ; plo = pl - - - - } - - - - if(pls<=pl) {return} - - - - - # Compute the criterion function if requested - - if(icrit==0) { return} - - else { # Compute Linearized CV score - - call shessh(bvcgr,sg,lg,ng,nbvcg, - nk,smhaz,rn,qwts,ngpt,hs0,hs1,hs2,hs3) - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) - - if(ierr > 0) {return} - - call scvhz(x,w,c,rn,n,crit,icrit,knot,nk, - grid,qwts,ngpt, - bvcd,nbvcd,sd,nd,ld, - bvcgr,nbvcg,sg,ng,lg, - coef,smhaz,cvlik,lev, - abd,p1ip,p2ip,ld4,ldnk, - wrk1,wrk2,wrk3) - } - - - - - - -return -end //GO.SYSIN DD bart/sslvhz.r echo bart/sslvlp.r 1>&2 sed >bart/sslvlp.r <<'//GO.SYSIN DD bart/sslvlp.r' 's/^-//' -subroutine sslvlp(f,z,u,lev,n, - knot,clev,ns, - sbd,ld, - spar,crit,icrit,info, - we,wf,wx) - -real f(n+2),z(n+2),u(n+2),lev(n+2), - knot(ns+6),clev(ns+2), - spar,crit, - we(n+2),wf(n+2),wx(n+2), - sbd(ld,n+2), - del, - rss,df - - -integer n,ns,icrit,ld,info,i - - - # Purpose : Solves the smoothing problem and computes the - # criterion function (OCV or GCV). - - - # Setup sbd and u - - del = 10.**(2.5*spar) - .995 - - for(i=1;i<=n+2;i=i+1){ - sbd(3,i) = 6.*del + f(i) ; u(i) = f(i)*z(i) } - sbd(3,1) = sbd(3,1)+5.*del ; sbd(3,n+2) = sbd(3,n+2)+5.*del - sbd(3,2) = sbd(3,2)+del ; sbd(3,n+1) = sbd(3,n+1)+del - - do i=2,n { sbd(2,i+1) = -4.*del } - sbd(2,2) = -6.*del ; sbd(2,n+2) = -6.*del - - do i=1,n { sbd(1,i+2) = del } - - # Solve system putting the result in u - - call spbfa(sbd,ld,n+2,2,info) - if(info.ne.0) {return} - call spbsl(sbd,ld,n+2,2,u) - - - - # Compute the criterion function if requested - - if(icrit==0) { return} - - else { # Ordinary or Generalized CV - - # Compute Approximate Leverages - - call lsamp1(f,n+2,sbd,ld,ns,lev,wx,we,wf,knot,clev) - - # Evaluate Criterion - - - if(icrit==1) { # Generalized CV - - rss = 0. ; df = 0. - do i=1,n { rss = rss + ((z(i+1)-u(i+1))**2)*f(i+1)} - do i=1,n { df = df + (1.-lev(i+1))} - crit = (rss/n)/((df/n)**2) } - - else { # Ordinary CV - - crit = 0. - do i=1,n { crit = crit + - f(i+1)*((z(i+1)-u(i+1))/(1.-lev(i+1)))**2 }} - - - return } - -end //GO.SYSIN DD bart/sslvlp.r echo bart/sslvrg.r 1>&2 sed >bart/sslvrg.r <<'//GO.SYSIN DD bart/sslvrg.r' 's/^-//' -subroutine sslvrg(x,y,w,n,knot,nk, - coef,sz,lev, - crit,icrit, - spar,lamr,dof, - xwy, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - abd,p1ip,p2ip,ld4,ldnk,info) - - -integer n,nk,icrit,ld4,ldnk,i,icoef,ileft,ilo,info,j,mflag - -real x(n),y(n),w(n), - knot(nk+4), - coef(nk),sz(n),lev(n), - crit, - spar,lamr,dof, - 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(15), - xv,bvalu,rss,df,d0,d1 - - - - ilo = 1 ; eps = .1e-10 - - # Purpose : Solves the smoothing problem and computes the - # criterion function (OCV or GCV). - - # Set lambda - d0 = .75 ; if(nk>60) {d1 = nk/5. -2. } - else {d1 = nk -2. } - d0 = 1./d0 ; d1 = 1./d1 - lambda = (dof-2.)*(d0*d1)/(spar*d1+(1.-spar)*d0) - - lambda = (lambda)**4 * lamr - - - - # The coeficients of estimated smooth - - - do i=1,nk { coef(i) = xwy(i) } - do i=1,nk { abd(4,i) = hs0(i)+lambda*sg0(i) } - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lambda*sg3(i) } - - - call spbfa(abd,ld4,nk,3,info) - if(info.ne.0) { return } - call spbsl(abd,ld4,nk,3,coef) - - - # Value of smooth at the data points - - icoef = 1 - do i=1,n { xv = x(i) - sz(i) = bvalu(knot,coef, - nk,4,0,xv,icoef,work(1)) } - - - # Compute the criterion function if requested - - if(icrit==0) { return} - - else { # Ordinary or Generalized CV - - # Get Leverages First - - call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) - - do i=1,n { xv = x(i) - call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) - if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } - if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } - j=ileft-3 - call bspvd(knot,4,1,xv,ileft,4,vnikx,work) - 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 } - - - - - - # Evaluate Criterion - - if(icrit==1) { # Generalized CV - - rss = 0. ; df = 0. - do i=1,n { rss = rss + ((y(i)-sz(i))*w(i))**2} - do i=1,n { df = df + 1.-lev(i)} - if((n-df)<2.) { - do i=1,n { lev(i) = (2./(n-df))*lev(i)} - df = n-2. } - - crit = (rss/n)/((df/n)**2) } - - else { # Ordinary CV - - crit = 0. - do i=1,n { crit = crit + - (((y(i)-sz(i))*w(i))/(1-lev(i)))**2 }} - - - return } - -end //GO.SYSIN DD bart/sslvrg.r echo bart/step.r 1>&2 sed >bart/step.r <<'//GO.SYSIN DD bart/step.r' 's/^-//' -subroutine step(alpha,tau,nco,diag,cmin,cmax,coef0,coef,dir,nk, - x,n,knot, - abd,ld4, - hs0,hs1,hs2,hs3, - sg0,sg1,sg2,sg3, - lambda) - - -real alpha,tau,nco,diag, - cmin,cmax,coef0(nk),coef(nk),dir(nk), - x(n),knot(nk+4), - abd(ld4,nk), - hs0(nk),hs1(nk),hs2(nk),hs3(nk), - sg0(nk),sg1(nk),sg2(nk),sg3(nk), - lambda - -integer nk,n,ld4,ierr,j - - - # Local Variables -real bvalu,corr,snrm2,sqrt,work(15) -integer i,icoef - - - # - # Evaluate the Levenberg-Marquart estimate for the given alpha value - # - - - - call scopy(nk,dir,1,coef,1) - - do i=1,nk { abd(4,i) = hs0(i)+2.*lambda*sg0(i)+diag*alpha} - do i=1,(nk-1) { abd(3,i+1) = hs1(i)+2.*lambda*sg1(i) } - do i=1,(nk-2) { abd(2,i+2) = hs2(i)+2.*lambda*sg2(i) } - do i=1,(nk-3) { abd(1,i+3) = hs3(i)+2.*lambda*sg3(i) } - - call spbfa(abd,ld4,nk,3,ierr) - call spbsl(abd,ld4,nk,3,coef) - - - # Adjust coef so that the constraints aren't violated - - do j=1,nk { coef(j) = coef0(j) - coef(j) } - - do j=1,nk { - if(coef(j)cmax+1.9) { coef(j) = cmax + 1.9 }} - - corr = 0. ; icoef = 1 - do i=1,n{ corr = corr + bvalu(knot,coef,nk,4,0,x(i),icoef,work) } - corr = corr/n - - do j=1,nk { coef(j) = coef(j) - corr - if(coef(j)cmax+1.9) { coef(j) = cmax + 1.9 }} - - corr = 0. - do i=1,n{ corr = corr + bvalu(knot,coef,nk,4,0,x(i),icoef,work) } - corr = corr/n - do j=1,nk { coef(j) = coef0(j) - coef(j) + corr } - - - - - # Compute Step Size - - tau = snrm2(nk,coef,1)/(sqrt(nk*1.)*(nco+1.)) - - - -return -end //GO.SYSIN DD bart/step.r echo bart/stepp.r 1>&2 sed >bart/stepp.r <<'//GO.SYSIN DD bart/stepp.r' 's/^-//' -subroutine stepp(alpha,tau,nco,diag,coef,dir,p, - svdhp,d,wd,ldp,wrk0) - - - -real alpha,tau,nco,diag, - coef(p),dir(p), - svdhp(ldp,p),d(p),wd(p),wrk0(p) - -integer p,ldp - - - # Local Variables -real snrm2,sqrt,sdot -integer j - - - # - # Evaluate the Levenberg-Marquart estimate for the given alpha value - # - - - - do j=1,p { wd(j) = d(j) + diag*alpha } - do j=1,p { coef(j) = sdot(p,dir,1,svdhp(1,j),1)/wd(j) } - call scopy(p,coef,1,wrk0,1) - do j=1,p { coef(j) = sdot(p,wrk0,1,svdhp(j,1),ldp) } - - - - # Compute Step Size - - tau = snrm2(p,coef,1)/(sqrt(p*1.)*(nco+1.)) - - - -return -end //GO.SYSIN DD bart/stepp.r echo bart/strap.r 1>&2 sed >bart/strap.r <<'//GO.SYSIN DD bart/strap.r' 's/^-//' - - # Approximates integrals by the Trapezoidal rule - - - - real function intg1(f,x,n) # Single function - - real f(n),x(n) - integer n,i - - intg1 = .5*f(1)*(x(2)-x(1)) - - do i=2,n-1 {intg1 = intg1 + .5*f(i)*(x(i+1)-x(i-1)) } - - intg1 = intg1 + .5*f(n)*(x(n)-x(n-1)) -return -end - - real function intg2(f,g,x,n) # 2 functions - - real f(n),g(n),x(n) - integer n,i - - intg2 = .5*f(1)*g(1)*(x(2)-x(1)) - - do i=2,n-1 {intg2 = intg2 + .5*f(i)*g(i)*(x(i+1)-x(i-1)) } - - intg2 = intg2 + .5*f(n)*g(n)*(x(n)-x(n-1)) -return -end - - real function intg3(f,g,h,x,n) # 3 functions - - real f(n),g(n),h(n),x(n) - integer n,i - - intg3 = .5*f(1)*g(1)*h(1)*(x(2)-x(1)) - - do i=2,n-1 {intg3 = intg3 + .5*f(i)*g(i)*h(i)*(x(i+1)-x(i-1)) } - - intg3 = intg3 + .5*f(n)*g(n)*h(n)*(x(n)-x(n-1)) -return -end - - real function intg4(f,g,h,k,x,n) # 4 functions - - real f(n),g(n),h(n),k(n),x(n) - integer n,i - - intg4 = .5*f(1)*g(1)*h(1)*k(1)*(x(2)-x(1)) - - do i=2,n-1 {intg4 = intg4 + .5*f(i)*g(i)*h(i)*k(i)*(x(i+1)-x(i-1)) } - - intg4 = intg4 + .5*f(n)*g(n)*h(n)*k(n)*(x(n)-x(n-1)) -return -end - - # Approximates integrals by the Gaussian rule - -# x are Gaussian Quadrature Weights - - real function gq1(f,x,n) # Single function - - integer n - - real f(n),x(n) - - - gq1 = 0. - do i=1,n { gq1 = gq1 + f(i)*x(i) } - - - -return -end - - real function gq2(f,g,x,n) # 2 functions - - integer n,i - real f(n),g(n),x(n) - - gq2 = 0. - do i=1,n { gq2 = gq2 + f(i)*g(i)*x(i) } - - -return -end - - real function gq3(f,g,h,x,n) # 3 functions - - integer n,i - real f(n),g(n),h(n),x(n) - - gq3 = 0. - do i=1,n { gq3 = gq3 + f(i)*g(i)*h(i)*x(i) } - - -return -end - - real function gq4(f,g,h,k,x,n) # 4 functions - - integer n,i - real f(n),g(n),h(n),k(n),x(n) - - gq4 = 0. - do i=1,n { gq4 = gq4 + f(i)*g(i)*h(i)*k(i)*x(i) } - -return -end //GO.SYSIN DD bart/strap.r echo bart/stxwx.r 1>&2 sed >bart/stxwx.r <<'//GO.SYSIN DD bart/stxwx.r' 's/^-//' -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(15) # local - -integer k,n, - - j,i,ilo,ileft,mflag # local - - - - - # Initialise the output vectors - - do i=1,n { y(i)=0. ; hs0(i)=0. ; hs1(i)=0. - hs2(i)=0. ; hs3(i)=0. } - - - # Compute X'WX -> hs0,hs1,hs2,hs3 and X'WZ -> y - - ilo=1 ; eps = .1e-9 - - do i=1,k { - - call intrv(xknot(1),(n+1),x(i),ilo,ileft,mflag) - if(mflag==-1) {write(6,'("Error in hess ",i2)')mflag;stop} - if(mflag== 1) {if(x(i)<=(xknot(ileft)+eps)){ileft=ileft-1} - else{write(6,'("Error in hess ",i2)')mflag;stop}} - - - call bspvd (xknot,4,1,x(i),ileft,4,vnikx,work) - - 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 } - -return -end //GO.SYSIN DD bart/stxwx.r