> # conclk.r September 2007 > # > # concentrated likelihood problem, Example 9.6 > # > x <- c(1,2,3,4,5,7) > y <- c(8.3, 10.3, 19.0, 16.0, 15.6, 19.8) > # construct error sum of squares function > sse <- function(t2) { + z <- 1 - exp(-t2*x) + t1 <- crossprod(y,z)/crossprod(z,z) # regression through origin + e <- y - t1*z # residuals + s <- crossprod(e,e) # error sum of squares + print("t1,t2,sse") + print(c(t1,t2,s)) + sse <- as.real(s) } # make it a scalar > # starting values? > sse(.1) [1] "t1,t2,sse" [1] 46.26346 0.10000 87.49828 > sse(1) [1] "t1,t2,sse" [1] 16.64283 1.00000 42.85488 > sse(5) [1] "t1,t2,sse" [1] 14.85742 5.00000 105.90782 > sse(10) [1] "t1,t2,sse" [1] 14.83350 10.00000 107.20453 > # minimize error sum of squares > that2 <- optimize(f=sse,lower=.1,upper=5) [1] "t1,t2,sse" [1] 15.376096 1.971633 80.639967 [1] "t1,t2,sse" [1] 14.993260 3.128367 98.750177 [1] "t1,t2,sse" [1] 16.099940 1.256733 55.200735 [1] "t1,t2,sse" [1] 17.2634313 0.8149004 33.9910427 [1] "t1,t2,sse" [1] 19.0324039 0.5418327 26.0078259 [1] "t1,t2,sse" [1] 21.5716334 0.3730676 31.0740041 [1] "t1,t2,sse" [1] 18.7704232 0.5693695 26.2033106 [1] "t1,t2,sse" [1] 19.0831309 0.5368292 25.9953196 [1] "t1,t2,sse" [1] 19.1488920 0.5304894 25.9903235 [1] "t1,t2,sse" [1] 19.1430332 0.5310477 25.9902676 [1] "t1,t2,sse" [1] 19.1425429 0.5310945 25.9902673 [1] "t1,t2,sse" [1] 19.1421164 0.5311352 25.9902676 [1] "t1,t2,sse" [1] 19.1425429 0.5310945 25.9902673 > # print it out > print("min sse") [1] "min sse" > print(that2) $minimum [1] 0.5310945 $objective [1] 25.99027 > # done > rm(list=ls()) > q() > # concll.r October 2007, 2008 > # > # concentrated likelihood problem, Example 9.6 > # redone as full three parameter problem > # > x <- c(1,2,3,4,5,7) > y <- c(8.3, 10.3, 19.0, 16.0, 15.6, 19.8) > # construct error sum of squares function > ull <- function(t) { + e <- y - t[1]*(1 - exp(-t[2]*x) ) + ee <- crossprod(e,e) + print(c(t,ee)) + ull <- (length(y)*log(t[3]) + crossprod(e,e)/t[3] ) / 2 } > # > # minimize unlikelihood > that2 <- nlm(ull, c(10,1,6), hessian=T, print.level=1) [1] 10.0000 1.0000 6.0000 263.4282 [1] 10.0000 1.0000 6.0000 263.4282 [1] 10.00001 1.00000 6.00000 263.42751 [1] 10.000000 1.000001 6.000000 263.428109 [1] 10.000000 1.000000 6.000006 263.428175 iteration = 0 Step: [1] 0 0 0 Parameter: [1] 10 1 6 Function Value [1] 27.32763 Gradient: [1] -5.534114 -5.501032 -3.158721 [1] 15.534114 6.501032 9.158721 109.822519 [1] 15.534129 6.501032 9.158721 109.822648 [1] 15.534114 6.501038 9.158721 109.822521 [1] 15.534114 6.501032 9.158730 109.822519 [1] 15.116230 6.501048 9.481728 107.384212 [1] 15.116245 6.501048 9.481728 107.384262 [1] 15.116230 6.501055 9.481728 107.384214 [1] 15.116230 6.501048 9.481737 107.384212 [1] 14.648473 6.492637 10.198840 107.136882 [1] 14.648488 6.492637 10.198840 107.136848 [1] 14.648473 6.492643 10.198840 107.136883 [1] 14.648473 6.492637 10.198850 107.136882 [1] 14.448676 6.476423 11.010802 107.827808 [1] 14.448690 6.476423 11.010802 107.827741 [1] 14.44868 6.47643 11.01080 107.82781 [1] 14.448676 6.476423 11.010813 107.827808 [1] 14.322166 6.432403 12.755313 108.503836 [1] 14.322180 6.432403 12.755313 108.503747 [1] 14.32217 6.43241 12.75531 108.50384 [1] 14.322166 6.432403 12.755326 108.503836 [1] 14.401734 6.372301 14.694728 108.031042 *** about 130 lines deleted *** [1] 19.1445791 0.5295839 4.3557543 25.9913521 [1] 19.1445983 0.5295839 4.3557543 25.9913501 [1] 19.1445791 0.5295849 4.3557543 25.9913505 [1] 19.1445791 0.5295839 4.3557586 25.9913521 [1] 19.1420509 0.5312487 4.3291098 25.9902760 [1] 19.1420700 0.5312487 4.3291098 25.9902761 [1] 19.1420509 0.5312497 4.3291098 25.9902761 [1] 19.1420509 0.5312487 4.3291141 25.9902760 [1] 19.1424321 0.5310991 4.3316842 25.9902673 [1] 19.1424512 0.5310991 4.3316842 25.9902673 [1] 19.1424321 0.5311001 4.3316842 25.9902673 [1] 19.1424321 0.5310991 4.3316885 25.9902673 [1] 19.1425303 0.5310933 4.3316855 25.9902673 [1] 19.1425495 0.5310933 4.3316855 25.9902673 [1] 19.1425303 0.5310943 4.3316855 25.9902673 [1] 19.1425303 0.5310933 4.3316899 25.9902673 [1] 19.142559 0.531092 4.331708 25.990267 [1] 19.142578 0.531092 4.331708 25.990267 [1] 19.142559 0.531093 4.331708 25.990267 [1] 19.142559 0.531092 4.331713 25.990267 iteration = 37 Parameter: [1] 19.142559 0.531092 4.331708 Function Value [1] 7.397888 Gradient: [1] 3.711848e-09 4.796163e-08 -9.308868e-08 Relative gradient close to zero. Current iterate is probably solution. [1] 19.144473 0.531092 4.331708 25.990281 [1] 19.142559 0.531192 4.331708 25.990273 [1] 19.142559 0.531092 4.332142 25.990267 [1] 19.146388 0.531092 4.331708 25.990323 [1] 19.144473 0.531192 4.331708 25.990302 [1] 19.144473 0.531092 4.332142 25.990281 [1] 19.142559 0.531292 4.331708 25.990290 [1] 19.142559 0.531192 4.332142 25.990273 [1] 19.142559 0.531092 4.332575 25.990267 > # print it out > print("minimize unlikelihood") [1] "minimize unlikelihood" > print(that2) $minimum [1] 7.397888 $estimate [1] 19.142559 0.531092 4.331708 $gradient [1] 3.711848e-09 4.796163e-08 -9.308868e-08 $hessian [,1] [,2] [,3] [1,] 0.8828751224 9.253518024 -0.0001931081 [2,] 9.2535180238 132.650704998 -0.0015164637 [3,] -0.0001931081 -0.001516464 0.1598194705 $code [1] 1 $iterations [1] 37 > # > # inverse hessian -- kiss > solve(that2$hessian) [,1] [,2] [,3] [1,] 4.212951713 -2.938893e-01 2.301865e-03 [2,] -0.293889287 2.803988e-02 -8.904394e-05 [3,] 0.002301865 -8.904394e-05 6.257062e+00 > # > warnings() Warning messages: 1: NaNs produced in: log(x) 2: NA/Inf replaced by maximum positive value 3: NaNs produced in: log(x) 4: NA/Inf replaced by maximum positive value 5: NaNs produced in: log(x) 6: NA/Inf replaced by maximum positive value 7: NaNs produced in: log(x) 8: NA/Inf replaced by maximum positive value 9: NaNs produced in: log(x) 10: NA/Inf replaced by maximum positive value 11: NaNs produced in: log(x) 12: NA/Inf replaced by maximum positive value 13: NaNs produced in: log(x) 14: NA/Inf replaced by maximum positive value 15: NaNs produced in: log(x) 16: NA/Inf replaced by maximum positive value 17: NaNs produced in: log(x) 18: NA/Inf replaced by maximum positive value > # > # > # now as nonlinear least squares problem > bwnlls <- list(y,x) # need list or data frame > that1 <- nls( y ~ th1-th1*exp(-th2*x), data=bwnlls, + start=list( th1=10, th2=1 ), trace=T ) # same start 263.4282 : 10 1 155.1979 : 14.1213295 0.4448702 26.15580 : 19.0199864 0.5566534 25.99084 : 19.1304520 0.5310812 25.99027 : 19.1425771 0.5310912 > print(that1) Nonlinear regression model model: y ~ th1 - th1 * exp(-th2 * x) data: bwnlls th1 th2 19.1425771 0.5310912 residual sum-of-squares: 25.99027 > summary(that1) # prints different stuff Formula: y ~ th1 - th1 * exp(-th2 * x) Parameters: Estimate Std. Error t value Pr(>|t|) th1 19.1426 2.4959 7.670 0.00155 ** th2 0.5311 0.2031 2.615 0.05910 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.549 on 4 degrees of freedom Correlation of Parameter Estimates: th1 th2 -0.8528 > # done > rm(list=ls()) > q()