R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), 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 a HTML browser interface to help. Type 'q()' to quit R. > # 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 [1] 14.401748 6.372301 14.694728 108.030967 [1] 14.401734 6.372307 14.694728 108.031044 [1] 14.401734 6.372301 14.694742 108.031042 [1] 14.592351 6.305808 16.449335 107.226793 [1] 14.592366 6.305808 16.449335 107.226750 [1] 14.592351 6.305814 16.449335 107.226796 [1] 14.592351 6.305808 16.449352 107.226793 [1] 14.750441 6.254439 17.493297 106.889154 [1] 14.750456 6.254439 17.493297 106.889138 [1] 14.750441 6.254445 17.493297 106.889156 [1] 14.750441 6.254439 17.493315 106.889154 [1] 14.819872 6.225081 17.855594 106.832289 [1] 14.819887 6.225081 17.855594 106.832285 [1] 14.819872 6.225087 17.855594 106.832291 [1] 14.819872 6.225081 17.855612 106.832289 [1] 14.83798 6.20838 17.90882 106.82334 [1] 14.83799 6.20838 17.90882 106.82334 [1] 14.837977 6.208386 17.908818 106.823347 [1] 14.83798 6.20838 17.90884 106.82334 [1] 14.861117 6.169496 17.933183 106.810320 [1] 14.861132 6.169496 17.933183 106.810323 [1] 14.861117 6.169502 17.933183 106.810322 [1] 14.861117 6.169496 17.933201 106.810320 [1] 14.901849 6.058179 17.915719 106.781813 [1] 14.901864 6.058179 17.915719 106.781823 [1] 14.901849 6.058185 17.915719 106.781815 [1] 14.901849 6.058179 17.915737 106.781813 [1] 14.984075 5.726317 17.800526 106.697751 [1] 14.984090 5.726317 17.800526 106.697776 [1] 14.984075 5.726323 17.800526 106.697755 [1] 14.984075 5.726317 17.800544 106.697751 [1] 15.786776 2.005278 16.574199 82.555002 [1] 15.786791 2.005278 16.574199 82.555080 [1] 15.786776 2.005280 16.574199 82.555058 [1] 15.786776 2.005278 16.574216 82.555002 [1] 99.87009 -379.60724 -104.01424 Inf [1] 2.419511e+01 -3.615597e+01 4.515356e+00 3.983046e+222 [1] 1.662761e+01 -1.810847e+00 1.536831e+01 2.832376e+13 [1] 15.870859 1.623666 16.453611 70.369241 [1] 15.870875 1.623666 16.453611 70.369282 [1] 15.870859 1.623667 16.453611 70.369302 [1] 15.870859 1.623666 16.453627 70.369241 [1] 132.9126 -529.4033 -151.6482 Inf [1] 27.5750379 -51.4790278 -0.3565678 Inf [1] 1.704128e+01 -3.686604e+00 1.477259e+01 7.551045e+24 [1] 15.987901 1.092639 16.285509 48.352928 [1] 15.987917 1.092639 16.285509 48.352858 [1] 15.987901 1.092640 16.285509 48.352971 [1] 15.987901 1.092639 16.285525 48.352928 [1] 145.6813 -587.1515 -171.3320 Inf [1] 28.957238 -57.731773 -2.476238 Inf [1] 1.728483e+01 -4.789802e+00 1.440933e+01 3.962021e+31 [1] 16.1175941 0.5043947 16.0978915 67.0033214 [1] 16.000870 1.033814 16.266747 46.073616 [1] 16.000886 1.033814 16.266747 46.073527 [1] 16.000870 1.033815 16.266747 46.073653 [1] 16.000870 1.033814 16.266764 46.073616 [1] 16.26687155 0.04628018 9.61381495 980.51098188 [1] 16.027470 0.935061 15.601454 42.821521 [1] 16.027486 0.935061 15.601454 42.821395 [1] 16.027470 0.935062 15.601454 42.821547 [1] 16.027470 0.935061 15.601470 42.821521 [1] 16.34764984 -0.05156214 1.66168069 2232.47431954 [1] 16.0594882 0.8363986 14.2074767 40.8507529 [1] 16.0595043 0.8363986 14.2074767 40.8505836 [1] 16.0594882 0.8363996 14.2074767 40.8507596 [1] 16.0594882 0.8363986 14.2074909 40.8507529 [1] 16.25185679 0.33035720 -0.07028211 146.29839733 [1] 16.0787251 0.7857945 12.7797008 40.6410672 [1] 16.0787411 0.7857945 12.7797008 40.6408728 [1] 16.0787251 0.7857955 12.7797008 40.6410594 [1] 16.0787251 0.7857945 12.7797136 40.6410672 [1] 16.2623617 0.4603549 -5.1792737 75.9488943 [1] 16.0970887 0.7532505 10.9838034 40.8457952 [1] 16.0971048 0.7532505 10.9838034 40.8455840 [1] 16.0970887 0.7532515 10.9838034 40.8457758 [1] 16.0970887 0.7532505 10.9838143 40.8457952 [1] 16.3614871 0.5981121 -6.2716268 45.7140811 [1] 16.1235286 0.7377367 9.2582603 40.8383651 [1] 16.1235447 0.7377367 9.2582603 40.8381478 [1] 16.1235286 0.7377377 9.2582603 40.8383400 [1] 16.1235286 0.7377367 9.2582696 40.8383651 [1] 16.4866383 0.6929185 -0.3146262 37.4801050 [1] 16.1598395 0.7332549 8.3009717 40.4655230 [1] 16.1598557 0.7332549 8.3009717 40.4653078 [1] 16.1598395 0.7332559 8.3009717 40.4654974 [1] 16.1598395 0.7332549 8.3009800 40.4655230 [1] 16.7163709 0.6970787 1.5509984 34.8430742 [1] 16.2154927 0.7296373 7.6259744 39.8258487 [1] 16.2155089 0.7296373 7.6259744 39.8256384 [1] 16.2154927 0.7296383 7.6259744 39.8258236 [1] 16.2154927 0.7296373 7.6259820 39.8258487 [1] 17.1535728 0.6723255 0.7938136 31.5018742 [1] 16.3093007 0.7239061 6.9427583 38.7788253 [1] 16.3093170 0.7239061 6.9427583 38.7786237 [1] 16.3093007 0.7239071 6.9427583 38.7788013 [1] 16.3093007 0.7239061 6.9427652 38.7788253 [1] 17.6453741 0.6424766 0.2910435 28.8554698 [1] 16.4429080 0.7157631 6.2775868 37.3736103 [1] 16.4429245 0.7157631 6.2775868 37.3734211 [1] 16.4429080 0.7157641 6.2775868 37.3735879 [1] 16.4429080 0.7157631 6.2775931 37.3736103 [1] 18.0330157 0.6173303 0.6491624 27.4729779 [1] 16.6019188 0.7059198 5.7147444 35.8334484 [1] 16.6019354 0.7059198 5.7147444 35.8332737 [1] 16.6019188 0.7059208 5.7147444 35.8334280 [1] 16.6019188 0.7059198 5.7147501 35.8334484 [1] 18.3618821 0.5940264 1.2614091 26.6998966 [1] 16.8753739 0.6885344 5.0228055 33.5064627 [1] 16.8753907 0.6885344 5.0228055 33.5063125 [1] 16.8753739 0.6885354 5.0228055 33.5064456 [1] 16.8753739 0.6885344 5.0228105 33.5064627 [1] 18.5254349 0.5785011 2.2639079 26.4201652 [1] 17.6254490 0.6385161 3.7686819 28.9782102 [1] 17.6254666 0.6385161 3.7686819 28.9781219 [1] 17.6254490 0.6385171 3.7686819 28.9782016 [1] 17.6254490 0.6385161 3.7686856 28.9782102 [1] 18.466896 0.572645 3.698242 26.539860 [1] 18.466915 0.572645 3.698242 26.539819 [1] 18.466896 0.572646 3.698242 26.539852 [1] 18.466896 0.572645 3.698246 26.539860 [1] 18.8236677 0.5785539 4.2784178 26.3856274 [1] 18.8236865 0.5785539 4.2784178 26.3856478 [1] 18.8236677 0.5785549 4.2784178 26.3856500 [1] 18.8236677 0.5785539 4.2784221 26.3856274 [1] 19.287903 0.492105 4.789499 26.567063 [1] 19.023997 0.541249 4.498963 26.006772 [1] 19.024016 0.541249 4.498963 26.006770 [1] 19.023997 0.541250 4.498963 26.006774 [1] 19.023997 0.541249 4.498967 26.006772 [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()