例題:
「表 1 に示すようなデータに曲線をあてはめなさい。」
x | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
y | 2.9 | 5.2 | 9.1 | 15.5 | 25 | 37.8 | 52.6 | 66.9 | 78.6 | 87 | 92.4 | 95.7 | 97.6 | 98.6 | 99.2 |
プログラム: # quartz(width=500/72, height=350/72) x <- 1:15 y <- c(2.9, 5.2, 9.1, 15.5, 25, 37.8, 52.6, 66.9, 78.6, 87, 92.4, 95.7, 97.6, 98.6, 99.2) dat <- data.frame(x=x, y=y) # http://aoki2.si.gunma-u.ac.jp/R/simplex.html の,シンプレックス法を使う model1 <- function(x, p) return(p[1]*p[2]^exp(-p[3]*x)) simplex(model1, c(100, 0.001, 0.1), x, y) # シンプレックス法で得られた解を初期値として, # eq. 1: y ~ a*b^exp(-c*x) であてはめる (ans <- nls(y ~ a*b^exp(-c*x), start=list(a=105, b=0.00016, c=0.37), data=dat)) summary(ans) plot(x, y, pch=19, col="red") x2 <- seq(1, 15, by=0.01) dat2 <- data.frame(x=x2) y2 <- predict(ans, newdata=dat2) lines(x2, y2, col="blue") # eq. 2: y ~ A*exp(-B*C^x) であてはめる(初期値の心配がない) ans2 <- nls(y ~ SSgompertz(x, A, B, C), data=dat) p <- ans2$m$getPars() p[1] # eq. 1 の a exp(-p[2]) # eq. 1 の b -log(p[3]) # eq. 1 の c 実行例: > x <- 1:15 > y <- c(2.9, 5.2, 9.1, 15.5, 25, 37.8, 52.6, 66.9, + 78.6, 87, 92.4, 95.7, 97.6, 98.6, 99.2) > dat <- data.frame(x=x, y=y) > # http://aoki2.si.gunma-u.ac.jp/R/simplex.html の,シンプレックス法を使う > model1 <- function(x, p) return(p[1]*p[2]^exp(-p[3]*x)) > simplex(model1, c(100, 0.001, 0.1), x, y) $converge [1] TRUE $parameters [1] 1.050828e+02 1.571673e-04 3.699247e-01 $residuals [1] 66.53056 > # シンプレックス法で得られた解を初期値として, > # eq. 1: y ~ a*b^exp(-c*x) であてはめる > (ans <- nls(y ~ a*b^exp(-c*x), start=list(a=105, b=0.00016, c=0.37), data=dat)) Nonlinear regression model model: y ~ a * b^exp(-c * x) data: dat a b c 1.051e+02 1.572e-04 3.699e-01 residual sum-of-squares: 66.53 Number of iterations to convergence: 4 Achieved convergence tolerance: 8.176e-06 > summary(ans) Formula: y ~ a * b^exp(-c * x) Parameters: Estimate Std. Error t value Pr(>|t|) a 1.051e+02 2.016e+00 52.124 1.63e-15 *** b 1.572e-04 1.633e-04 0.963 0.355 c 3.699e-01 2.170e-02 17.050 8.91e-10 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.355 on 12 degrees of freedom Number of iterations to convergence: 4 Achieved convergence tolerance: 8.176e-06 > plot(x, y, pch=19, col="red") > x2 <- seq(1, 15, by=0.01) > dat2 <- data.frame(x=x2) > y2 <- predict(ans, newdata=dat2) > lines(x2, y2, col="blue") > # eq. 2: y ~ A*exp(-B*C^x) であてはめる(初期値の心配がない) > ans2 <- nls(y ~ SSgompertz(x, A, B, C), data=dat) > p <- ans2$m$getPars() > p[1] # y = a*b^exp(-c*x) の a A 105.0826 > exp(-p[2]) # y = a*b^exp(-c*x) の b B 0.0001571747 > -log(p[3]) # y = a*b^exp(-c*x) の c C 0.3699248
演習問題:
応用問題: