R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # splsmu.r smoothing cubic splines > # > all <- read.table("eppright.dat") > x <- all[,1] > y <- all[,2] > # direct calculation > n <- length(x) > nm1 <- n-1 > h <- diff(x) > t <- 72*1.38 > Estar <- diag(2*(h[-1]+h[-nm1]),n-2) > diag(Estar[2:(n-2),1:(n-3)]) <- h[2:(n-2)] # subdiagonal > diag(Estar[1:(n-3),2:(n-2)]) <- h[2:(n-2)] # superdiagonal > Estar[1:6,1:6] [,1] [,2] [,3] [,4] [,5] [,6] [1,] 4 1 0 0 0 0 [2,] 1 4 1 0 0 0 [3,] 0 1 4 1 0 0 [4,] 0 0 1 4 1 0 [5,] 0 0 0 1 4 1 [6,] 0 0 0 0 1 4 > C <- matrix(0,n-2,n) > diag(C[,1:(n-2)]) <- 1/h[-nm1] > diag(C[,3:n]) <- 1/h[-1] > diag(C[,2:nm1]) <- -(1/h[-1]+1/h[-nm1]) > C[1:6,1:8] [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 1 -2 1 0 0 0 0 0 [2,] 0 1 -2 1 0 0 0 0 [3,] 0 0 1 -2 1 0 0 0 [4,] 0 0 0 1 -2 1 0 0 [5,] 0 0 0 0 1 -2 1 0 [6,] 0 0 0 0 0 1 -2 1 > Lt <- chol(Estar+6*t*C%*%t(C)) > Lt[1:6,1:6] [,1] [,2] [,3] [,4] [,5] [,6] [1,] 59.84112 -39.83281 9.96238 0.00000 0.00000 0.00000 [2,] 0.00000 44.65767 -44.48979 13.34955 0.00000 0.00000 [3,] 0.00000 0.00000 38.76041 -46.17395 15.38064 0.00000 [4,] 0.00000 0.00000 0.00000 35.64710 -46.94505 16.72394 [5,] 0.00000 0.00000 0.00000 0.00000 33.77215 -47.33290 [6,] 0.00000 0.00000 0.00000 0.00000 0.00000 32.57095 > mstar <- backsolve(Lt,forwardsolve(t(Lt),6*C%*%y)) > z <- y - t*t(C)%*%mstar > smuspl <- smooth.spline(x,y,spar=0) > smuspl$lambda [1] 5.205431e-09 > smuspl$lev [1] 0.9970404 0.9804708 0.9391799 0.7495909 0.5627758 0.7519262 0.9096992 [8] 0.7487089 0.5577719 0.7487420 0.9092659 0.7487089 0.5577242 0.7487090 [15] 0.9092645 0.7487319 0.5577571 0.7487088 0.9095706 0.7509532 0.5610125 [22] 0.7487040 0.9368881 0.9368881 0.7487040 0.5610125 0.7509532 0.9095706 [29] 0.7487088 0.5577571 0.7487319 0.9092645 0.7487087 0.5577237 0.7487089 [36] 0.9092613 0.7487089 0.5577237 0.7487087 0.9092645 0.7487319 0.5577571 [43] 0.7487088 0.9095706 0.7509532 0.5610125 0.7487040 0.9368881 0.9368881 [50] 0.7487040 0.5610125 0.7509532 0.9095706 0.7487088 0.5577571 0.7487319 [57] 0.9092645 0.7487089 0.5577240 0.7487089 0.9092644 0.7487310 0.5577558 [64] 0.7487088 0.9095586 0.7508659 0.5608849 0.7487318 0.9364633 0.9576007 [71] 0.9644304 0.9946130 > smuspl <- smooth.spline(x,y,spar=.654278) > smuspl$lambda [1] 0.000277608 > smuspl$lev [1] 0.3610481 0.2366454 0.1684313 0.1353979 0.1213059 0.1175706 0.1168150 [8] 0.1168140 0.1160152 0.1159138 0.1150429 0.1141586 0.1128514 0.1128323 [15] 0.1124466 0.1122112 0.1115292 0.1120187 0.1119927 0.1119849 0.1114337 [22] 0.1119839 0.1119832 0.1119820 0.1119806 0.1114288 0.1119770 0.1119755 [29] 0.1119748 0.1114217 0.1119740 0.1119736 0.1119738 0.1114212 0.1119738 [36] 0.1119736 0.1119738 0.1114212 0.1119738 0.1119736 0.1119739 0.1114215 [43] 0.1119744 0.1119748 0.1119759 0.1114274 0.1119789 0.1119803 0.1119820 [50] 0.1119834 0.1114337 0.1119834 0.1119844 0.1119927 0.1114661 0.1120816 [57] 0.1122109 0.1124468 0.1122808 0.1134020 0.1141582 0.1150430 0.1153669 [64] 0.1165616 0.1168137 0.1168150 0.1170254 0.1218392 0.1353905 0.1684177 [71] 0.2363362 0.3610144 > cbind(x,y,smuspl$y,z,c(0,mstar,0)) x y [1,] 0.5 0.46 0.4847770 0.4847772 0.000000e+00 [2,] 1.5 0.47 0.5196523 0.5196533 -2.493685e-04 [3,] 2.5 0.56 0.5541952 0.5541967 -9.984684e-04 [4,] 3.5 0.61 0.5877492 0.5877514 -1.689162e-03 [5,] 4.5 0.61 0.6196583 0.6196543 -2.155937e-03 [6,] 5.5 0.67 0.6493827 0.6493850 -2.719876e-03 [7,] 6.5 0.68 0.6764276 0.6764304 -3.076337e-03 [8,] 7.5 0.78 0.7004022 0.7004054 -3.396872e-03 [9,] 8.5 0.69 0.7211344 0.7211172 -2.916334e-03 [10,] 9.5 0.74 0.7388557 0.7388603 -2.748973e-03 [11,] 10.5 0.77 0.7538519 0.7538564 -2.570141e-03 [12,] 11.5 0.78 0.7663050 0.7663095 -2.228834e-03 [13,] 12.5 0.75 0.7765698 0.7765567 -1.749739e-03 [14,] 13.5 0.80 0.7850053 0.7850096 -1.537923e-03 [15,] 14.5 0.78 0.7919461 0.7919497 -1.175236e-03 [16,] 15.5 0.82 0.7976918 0.7976945 -9.328158e-04 [17,] 16.5 0.77 0.8025651 0.8025440 -4.659043e-04 [18,] 17.5 0.80 0.8068717 0.8068729 -3.265286e-04 [19,] 18.5 0.81 0.8108638 0.8108638 -2.563248e-04 [20,] 19.5 0.78 0.8145980 0.8145969 -1.948146e-04 [21,] 20.5 0.87 0.8180422 0.8180772 -4.815019e-04 [22,] 21.5 0.80 0.8211652 0.8211630 -2.456166e-04 [23,] 22.5 0.83 0.8239699 0.8239678 -2.227249e-04 [24,] 23.5 0.81 0.8265617 0.8265599 -1.391224e-04 [25,] 24.5 0.88 0.8289867 0.8289852 -2.221858e-04 [26,] 25.5 0.81 0.8312896 0.8312738 2.081853e-04 [27,] 26.5 0.83 0.8337352 0.8337349 4.244482e-04 [28,] 27.5 0.82 0.8366144 0.8366142 6.031214e-04 [29,] 28.5 0.82 0.8400690 0.8400688 6.145824e-04 [30,] 29.5 0.86 0.8440931 0.8441043 4.240629e-04 [31,] 30.5 0.82 0.8485909 0.8485905 3.935247e-04 [32,] 31.5 0.85 0.8534226 0.8534222 7.524038e-05 [33,] 32.5 0.88 0.8583238 0.8583235 -2.774868e-04 [34,] 33.5 0.86 0.8629863 0.8629837 -4.120528e-04 [35,] 34.5 0.91 0.8672272 0.8672267 -5.766475e-04 [36,] 35.5 0.87 0.8709652 0.8709649 -3.107546e-04 [37,] 36.5 0.87 0.8743909 0.8743908 -5.457314e-05 [38,] 37.5 0.87 0.8777602 0.8777546 1.574180e-04 [39,] 38.5 0.85 0.8812627 0.8812629 2.913633e-04 [40,] 39.5 0.90 0.8850101 0.8850101 1.106658e-04 [41,] 40.5 0.87 0.8888932 0.8888931 8.083247e-05 [42,] 41.5 0.91 0.8928130 0.8928253 -1.391493e-04 [43,] 42.5 0.90 0.8966473 0.8966471 -1.862779e-04 [44,] 43.5 0.93 0.9002881 0.9002883 -1.996619e-04 [45,] 44.5 0.89 0.9037790 0.9037796 8.598494e-05 [46,] 45.5 0.89 0.9073453 0.9073339 2.329477e-04 [47,] 46.5 0.92 0.9110907 0.9110919 2.054553e-04 [48,] 47.5 0.89 0.9150690 0.9150704 2.676172e-04 [49,] 48.5 0.92 0.9192730 0.9192745 7.745996e-05 [50,] 49.5 0.96 0.9235557 0.9235572 -1.053954e-04 [51,] 50.5 0.92 0.9277993 0.9277957 7.852447e-05 [52,] 51.5 0.91 0.9320977 0.9320996 1.839855e-04 [53,] 52.5 0.95 0.9365486 0.9365504 6.702732e-05 [54,] 53.5 0.93 0.9410891 0.9410908 8.543143e-05 [55,] 54.5 0.93 0.9457076 0.9456981 -7.786981e-06 [56,] 55.5 0.98 0.9502701 0.9502712 -2.589970e-04 [57,] 56.5 0.95 0.9546343 0.9546352 -2.110039e-04 [58,] 57.5 0.97 0.9587796 0.9587804 -2.096610e-04 [59,] 58.5 0.97 0.9627282 0.9627347 -9.539908e-05 [60,] 59.5 0.96 0.9666047 0.9666059 9.198337e-05 [61,] 60.5 0.97 0.9705565 0.9705580 2.128813e-04 [62,] 61.5 0.94 0.9747203 0.9747220 3.281637e-04 [63,] 62.5 0.96 0.9791679 0.9791559 9.398990e-05 [64,] 63.5 1.03 0.9836505 0.9836517 -3.329767e-04 [65,] 64.5 0.99 0.9878909 0.9878922 -2.934747e-04 [66,] 65.5 1.01 0.9918414 0.9918428 -2.327592e-04 [67,] 66.5 0.99 0.9955926 0.9955912 1.069740e-05 [68,] 67.5 0.99 0.9993384 0.9993408 1.978823e-04 [69,] 68.5 0.99 1.0032700 1.0032726 2.910577e-04 [70,] 69.5 1.01 1.0074704 1.0074733 2.506519e-04 [71,] 70.5 0.99 1.0119375 1.0119288 2.356762e-04 [72,] 71.5 1.04 1.0165801 1.0165832 0.000000e+00 > # clean up & go > rm(list=ls()) > q()