| 1 | # This computes the univariate smoothing spline, automatically choosing |
|---|
| 2 | # the smoothing parameter by minimizing generalized cross validation. It |
|---|
| 3 | # was written by Finbarr O'Sullivan using the scheme described in |
|---|
| 4 | # "Comments on Dr. Silverman's Paper", J. Royal Statistical Society B |
|---|
| 5 | # (1985) 47, pp.39-40. |
|---|
| 6 | # This version was installed in netlib 14 Mar 1987. |
|---|
| 7 | # Remember to check whether your program needs to square or sqrt the |
|---|
| 8 | # weights before calling this (or any other) least squares routine. |
|---|
| 9 | # The code has been slightly modified by Trevor Hastie and Eric Grosse |
|---|
| 10 | # to use de Boor's spline routines from netlib/pppack. Also, auxiliary |
|---|
| 11 | # routines "setreg" (for sorting and standardizing the range of the data) |
|---|
| 12 | # and "sknotl" (for filling the knot array) have been appended. |
|---|
| 13 | |
|---|
| 14 | |
|---|
| 15 | subroutine sbart(xs,ys,ws,n,knot,nk, |
|---|
| 16 | coef,sz,lev, |
|---|
| 17 | crit,icrit,spar,ispar,lspar,uspar,tol, |
|---|
| 18 | isetup, |
|---|
| 19 | xwy, |
|---|
| 20 | hs0,hs1,hs2,hs3, |
|---|
| 21 | sg0,sg1,sg2,sg3, |
|---|
| 22 | abd,p1ip,p2ip,ld4,ldnk,ier) |
|---|
| 23 | |
|---|
| 24 | |
|---|
| 25 | # A Cubic B-spline Smoothing routine. |
|---|
| 26 | |
|---|
| 27 | # |
|---|
| 28 | # The algorithm minimises: |
|---|
| 29 | # |
|---|
| 30 | # (1/n) * sum ws(i)**2 * (ys(i)-sz(i))**2 + lambda* int ( sz"(xs) )**2 dxs |
|---|
| 31 | # |
|---|
| 32 | # lambda is a function of the spar which is assumed to be between |
|---|
| 33 | # 0 and 1 |
|---|
| 34 | |
|---|
| 35 | |
|---|
| 36 | # Input |
|---|
| 37 | |
|---|
| 38 | # n number of data points |
|---|
| 39 | # ys(n) vector of length n containing the observations |
|---|
| 40 | # ws(n) vector containing the weights given to each data point |
|---|
| 41 | # xs(n) vector containing the ordinates of the observations |
|---|
| 42 | |
|---|
| 43 | |
|---|
| 44 | # nk number of b-spline coefficients to be estimated |
|---|
| 45 | # nk <= n+2 |
|---|
| 46 | # knot(nk+4) vector of knot points defining the cubic b-spline basis. |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | # spar penalised likelihood smoothing parameter |
|---|
| 50 | # ispar indicator saying if spar is supplied or to be estimated |
|---|
| 51 | # lspar, uspar lower and upper values for spar 0.,1. are good values |
|---|
| 52 | # tol used in Golden Search routine |
|---|
| 53 | |
|---|
| 54 | # isetup setup indicator |
|---|
| 55 | |
|---|
| 56 | # icrit indicator saying which cross validation score |
|---|
| 57 | # is to be computed |
|---|
| 58 | |
|---|
| 59 | # ld4 the leading dimension of abd (ie ld4=4) |
|---|
| 60 | # ldnk the leading dimension of p2ip (not referenced) |
|---|
| 61 | |
|---|
| 62 | |
|---|
| 63 | # Output |
|---|
| 64 | |
|---|
| 65 | # coef(nk) vector of spline coefficients |
|---|
| 66 | # sz(n) vector of smoothed z-values |
|---|
| 67 | # lev(n) vector of leverages |
|---|
| 68 | # crit either ordinary of generalized CV score |
|---|
| 69 | # ier error indicator |
|---|
| 70 | # ier = 0 ___ everything fine |
|---|
| 71 | # ier = 1 ___ spar too small or too big |
|---|
| 72 | # problem in cholesky decomposition |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | |
|---|
| 76 | # Working arrays/matrix |
|---|
| 77 | # xwy X'Wy |
|---|
| 78 | # hs0,hs1,hs2,hs3 the diagonals of the X'WX matrix |
|---|
| 79 | # sg0,sg1,sg2,sg3 the diagonals of the Gram matrix |
|---|
| 80 | # abd(ld4,nk) [ X'WX+lambda*SIGMA] in diagonal form |
|---|
| 81 | # p1ip(ld4,nk) inner products between columns of L inverse |
|---|
| 82 | # p2ip(ldnk,nk) all inner products between columns of L inverse |
|---|
| 83 | # L'L = [X'WX+lambdaSIGMA] NOT REFERENCED |
|---|
| 84 | |
|---|
| 85 | real xs(n),ys(n),ws(n), |
|---|
| 86 | knot(nk+4), |
|---|
| 87 | coef(nk),sz(n),lev(n), |
|---|
| 88 | crit,spar,lspar,uspar,tol, |
|---|
| 89 | xwy(nk), |
|---|
| 90 | hs0(nk),hs1(nk),hs2(nk),hs3(nk), |
|---|
| 91 | sg0(nk),sg1(nk),sg2(nk),sg3(nk), |
|---|
| 92 | abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk) |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | integer n,nk,isetup,icrit,ispar,ld4,ldnk,ier |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | |
|---|
| 100 | # Local variables |
|---|
| 101 | |
|---|
| 102 | real t1,t2,ratio, |
|---|
| 103 | a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w, |
|---|
| 104 | fu,fv,fw,fx,x, |
|---|
| 105 | ax,bx |
|---|
| 106 | |
|---|
| 107 | integer i |
|---|
| 108 | |
|---|
| 109 | |
|---|
| 110 | |
|---|
| 111 | |
|---|
| 112 | |
|---|
| 113 | |
|---|
| 114 | |
|---|
| 115 | |
|---|
| 116 | # Compute SIGMA, X' W**2 X, X' W**2 z, trace ratio, s0, s1. |
|---|
| 117 | |
|---|
| 118 | # SIGMA -> sg0,sg1,sg2,sg3 |
|---|
| 119 | # X' W**2 X -> hs0,hs1,hs2,hs3 |
|---|
| 120 | # X' W**2 Z -> xwy |
|---|
| 121 | |
|---|
| 122 | if(isetup==0){ |
|---|
| 123 | call sgram(sg0,sg1,sg2,sg3,knot,nk) |
|---|
| 124 | call stxwx(xs,ys,ws,n,knot,nk,xwy,hs0,hs1,hs2,hs3) |
|---|
| 125 | # check |
|---|
| 126 | # print 999;999 format(" got through check ") |
|---|
| 127 | |
|---|
| 128 | |
|---|
| 129 | t1=0. ; t2=0. |
|---|
| 130 | do i=3,nk-3 { t1 = t1 + hs0(i) } |
|---|
| 131 | do i=3,nk-3 { t2 = t2 + sg0(i) } |
|---|
| 132 | ratio = t1/t2 |
|---|
| 133 | |
|---|
| 134 | isetup = 1 } |
|---|
| 135 | |
|---|
| 136 | # check 1 |
|---|
| 137 | # print 1999;1999 format(" got through check 1") |
|---|
| 138 | |
|---|
| 139 | |
|---|
| 140 | # Compute estimate |
|---|
| 141 | |
|---|
| 142 | |
|---|
| 143 | if(ispar==1) { # Value of spar supplied |
|---|
| 144 | |
|---|
| 145 | call sslvrg(xs,ys,ws,n,knot,nk, |
|---|
| 146 | coef,sz,lev,crit,icrit, |
|---|
| 147 | spar,ratio, |
|---|
| 148 | xwy, |
|---|
| 149 | hs0,hs1,hs2,hs3, |
|---|
| 150 | sg0,sg1,sg2,sg3, |
|---|
| 151 | abd,p1ip,p2ip,ld4,ldnk,ier) |
|---|
| 152 | # check 2 |
|---|
| 153 | # print 2999;2999 format(" got through check 2") |
|---|
| 154 | |
|---|
| 155 | return } |
|---|
| 156 | |
|---|
| 157 | else { |
|---|
| 158 | |
|---|
| 159 | # Use Forsythe Malcom and Moler routine to minimise criterion |
|---|
| 160 | |
|---|
| 161 | ax=lspar ; bx=uspar # f denotes the value of the criterion |
|---|
| 162 | |
|---|
| 163 | # |
|---|
| 164 | # an approximation x to the point where f attains a minimum on |
|---|
| 165 | # the interval (ax,bx) is determined. |
|---|
| 166 | # |
|---|
| 167 | # |
|---|
| 168 | # input.. |
|---|
| 169 | # |
|---|
| 170 | # ax left endpoint of initial interval |
|---|
| 171 | # bx right endpoint of initial interval |
|---|
| 172 | # f function subprogram which evaluates f(x) for any x |
|---|
| 173 | # in the interval (ax,bx) |
|---|
| 174 | # tol desired length of the interval of uncertainty of the final |
|---|
| 175 | # result ( .ge. 0.0) |
|---|
| 176 | # |
|---|
| 177 | # |
|---|
| 178 | # output.. |
|---|
| 179 | # |
|---|
| 180 | # fmin abcissa approximating the point where f attains a minimum |
|---|
| 181 | # |
|---|
| 182 | # |
|---|
| 183 | # the method used is a combination of golden section search and |
|---|
| 184 | # successive parabolic interpolation. convergence is never much slower |
|---|
| 185 | # than that for a fibonacci search. if f has a continuous second |
|---|
| 186 | # derivative which is positive at the minimum (which is not at ax or |
|---|
| 187 | # bx), then convergence is superlinear, and usually of the order of |
|---|
| 188 | # about 1.324.... |
|---|
| 189 | # the function f is never evaluated at two points closer together |
|---|
| 190 | # than eps*abs(fmin) + (tol/3), where eps is approximately the square |
|---|
| 191 | # root of the relative machine precision. if f is a unimodal |
|---|
| 192 | # function and the computed values of f are always unimodal when |
|---|
| 193 | # separated by at least eps*abs(x) + (tol/3), then fmin approximates |
|---|
| 194 | # the abcissa of the global minimum of f on the interval ax,bx with |
|---|
| 195 | # an error less than 3*eps*abs(fmin) + tol. if f is not unimodal, |
|---|
| 196 | # then fmin may approximate a local, but perhaps non-global, minimum to |
|---|
| 197 | # the same accuracy. |
|---|
| 198 | # this function subprogram is a slightly modified version of the |
|---|
| 199 | # algol 60 procedure localmin given in richard brent, algorithms for |
|---|
| 200 | # minimization without derivatives, prentice - hall, inc. (1973). |
|---|
| 201 | # |
|---|
| 202 | # |
|---|
| 203 | # real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w |
|---|
| 204 | # real fu,fv,fw,fx,x |
|---|
| 205 | # |
|---|
| 206 | # c is the squared inverse of the golden ratio |
|---|
| 207 | # |
|---|
| 208 | c = 0.5*(3. - sqrt(5.0)) |
|---|
| 209 | # |
|---|
| 210 | # eps is approximately the square root of the relative machine |
|---|
| 211 | # precision. |
|---|
| 212 | # |
|---|
| 213 | eps = 1.00 |
|---|
| 214 | 10 eps = eps/2.00 |
|---|
| 215 | tol1 = 1.0 + eps |
|---|
| 216 | if (tol1 .gt. 1.00) go to 10 |
|---|
| 217 | eps = sqrt(eps) |
|---|
| 218 | # |
|---|
| 219 | # initialization |
|---|
| 220 | # |
|---|
| 221 | a = ax |
|---|
| 222 | b = bx |
|---|
| 223 | v = a + c*(b - a) |
|---|
| 224 | w = v |
|---|
| 225 | x = v |
|---|
| 226 | e = 0.0 |
|---|
| 227 | |
|---|
| 228 | spar = x |
|---|
| 229 | call sslvrg(xs,ys,ws,n,knot,nk, |
|---|
| 230 | coef,sz,lev,crit,icrit, |
|---|
| 231 | spar,ratio, |
|---|
| 232 | xwy, |
|---|
| 233 | hs0,hs1,hs2,hs3, |
|---|
| 234 | sg0,sg1,sg2,sg3, |
|---|
| 235 | abd,p1ip,p2ip,ld4,ldnk,ier) |
|---|
| 236 | |
|---|
| 237 | fx = crit |
|---|
| 238 | fv = fx |
|---|
| 239 | fw = fx |
|---|
| 240 | # |
|---|
| 241 | # main loop starts here |
|---|
| 242 | # |
|---|
| 243 | 20 xm = 0.5*(a + b) |
|---|
| 244 | tol1 = eps*abs(x) + tol/3.0 |
|---|
| 245 | tol2 = 2.0*tol1 |
|---|
| 246 | # |
|---|
| 247 | # check stopping criterion |
|---|
| 248 | # |
|---|
| 249 | if (abs(x - xm) .le. (tol2 - 0.5*(b - a))) go to 90 |
|---|
| 250 | # |
|---|
| 251 | # is golden-section necessary |
|---|
| 252 | # |
|---|
| 253 | if (abs(e) .le. tol1) go to 40 |
|---|
| 254 | # |
|---|
| 255 | # fit parabola |
|---|
| 256 | # |
|---|
| 257 | r = (x - w)*(fx - fv) |
|---|
| 258 | q = (x - v)*(fx - fw) |
|---|
| 259 | p = (x - v)*q - (x - w)*r |
|---|
| 260 | q = 2.00*(q - r) |
|---|
| 261 | if (q .gt. 0.0) p = -p |
|---|
| 262 | q = abs(q) |
|---|
| 263 | r = e |
|---|
| 264 | e = d |
|---|
| 265 | # |
|---|
| 266 | # is parabola acceptable |
|---|
| 267 | # |
|---|
| 268 | 30 if (abs(p) .ge. abs(0.5*q*r)) go to 40 |
|---|
| 269 | if (p .le. q*(a - x)) go to 40 |
|---|
| 270 | if (p .ge. q*(b - x)) go to 40 |
|---|
| 271 | # |
|---|
| 272 | # a parabolic interpolation step |
|---|
| 273 | # |
|---|
| 274 | d = p/q |
|---|
| 275 | u = x + d |
|---|
| 276 | # |
|---|
| 277 | # f must not be evaluated too close to ax or bx |
|---|
| 278 | # |
|---|
| 279 | if ((u - a) .lt. tol2) d = sign(tol1, xm - x) |
|---|
| 280 | if ((b - u) .lt. tol2) d = sign(tol1, xm - x) |
|---|
| 281 | go to 50 |
|---|
| 282 | # |
|---|
| 283 | # a golden-section step |
|---|
| 284 | # |
|---|
| 285 | 40 if (x .ge. xm) e = a - x |
|---|
| 286 | if (x .lt. xm) e = b - x |
|---|
| 287 | d = c*e |
|---|
| 288 | # |
|---|
| 289 | # f must not be evaluated too close to x |
|---|
| 290 | # |
|---|
| 291 | 50 if (abs(d) .ge. tol1) u = x + d |
|---|
| 292 | if (abs(d) .lt. tol1) u = x + sign(tol1, d) |
|---|
| 293 | |
|---|
| 294 | spar = u |
|---|
| 295 | call sslvrg(xs,ys,ws,n,knot,nk, |
|---|
| 296 | coef,sz,lev,crit,icrit, |
|---|
| 297 | spar,ratio, |
|---|
| 298 | xwy, |
|---|
| 299 | hs0,hs1,hs2,hs3, |
|---|
| 300 | sg0,sg1,sg2,sg3, |
|---|
| 301 | abd,p1ip,p2ip,ld4,ldnk,ier) |
|---|
| 302 | |
|---|
| 303 | fu = crit |
|---|
| 304 | # |
|---|
| 305 | # update a, b, v, w, and x |
|---|
| 306 | # |
|---|
| 307 | if (fu .gt. fx) go to 60 |
|---|
| 308 | if (u .ge. x) a = x |
|---|
| 309 | if (u .lt. x) b = x |
|---|
| 310 | v = w |
|---|
| 311 | fv = fw |
|---|
| 312 | w = x |
|---|
| 313 | fw = fx |
|---|
| 314 | x = u |
|---|
| 315 | fx = fu |
|---|
| 316 | go to 20 |
|---|
| 317 | 60 if (u .lt. x) a = u |
|---|
| 318 | if (u .ge. x) b = u |
|---|
| 319 | if (fu .le. fw) go to 70 |
|---|
| 320 | if (w .eq. x) go to 70 |
|---|
| 321 | if (fu .le. fv) go to 80 |
|---|
| 322 | if (v .eq. x) go to 80 |
|---|
| 323 | if (v .eq. w) go to 80 |
|---|
| 324 | go to 20 |
|---|
| 325 | 70 v = w |
|---|
| 326 | fv = fw |
|---|
| 327 | w = u |
|---|
| 328 | fw = fu |
|---|
| 329 | go to 20 |
|---|
| 330 | 80 v = u |
|---|
| 331 | fv = fu |
|---|
| 332 | go to 20 |
|---|
| 333 | # |
|---|
| 334 | # end of main loop |
|---|
| 335 | # |
|---|
| 336 | 90 continue ; spar = x ; crit = fx |
|---|
| 337 | return } |
|---|
| 338 | |
|---|
| 339 | |
|---|
| 340 | |
|---|
| 341 | |
|---|
| 342 | |
|---|
| 343 | return |
|---|
| 344 | end |
|---|
| 345 | subroutine sgram(sg0,sg1,sg2,sg3,tb,nb) |
|---|
| 346 | |
|---|
| 347 | real sg0(nb),sg1(nb),sg2(nb),sg3(nb),tb(nb+4), |
|---|
| 348 | vnikx(4,3),work(16),yw1(4),yw2(4), |
|---|
| 349 | wpt |
|---|
| 350 | |
|---|
| 351 | integer nb,ileft,ilo,mflag, |
|---|
| 352 | i,ii,jj # indices |
|---|
| 353 | |
|---|
| 354 | |
|---|
| 355 | #PURPOSE |
|---|
| 356 | |
|---|
| 357 | # Calculation of the cubic B-spline smoothness prior |
|---|
| 358 | # for "usual" interior knot setup. |
|---|
| 359 | |
|---|
| 360 | |
|---|
| 361 | # Uses BSPVD and INTRV in the CMLIB |
|---|
| 362 | |
|---|
| 363 | |
|---|
| 364 | |
|---|
| 365 | |
|---|
| 366 | # sgm[0-3](nb) Symmetric matrix |
|---|
| 367 | # whose (i,j)'th element contains the integral of |
|---|
| 368 | # B''(i,.) B''(j,.) , i=1,2 ... nb and j=i,...nb. |
|---|
| 369 | # Only the upper four diagonals are computed. |
|---|
| 370 | |
|---|
| 371 | |
|---|
| 372 | #Initialise the sigma vectors |
|---|
| 373 | do i=1,nb{ sg0(i)=0.;sg1(i)=0.;sg2(i)=0.;sg3(i)=0.} |
|---|
| 374 | |
|---|
| 375 | ilo = 1 |
|---|
| 376 | |
|---|
| 377 | do i=1,nb { |
|---|
| 378 | |
|---|
| 379 | # Calculate a linear approximation to the |
|---|
| 380 | # second derivative of the non-zero B-splines |
|---|
| 381 | # over the interval [tb(i),tb(i+1)]. |
|---|
| 382 | |
|---|
| 383 | #call intrv(tb(1),(nb+1),tb(i),ilo,ileft,mflag) |
|---|
| 384 | call interv(tb(1),(nb+1),tb(i),ileft,mflag) |
|---|
| 385 | |
|---|
| 386 | |
|---|
| 387 | #Left end second derivatives |
|---|
| 388 | #call bspvd (tb,4,3,tb(i),ileft,4,vnikx,work) |
|---|
| 389 | call bsplvd (tb,4,tb(i),ileft,work,vnikx,3) |
|---|
| 390 | |
|---|
| 391 | # Put values into yw1 |
|---|
| 392 | do ii=1,4 { yw1(ii) = vnikx(ii,3) } |
|---|
| 393 | |
|---|
| 394 | |
|---|
| 395 | #Right end second derivatives |
|---|
| 396 | #call bspvd (tb,4,3,tb(i+1),ileft,4,vnikx,work) |
|---|
| 397 | call bsplvd (tb,4,tb(i+1),ileft,work,vnikx,3) |
|---|
| 398 | |
|---|
| 399 | |
|---|
| 400 | # Slope*(length of interval) in Linear Approximation to B'' |
|---|
| 401 | do ii=1,4 { yw2(ii) = vnikx(ii,3) - yw1(ii) } |
|---|
| 402 | wpt = tb(i+1) - tb(i) |
|---|
| 403 | |
|---|
| 404 | |
|---|
| 405 | # Calculate Contributions to the simga vectors |
|---|
| 406 | |
|---|
| 407 | if(ileft>=4){ |
|---|
| 408 | do ii=1,4 { |
|---|
| 409 | |
|---|
| 410 | jj=ii |
|---|
| 411 | sg0(ileft-4+ii) = sg0(ileft-4+ii) + |
|---|
| 412 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 413 | yw2(ii)*yw2(jj)*.3330 ) |
|---|
| 414 | jj=ii+1 |
|---|
| 415 | if(jj<=4) {sg1(ileft+ii-4) = sg1(ileft+ii-4) + |
|---|
| 416 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 417 | yw2(ii)*yw2(jj)*.3330 )} |
|---|
| 418 | jj=ii+2 |
|---|
| 419 | if(jj<=4) {sg2(ileft+ii-4) = sg2(ileft+ii-4) + |
|---|
| 420 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 421 | yw2(ii)*yw2(jj)*.3330 )} |
|---|
| 422 | jj=ii+3 |
|---|
| 423 | if(jj<=4) {sg3(ileft+ii-4) = sg3(ileft+ii-4) + |
|---|
| 424 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 425 | yw2(ii)*yw2(jj)*.3330 )} |
|---|
| 426 | |
|---|
| 427 | } |
|---|
| 428 | } |
|---|
| 429 | |
|---|
| 430 | else if(ileft==3){ |
|---|
| 431 | do ii=1,3 { |
|---|
| 432 | |
|---|
| 433 | jj=ii |
|---|
| 434 | sg0(ileft-3+ii) = sg0(ileft-3+ii) + |
|---|
| 435 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 436 | yw2(ii)*yw2(jj)*.3330 ) |
|---|
| 437 | jj=ii+1 |
|---|
| 438 | if(jj<=3) {sg1(ileft+ii-3) = sg1(ileft+ii-3) + |
|---|
| 439 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 440 | yw2(ii)*yw2(jj)*.3330 )} |
|---|
| 441 | jj=ii+2 |
|---|
| 442 | if(jj<=3) {sg2(ileft+ii-3) = sg2(ileft+ii-3) + |
|---|
| 443 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 444 | yw2(ii)*yw2(jj)*.3330 )} |
|---|
| 445 | |
|---|
| 446 | } |
|---|
| 447 | } |
|---|
| 448 | |
|---|
| 449 | else if(ileft==2){ |
|---|
| 450 | do ii=1,2 { |
|---|
| 451 | |
|---|
| 452 | jj=ii |
|---|
| 453 | sg0(ileft-2+ii) = sg0(ileft-2+ii) + |
|---|
| 454 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 455 | yw2(ii)*yw2(jj)*.3330 ) |
|---|
| 456 | jj=ii+1 |
|---|
| 457 | if(jj<=2) {sg1(ileft+ii-2) = sg1(ileft+ii-2) + |
|---|
| 458 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 459 | yw2(ii)*yw2(jj)*.3330 )} |
|---|
| 460 | |
|---|
| 461 | } |
|---|
| 462 | } |
|---|
| 463 | |
|---|
| 464 | else if(ileft==1){ |
|---|
| 465 | do ii=1,1 { |
|---|
| 466 | |
|---|
| 467 | jj=ii |
|---|
| 468 | sg0(ileft-1+ii) = sg0(ileft-1+ii) + |
|---|
| 469 | wpt* (yw1(ii)*yw1(jj) + (yw2(ii)*yw1(jj) + yw2(jj)*yw1(ii))*.50 + |
|---|
| 470 | yw2(ii)*yw2(jj)*.3330 ) |
|---|
| 471 | |
|---|
| 472 | } |
|---|
| 473 | }} |
|---|
| 474 | |
|---|
| 475 | |
|---|
| 476 | return |
|---|
| 477 | end |
|---|
| 478 | subroutine sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,flag) |
|---|
| 479 | |
|---|
| 480 | real abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk), |
|---|
| 481 | wjm3(3),wjm2(2),wjm1(1),c0,c1,c2,c3 |
|---|
| 482 | |
|---|
| 483 | integer flag,ld4,nk,ldnk,i,j,k |
|---|
| 484 | |
|---|
| 485 | |
|---|
| 486 | |
|---|
| 487 | |
|---|
| 488 | # Purpose : Computes Inner Products between columns of L(-1) |
|---|
| 489 | # where L = abd is a Banded Matrix with 3 subdiagonals |
|---|
| 490 | |
|---|
| 491 | # The algorithm works in two passes: |
|---|
| 492 | # |
|---|
| 493 | # Pass 1 computes (cj,ck) k=j,j-1,j-2,j-3 ,j=nk, .. 1 |
|---|
| 494 | # Pass 2 computes (cj,ck) k<=j-4 (If flag == 1 ). |
|---|
| 495 | # |
|---|
| 496 | # A refinement of Elden's trick is used. |
|---|
| 497 | # |
|---|
| 498 | |
|---|
| 499 | |
|---|
| 500 | |
|---|
| 501 | |
|---|
| 502 | # Pass 1 |
|---|
| 503 | |
|---|
| 504 | |
|---|
| 505 | wjm3(1)=0. ; wjm3(2)=0. ; wjm3(1)=0. |
|---|
| 506 | wjm2(1)=0. ; wjm2(2)=0. |
|---|
| 507 | wjm1(1)=0. |
|---|
| 508 | |
|---|
| 509 | do i=1,nk { j=nk-i+1 |
|---|
| 510 | |
|---|
| 511 | c0 = 1./abd(4,j) |
|---|
| 512 | if(j<=nk-3) {c1 = abd(1,j+3)*c0 |
|---|
| 513 | c2 = abd(2,j+2)*c0 |
|---|
| 514 | c3 = abd(3,j+1)*c0 } |
|---|
| 515 | else if(j==nk-2) {c1 = 0. |
|---|
| 516 | c2 = abd(2,j+2)*c0 |
|---|
| 517 | c3 = abd(3,j+1)*c0 } |
|---|
| 518 | else if(j==nk-1) {c1 = 0. |
|---|
| 519 | c2 = 0. |
|---|
| 520 | c3 = abd(3,j+1)*c0 } |
|---|
| 521 | else if(j==nk) {c1 = 0. |
|---|
| 522 | c2 = 0. |
|---|
| 523 | c3 = 0.} |
|---|
| 524 | |
|---|
| 525 | |
|---|
| 526 | p1ip(1,j) = 0. - (c1*wjm3(1)+c2*wjm3(2)+c3*wjm3(3)) |
|---|
| 527 | p1ip(2,j) = 0. - (c1*wjm3(2)+c2*wjm2(1)+c3*wjm2(2)) |
|---|
| 528 | p1ip(3,j) = 0. - (c1*wjm3(3)+c2*wjm2(2)+c3*wjm1(1)) |
|---|
| 529 | |
|---|
| 530 | p1ip(4,j) = c0**2 + |
|---|
| 531 | c1**2*wjm3(1)+2.*c1*c2*wjm3(2)+2.*c1*c3*wjm3(3) + |
|---|
| 532 | c2**2*wjm2(1)+2.*c2*c3*wjm2(2) + |
|---|
| 533 | c3**2*wjm1(1) |
|---|
| 534 | |
|---|
| 535 | wjm3(1)=wjm2(1) ; wjm3(2)=wjm2(2) ; wjm3(3)=p1ip(2,j) |
|---|
| 536 | wjm2(1)=wjm1(1) ; wjm2(2)=p1ip(3,j) |
|---|
| 537 | wjm1(1)=p1ip(4,j) |
|---|
| 538 | |
|---|
| 539 | } |
|---|
| 540 | |
|---|
| 541 | |
|---|
| 542 | if(flag==0) {return} |
|---|
| 543 | |
|---|
| 544 | |
|---|
| 545 | # Pass 2 |
|---|
| 546 | |
|---|
| 547 | |
|---|
| 548 | else { # Compute p2ip |
|---|
| 549 | |
|---|
| 550 | do i=1,nk { j=nk-i+1 |
|---|
| 551 | for(k=1;k<=4 & j+k-1<=nk;k=k+1) |
|---|
| 552 | { p2ip(j,j+k-1) = p1ip(5-k,j) }} |
|---|
| 553 | |
|---|
| 554 | |
|---|
| 555 | do i=1,nk { j=nk-i+1 |
|---|
| 556 | for(k=j-4;k>=1;k=k-1){ |
|---|
| 557 | |
|---|
| 558 | c0 = 1./abd(4,k) ; c1 = abd(1,k+3)*c0 |
|---|
| 559 | c2 = abd(2,k+2)*c0 ; c3 = abd(3,k+1)*c0 |
|---|
| 560 | |
|---|
| 561 | p2ip(k,j) = 0. - ( c1*p2ip(k+3,j) + |
|---|
| 562 | c2*p2ip(k+2,j) + |
|---|
| 563 | c3*p2ip(k+1,j) ) } |
|---|
| 564 | } |
|---|
| 565 | |
|---|
| 566 | return |
|---|
| 567 | |
|---|
| 568 | } |
|---|
| 569 | |
|---|
| 570 | |
|---|
| 571 | end |
|---|
| 572 | subroutine sslvrg(x,y,w,n,knot,nk, |
|---|
| 573 | coef,sz,lev, |
|---|
| 574 | crit,icrit, |
|---|
| 575 | spar,ratio, |
|---|
| 576 | xwy, |
|---|
| 577 | hs0,hs1,hs2,hs3, |
|---|
| 578 | sg0,sg1,sg2,sg3, |
|---|
| 579 | abd,p1ip,p2ip,ld4,ldnk,info) |
|---|
| 580 | |
|---|
| 581 | |
|---|
| 582 | real x(n),y(n),w(n), |
|---|
| 583 | knot(nk+4), |
|---|
| 584 | coef(nk),sz(n),lev(n), |
|---|
| 585 | crit, |
|---|
| 586 | ratio,spar, |
|---|
| 587 | xwy(nk), |
|---|
| 588 | hs0(nk),hs1(nk),hs2(nk),hs3(nk), |
|---|
| 589 | sg0(nk),sg1(nk),sg2(nk),sg3(nk), |
|---|
| 590 | abd(ld4,nk),p1ip(ld4,nk),p2ip(ldnk,nk), |
|---|
| 591 | lambda, |
|---|
| 592 | b0,b1,b2,b3,eps, |
|---|
| 593 | vnikx(4,1),work(16), |
|---|
| 594 | # xv,bvalu,rss,df |
|---|
| 595 | xv,bvalu2,rss,df |
|---|
| 596 | |
|---|
| 597 | |
|---|
| 598 | |
|---|
| 599 | integer n,nk,icrit,ld4,ldnk,i,icoef,ileft,ilo,info,j,mflag |
|---|
| 600 | |
|---|
| 601 | ilo = 1 ; eps = .1e-10 |
|---|
| 602 | |
|---|
| 603 | # Purpose : Solves the smoothing problem and computes the |
|---|
| 604 | # criterion function (OCV or GCV). |
|---|
| 605 | |
|---|
| 606 | |
|---|
| 607 | # The coeficients of estimated smooth |
|---|
| 608 | |
|---|
| 609 | lambda = ratio*16.**(-2. + spar*(6.)) |
|---|
| 610 | |
|---|
| 611 | do i=1,nk { coef(i) = xwy(i) } |
|---|
| 612 | do i=1,nk { abd(4,i) = hs0(i)+lambda*sg0(i) } |
|---|
| 613 | do i=1,(nk-1) { abd(3,i+1) = hs1(i)+lambda*sg1(i) } |
|---|
| 614 | do i=1,(nk-2) { abd(2,i+2) = hs2(i)+lambda*sg2(i) } |
|---|
| 615 | do i=1,(nk-3) { abd(1,i+3) = hs3(i)+lambda*sg3(i) } |
|---|
| 616 | |
|---|
| 617 | |
|---|
| 618 | call spbfa(abd,ld4,nk,3,info) |
|---|
| 619 | if(info.ne.0) { return } |
|---|
| 620 | call spbsl(abd,ld4,nk,3,coef) |
|---|
| 621 | |
|---|
| 622 | |
|---|
| 623 | # Value of smooth at the data points |
|---|
| 624 | |
|---|
| 625 | icoef = 1 |
|---|
| 626 | do i=1,n { xv = x(i) |
|---|
| 627 | #sz(i) = bvalu(knot,coef, |
|---|
| 628 | # nk,4,0,xv,icoef,work(1)) } |
|---|
| 629 | sz(i) = bvalu2(knot,coef, |
|---|
| 630 | nk,4,xv,0) |
|---|
| 631 | } |
|---|
| 632 | |
|---|
| 633 | |
|---|
| 634 | # Compute the criterion function if requested |
|---|
| 635 | |
|---|
| 636 | if(icrit==0) { return} |
|---|
| 637 | |
|---|
| 638 | else { # Ordinary or Generalized CV |
|---|
| 639 | |
|---|
| 640 | # Get Leverages First |
|---|
| 641 | |
|---|
| 642 | call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0) |
|---|
| 643 | |
|---|
| 644 | do i=1,n { xv = x(i) |
|---|
| 645 | #call intrv(knot(1),(nk+1),xv,ilo,ileft,mflag) |
|---|
| 646 | call interv(knot(1),(nk+1),xv,ileft,mflag) |
|---|
| 647 | |
|---|
| 648 | if(mflag==-1) { ileft = 4 ; xv = knot(4)+eps } |
|---|
| 649 | if(mflag==1) { ileft = nk ; xv = knot(nk+1)-eps } |
|---|
| 650 | j=ileft-3 |
|---|
| 651 | #call bspvd(knot,4,1,xv,ileft,4,vnikx,work) |
|---|
| 652 | call bsplvd(knot,4,xv,ileft,work,vnikx,1) |
|---|
| 653 | |
|---|
| 654 | b0=vnikx(1,1);b1=vnikx(2,1);b2=vnikx(3,1);b3=vnikx(4,1) |
|---|
| 655 | |
|---|
| 656 | lev(i) = (p1ip(4,j)*b0**2 + 2.*p1ip(3,j)*b0*b1 + |
|---|
| 657 | 2.*p1ip(2,j)*b0*b2 + 2.*p1ip(1,j)*b0*b3 + |
|---|
| 658 | p1ip(4,j+1)*b1**2 + 2.*p1ip(3,j+1)*b1*b2 + |
|---|
| 659 | 2.*p1ip(2,j+1)*b1*b3 + |
|---|
| 660 | p1ip(4,j+2)*b2**2 + 2.*p1ip(3,j+2)*b2*b3 + |
|---|
| 661 | p1ip(4,j+3)*b3**2 )*w(i)**2 } |
|---|
| 662 | |
|---|
| 663 | |
|---|
| 664 | |
|---|
| 665 | |
|---|
| 666 | |
|---|
| 667 | # Evaluate Criterion |
|---|
| 668 | |
|---|
| 669 | if(icrit==1) { # Generalized CV |
|---|
| 670 | |
|---|
| 671 | rss = 0. ; df = 0. |
|---|
| 672 | do i=1,n { rss = rss + ((y(i)-sz(i))*w(i))**2} |
|---|
| 673 | do i=1,n { df = df + 1.-lev(i)} |
|---|
| 674 | crit = (rss/n)/((df/n)**2) } |
|---|
| 675 | |
|---|
| 676 | else { # Ordinary CV |
|---|
| 677 | |
|---|
| 678 | crit = 0. |
|---|
| 679 | do i=1,n { crit = crit + |
|---|
| 680 | (((y(i)-sz(i))*w(i))/(1-lev(i)))**2 }} |
|---|
| 681 | |
|---|
| 682 | |
|---|
| 683 | return } |
|---|
| 684 | |
|---|
| 685 | end |
|---|
| 686 | subroutine stxwx(x,z,w,k,xknot,n,y,hs0,hs1,hs2,hs3) |
|---|
| 687 | |
|---|
| 688 | real z(k),w(k),x(k), |
|---|
| 689 | xknot(n+4), |
|---|
| 690 | y(n), |
|---|
| 691 | hs0(n),hs1(n),hs2(n),hs3(n), |
|---|
| 692 | |
|---|
| 693 | eps,vnikx(4,1),work(16) # local |
|---|
| 694 | |
|---|
| 695 | integer k,n, |
|---|
| 696 | |
|---|
| 697 | j,i,ilo,ileft,mflag # local |
|---|
| 698 | |
|---|
| 699 | |
|---|
| 700 | |
|---|
| 701 | |
|---|
| 702 | # Initialise the output vectors |
|---|
| 703 | |
|---|
| 704 | do i=1,n { y(i)=0. ; hs0(i)=0. ; hs1(i)=0. |
|---|
| 705 | hs2(i)=0. ; hs3(i)=0. } |
|---|
| 706 | |
|---|
| 707 | |
|---|
| 708 | # Compute X'WX -> hs0,hs1,hs2,hs3 and X'WZ -> y |
|---|
| 709 | |
|---|
| 710 | ilo=1 ; eps = .1e-9 |
|---|
| 711 | |
|---|
| 712 | do i=1,k { |
|---|
| 713 | |
|---|
| 714 | #call intrv(xknot(1),(n+1),x(i),ilo,ileft,mflag) |
|---|
| 715 | call interv(xknot(1),(n+1),x(i),ileft,mflag) |
|---|
| 716 | # check |
|---|
| 717 | # print 999;999 format(" got through check stxwx1") |
|---|
| 718 | |
|---|
| 719 | |
|---|
| 720 | if(mflag==-1) {write(6,'("Error in hess ",i2)')mflag;stop} |
|---|
| 721 | if(mflag== 1) {if(x(i)<=(xknot(ileft)+eps)){ileft=ileft-1} |
|---|
| 722 | else{write(6,'("Error in hess ",i2)')mflag;stop}} |
|---|
| 723 | |
|---|
| 724 | |
|---|
| 725 | #call bspvd (xknot,4,1,x(i),ileft,4,vnikx,work) |
|---|
| 726 | call bsplvd (xknot,4,x(i),ileft,work,vnikx,1) |
|---|
| 727 | # check 2 |
|---|
| 728 | # print 2999;2999 format(" got through check2 stxwx ") |
|---|
| 729 | |
|---|
| 730 | |
|---|
| 731 | |
|---|
| 732 | j= ileft-4+1 |
|---|
| 733 | y(j) = y(j)+w(i)**2*z(i)*vnikx(1,1) |
|---|
| 734 | hs0(j)=hs0(j)+w(i)**2*vnikx(1,1)**2 |
|---|
| 735 | hs1(j)=hs1(j)+w(i)**2*vnikx(1,1)*vnikx(2,1) |
|---|
| 736 | hs2(j)=hs2(j)+w(i)**2*vnikx(1,1)*vnikx(3,1) |
|---|
| 737 | hs3(j)=hs3(j)+w(i)**2*vnikx(1,1)*vnikx(4,1) |
|---|
| 738 | |
|---|
| 739 | j= ileft-4+2 |
|---|
| 740 | y(j) = y(j)+w(i)**2*z(i)*vnikx(2,1) |
|---|
| 741 | hs0(j)=hs0(j)+w(i)**2*vnikx(2,1)**2 |
|---|
| 742 | hs1(j)=hs1(j)+w(i)**2*vnikx(2,1)*vnikx(3,1) |
|---|
| 743 | hs2(j)=hs2(j)+w(i)**2*vnikx(2,1)*vnikx(4,1) |
|---|
| 744 | |
|---|
| 745 | j= ileft-4+3 |
|---|
| 746 | y(j) = y(j)+w(i)**2*z(i)*vnikx(3,1) |
|---|
| 747 | hs0(j)=hs0(j)+w(i)**2*vnikx(3,1)**2 |
|---|
| 748 | hs1(j)=hs1(j)+w(i)**2*vnikx(3,1)*vnikx(4,1) |
|---|
| 749 | |
|---|
| 750 | j= ileft-4+4 |
|---|
| 751 | y(j) = y(j)+w(i)**2*z(i)*vnikx(4,1) |
|---|
| 752 | hs0(j)=hs0(j)+w(i)**2*vnikx(4,1)**2 } |
|---|
| 753 | |
|---|
| 754 | return |
|---|
| 755 | end |
|---|
| 756 | subroutine sknotl(x,n,knot,k) |
|---|
| 757 | |
|---|
| 758 | real x(n),knot(n+6),a1,a2,a3,a4 |
|---|
| 759 | integer n,k,ndk,j |
|---|
| 760 | |
|---|
| 761 | |
|---|
| 762 | # Allocate knots acording to the cutoffs given below |
|---|
| 763 | |
|---|
| 764 | |
|---|
| 765 | # Cutoff constants |
|---|
| 766 | |
|---|
| 767 | a1 = log(50.)/log(2.) ; a2 = log(100.)/log(2.) |
|---|
| 768 | a3 = log(140.)/log(2.) ; a4 = log(200.)/log(2.) |
|---|
| 769 | |
|---|
| 770 | # Cutoff Criteria |
|---|
| 771 | |
|---|
| 772 | if(n<50) { ndk = n } |
|---|
| 773 | else if (n>=50 & n<200) { ndk = 2.**(a1+(a2-a1)*(n-50.)/150.) } |
|---|
| 774 | else if (n>=200 & n<800) { ndk = 2.**(a2+(a3-a2)*(n-200.)/600.) } |
|---|
| 775 | else if (n>=800 & n<3200) { ndk = 2.**(a3+(a4-a3)*(n-800.)/2400.) } |
|---|
| 776 | else if (n>=3200) { ndk = 200. + (n-3200)**.2 } |
|---|
| 777 | |
|---|
| 778 | k = ndk + 6 |
|---|
| 779 | |
|---|
| 780 | |
|---|
| 781 | # Allocate Knots ( note no account is taken of any weighting vector ) |
|---|
| 782 | |
|---|
| 783 | do j=1,3 { knot(j) = x(1) } |
|---|
| 784 | do j=1,ndk { knot(j+3) = x( 1 + (j-1)*(n-1)/(ndk-1) ) } |
|---|
| 785 | do j=1,3 { knot(ndk+3+j) = x(n) } |
|---|
| 786 | |
|---|
| 787 | return |
|---|
| 788 | end |
|---|
| 789 | subroutine setreg(x,y,w,n,xw,nx,min,range,knot,nk) |
|---|
| 790 | |
|---|
| 791 | # sbart uses the square of the weights; the following rectifies that. |
|---|
| 792 | # Also, the data is sorted (by a routine you may need to change for |
|---|
| 793 | # your local machine) and standardized so that spar=1 gives a smooth fit |
|---|
| 794 | # with sum(leverage values) about 2-4. |
|---|
| 795 | |
|---|
| 796 | real x(n),y(n),w(n),xw(n), |
|---|
| 797 | min,range,knot(n+6) |
|---|
| 798 | |
|---|
| 799 | integer n,nk,nx |
|---|
| 800 | |
|---|
| 801 | |
|---|
| 802 | # Local |
|---|
| 803 | |
|---|
| 804 | real eps |
|---|
| 805 | integer ycnt,i,k |
|---|
| 806 | |
|---|
| 807 | |
|---|
| 808 | call scopy(n,x,1,xw,1) |
|---|
| 809 | #call ssort(x,w,n,2) |
|---|
| 810 | #call ssort(xw,y,n,2) |
|---|
| 811 | call sortpr(x,n,w) |
|---|
| 812 | call sortpr(xw,n,y) |
|---|
| 813 | |
|---|
| 814 | range = x(n)-x(1) ; min = x(1) ; eps = .1e-9 |
|---|
| 815 | do i=1,n { x(i) = (x(i)-min)/range } |
|---|
| 816 | call scopy(n,x,1,xw,1) |
|---|
| 817 | |
|---|
| 818 | |
|---|
| 819 | nx = 1 ; x(nx) = x(1) ; w(nx) = w(1) ; y(nx) = y(1)*w(1) ; |
|---|
| 820 | |
|---|
| 821 | for(i=2;i<=n;i=i+1) |
|---|
| 822 | |
|---|
| 823 | { if(xw(i)>x(nx)+eps) |
|---|
| 824 | { if(w(nx)>0.0) y(nx) = y(nx)/w(nx) |
|---|
| 825 | nx = nx + 1 |
|---|
| 826 | x(nx) = x(i) |
|---|
| 827 | y(nx) = y(i)*w(i) |
|---|
| 828 | w(nx) = w(i) } |
|---|
| 829 | else |
|---|
| 830 | { y(nx) = y(nx)+y(i)*w(i) |
|---|
| 831 | w(nx) = w(nx) + w(i) } |
|---|
| 832 | } |
|---|
| 833 | if(w(nx)>0.0) y(nx) = y(nx)/w(nx) |
|---|
| 834 | |
|---|
| 835 | for(i=1;i<=nx;i=i+1) |
|---|
| 836 | { if (w(i)>0) w(i)=sqrt(w(i)) } |
|---|
| 837 | |
|---|
| 838 | |
|---|
| 839 | call sknotl(x,nx,knot,k) ; nk = k-4 |
|---|
| 840 | |
|---|
| 841 | |
|---|
| 842 | |
|---|
| 843 | |
|---|
| 844 | return |
|---|
| 845 | end |
|---|
| 846 | |
|---|
| 847 | |
|---|
