指数・対数を用いた回帰モデル
UPDATE: 2022-12-19 20:46:03
はじめに
ここでは指数関数\(y = \alpha \cdot e^{ \beta \cdot x }\)の対数線型化について内容をまとめておく。ちょいと調べることがあって、そのまとめ。
指数関数の変化量
指数関数\(y = \alpha \cdot e^{ \beta \cdot x }\)は対数線型化すると\(\ln y = ln \alpha + \beta \cdot x\)に変換できる。いわゆる「半対数モデル」に変換できる。左辺にだけ対数が使用されているのでセミログモデルともよばれる。半対数モデルは\(\ln y\)と\(x\)の直線的な関係をみており、\(x\)が1単位増えると、\(\ln y\)は\(\beta\)％だけ変化する。要するに直感的ではないのでわかりにくい・・・私が頭良くないので。この解釈を深ぼってみた。
対数差分
\(\beta\)の意味は\(x\)が1単位増えるときの\(\ln y\)の変化量 ということなので、下記のように差分を計算すると、\(\beta\)というのは変化率の近似値ということがわかる。
Rで実装
指数関数的なインドの人口データをサンプルとして利用する。5年毎の人口数を記録しているデータ。
x <- seq(1950, 2010, length.out = 13)
y <- c(371857, 406374, 447844, 496400, 553874, 622097, 700059,
784491, 873785, 964486, 1053898, 1140043, 1224614)
df <- data.frame(x, y)下記の指数関数でフィットできそうなので、この指数関数を対数化してlm()でパラメタを推定する。もちろん、指数関数を定義してnls()で計算してもよいが、nls()は色々面倒なので、ここではパス。
\[ \begin{eqnarray} y = \alpha \cdot e^{ \beta \cdot x } \end{eqnarray} \]
# y = a*exp(b*x) ⇔ ln(y) ~ ln(a) + b*x
fit <- lm(log(y) ~ x, data = df)
fit##
## Call:
## lm(formula = log(y) ~ x, data = df)
##
## Coefficients:
## (Intercept) x
## -27.61171 0.02073\[ \begin{eqnarray} \hat{ y } = -27.61171 + 0.02073 \cdot x \end{eqnarray} \]
ということなので、-27.61171は\(\ln a\)のことなので、\(a\)を計算するために\(e^{a}\)とした上で、当初の指数関数で予測値を算出する。指数表記すると\(10^{-12} * 1.0195\)だと思うが、ここでは書き下しておく。
\[ \begin{eqnarray} y = 0.000000000001019499 \cdot e^{ 0.02073 \cdot x } \end{eqnarray} \]
# ln(a)=-27.61171 → a = exp(a)
a <- exp(-27.61171)
# a [1] 1.019499e-12
b <- 0.02073
df$fitted <- a * exp(b * x)回帰係数\(b = 0.02073\)の解釈については、
\[ \begin{eqnarray} \frac{ y_{i} - y_{i-1} }{ y_{i-1} } = e^{b \cdot x} -1 \end{eqnarray} \]
ということなので、今回のデータは5年ごとの人口数の推移なので、
\[ \begin{eqnarray} \frac{ y_{i} - y_{i-1} }{ y_{i-1} } = e^{0.02073 \cdot 5} -1 = 0.1107 \end{eqnarray} \]
となり、5年あたりの変化率(人口成長率)は11.07%となり、これが1年であれば、
\[ \begin{eqnarray} \frac{ y_{i} - y_{i-1} }{ y_{i-1} } = e^{0.02073 \cdot 1} - 1 =0.0212 \end{eqnarray} \]
なので、1年あたりの変化率(人口成長率)は2.12%となる。可視化するとこうなる。
library(ggplot2)
ggplot(df)+
geom_point(aes(x, y)) +
geom_line(aes(x, fitted)) +
theme_classic()
nlsパッケージの代わりのminpack.lmパッケージ
こんな感じでコネコネ対数化しなくても、nlsパッケージで指数関数をそのまま渡して、フィッティングしても問題ないと思う・・・が、nlsパッケージはなにかとエラーが出やすいので、私のようなトーシロには扱いにくい。そこでminpack.lmパッケージを使う。minpack.lmパッケージでは、非線形関数に対する最小値問題の解法として、Levenberg-Marquerdt法でパラメタを計算する。Levenberg-Marquerdt法は最初は再急降下法、解に近づいたらガウス・ニュートン法に切り替えて、目的関数を繰り返し計算することで小さくしていく方法。なので、nls.lm()では、目的関数を最小化するように引数を渡すことになる。
library(minpack.lm)
x <- seq(1950, 2010, length.out = 13)
y <- c(371857, 406374, 447844, 496400, 553874, 622097, 700059,
784491, 873785, 964486, 1053898, 1140043, 1224614)
df <- data.frame(x, y)
# 初期パラメタ設定のカンニング
res_lm <- lm(formula = log(y) ~ x, data = df)
init_a <- exp(res_lm$coefficients[1])
init_b <- res_lm$coefficients[2]
x <- df$x
obs <- df$y
pred <- function(p, xx) p$a * exp(p$b * xx)
resid <- function(p, observed, xx) observed - pred(p, xx)
init_par <- list(a = init_a, b = init_b)
fit <- nls.lm(
xx = x,
observed = obs,
fn = resid,
par = init_par,
control = nls.lm.control(maxiter = 1024, nprint = 1)
)## It. 0, RSS = 5.83412e+09, Par. = 1.0195e-12 0.0207347
## It. 1, RSS = 5.69381e+09, Par. = 1.05889e-12 0.0207136
## It. 2, RSS = 5.65148e+09, Par. = 1.13638e-12 0.0206773
## It. 3, RSS = 5.57466e+09, Par. = 1.22425e-12 0.0206398
## It. 4, RSS = 5.47379e+09, Par. = 1.30633e-12 0.0206075
## It. 5, RSS = 5.41075e+09, Par. = 1.3991e-12 0.020573
## It. 6, RSS = 5.33862e+09, Par. = 1.49427e-12 0.0205401
## It. 7, RSS = 5.27435e+09, Par. = 1.5914e-12 0.0205086
## It. 8, RSS = 5.21752e+09, Par. = 1.69028e-12 0.0204784
## It. 9, RSS = 5.16736e+09, Par. = 1.79072e-12 0.0204495
## It. 10, RSS = 5.12316e+09, Par. = 1.89255e-12 0.0204218
## It. 11, RSS = 5.0843e+09, Par. = 1.99562e-12 0.0203953
## It. 12, RSS = 5.05019e+09, Par. = 2.0998e-12 0.0203698
## It. 13, RSS = 5.02034e+09, Par. = 2.20496e-12 0.0203454
## It. 14, RSS = 4.99431e+09, Par. = 2.31102e-12 0.0203218
## It. 15, RSS = 4.9717e+09, Par. = 2.41787e-12 0.0202992
## It. 16, RSS = 4.95215e+09, Par. = 2.52543e-12 0.0202774
## It. 17, RSS = 4.93537e+09, Par. = 2.63364e-12 0.0202564
## It. 18, RSS = 4.92109e+09, Par. = 2.74243e-12 0.0202361
## It. 19, RSS = 4.90905e+09, Par. = 2.85174e-12 0.0202165
## It. 20, RSS = 4.89906e+09, Par. = 2.96152e-12 0.0201976
## It. 21, RSS = 4.89091e+09, Par. = 3.07174e-12 0.0201793
## It. 22, RSS = 4.88444e+09, Par. = 3.18235e-12 0.0201615
## It. 23, RSS = 4.87951e+09, Par. = 3.29332e-12 0.0201444
## It. 24, RSS = 4.87596e+09, Par. = 3.40462e-12 0.0201277
## It. 25, RSS = 4.87368e+09, Par. = 3.51621e-12 0.0201115
## It. 26, RSS = 4.87257e+09, Par. = 3.62808e-12 0.0200958
## It. 27, RSS = 4.87055e+09, Par. = 3.66512e-12 0.020091
## It. 28, RSS = 4.87053e+09, Par. = 3.66455e-12 0.0200911
## It. 29, RSS = 4.87053e+09, Par. = 3.66457e-12 0.0200911summary(fit)##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a.(Intercept) 3.665e-12 3.364e-12 1.089 0.299
## b.x 2.009e-02 4.606e-04 43.615 1.12e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21040 on 11 degrees of freedom
## Number of iterations to termination: 29
## Reason for termination: Relative error in the sum of squares is at most `ftol'.nls_a <- fit$par$a
nls_b <- fit$par$b
df$nls_fitted <- pred(p = list(a = nls_a, b = nls_b), xx = x)
library(ggplot2)
ggplot(df) +
geom_point(aes(x, y)) +
geom_line(aes(x, nls_fitted)) +
theme_classic()
おまけ1
30%が39%になったとき、39/30=1.3なので「9%」の増加ではなく「3%」の増加。これを「9ポイント」の増加と表現する。
おまけ2
両対数モデル対数線型化すると\(\ln y = ln \alpha + \beta \cdot ln x\)は「両対数モデル」。両対数モデルは\(x\)が1％単位増えると、\(\ln y\)は\(\beta\)％だけ変化する。
library(tidyverse)
library(scales)
df <- structure(list(x = c(37.02, 39.51, 47.73, 48.28, 51.27, 56.39,
59.64, 59.71, 60.61, 61.84, 62.14, 62.84, 66.06, 68.67, 69.09,
69.18, 70.08, 71.96, 72.86, 73.48, 74.13, 74.3, 74.67, 75.79,
75.8, 76.34, 76.5, 77.51, 77.61, 79.91, 80.38, 80.39, 80.43,
80.86, 82.51, 86.72, 87.08, 87.36, 89.47, 90.57, 91.97, 92.01,
92.13, 92.48, 94.92, 97.7, 100.76, 101.71, 102.57, 103.2, 104.17,
107.82, 108.04, 109.72, 109.93, 110.8, 111.77, 111.95, 112.29,
113.68, 113.93, 114.16, 115.11, 115.45, 116.13, 116.32, 116.5,
116.93, 119.16, 119.71, 120.03, 120.69, 121.02, 121.67, 122.39,
122.96, 123.59, 124.63, 125.2, 126.13, 128.2, 129.28, 129.71,
130.96, 131.35, 132.23, 132.63, 132.66, 134.41, 135.29, 135.3,
136.01, 136.38, 136.93, 137.24, 138.18, 138.32, 140.77, 144.03,
144.96, 145.15, 149.93, 151.32, 151.67, 152.52, 157.65, 160.26,
160.52, 161.44, 161.58, 161.96, 162.31, 163.3, 164.91, 166.4,
166.43, 166.86, 168.35, 170.67, 170.89, 174.2, 174.6, 174.73,
174.89, 174.93, 176.41, 176.43, 176.53, 177.04, 178.25, 180.31,
181.7, 183.17, 183.3, 183.32, 184.26, 184.35, 185.07, 186.58,
189, 189.17, 189.26, 190.96, 191.45, 192.03, 192.04, 193.1, 193.29,
194.84, 196.03, 196.11, 197.46, 197.75, 198.5, 200.14, 200.21,
201, 202.36, 203.71, 204.43, 205.18, 208.11, 208.93, 209.25,
210.83, 211.69, 212.1, 212.18, 213.36, 214.17, 217.43, 219.82,
220.69, 221.61, 222.18, 224.6, 226.2, 227.49, 227.89, 228.36,
229.47, 230.13, 231.9, 235.72, 236.5, 244.46, 244.51, 244.81,
246.81, 247.07, 249.9, 251.02, 253.3, 254.73, 254.86, 255.57,
255.75, 256.87, 257.86, 258.54), y = c(2222.32, 2202.38, 2590.29,
2281.37, 2370.36, 2348.79, 2463.7, 2511.71, 2385.07, 2339.25,
2555.6, 2547.69, 2438.93, 2546.07, 2785.52, 2665.16, 2494.81,
2689.2, 2603.23, 2410.28, 2467.54, 2420.38, 2549.48, 2737.52,
2539.91, 2349.46, 2803.46, 2636.5, 2528.86, 2617.78, 2788.7,
2564.32, 2825.09, 2623.81, 2637.92, 2727.83, 2910.96, 2753.95,
2724.13, 3048.19, 2674.23, 2685.19, 2709.63, 2653.29, 2683.19,
2819.6, 2908.35, 2773.82, 2844.26, 2645.23, 2873.45, 2685.92,
2927.15, 2671.7, 2585.38, 2594.85, 2983.36, 2992.01, 3179.27,
2894.2, 2807.12, 3058.08, 2836.37, 2850.06, 2909.13, 3015.87,
2888.5, 2873.56, 2975.6, 2883.41, 2852.07, 3064.73, 2835.07,
2859.48, 2977.81, 2850.68, 2948.12, 2847.23, 2892.27, 2828.08,
2897.52, 2977.36, 2914.79, 2856.47, 2990.32, 2966.5, 3111.06,
2606.25, 2781.62, 3001.49, 3145.2, 2821.24, 2742.66, 3097.45,
3000.8, 3115.1, 3058.95, 2796.17, 2927.21, 3153.79, 3033.11,
3192.74, 2916.45, 3104.92, 2927.72, 3022.24, 3004.02, 3137.29,
3096.55, 2965.18, 3157.04, 3285.46, 3068.87, 3104.39, 3199.4,
3114.7, 3214.57, 2849.37, 3363.34, 2847.73, 3074.9, 3018.01,
3095.36, 3060.31, 3103.62, 3131.26, 2949.09, 3278.77, 2969.07,
3233.33, 3073.76, 3126.97, 3236.47, 3355.7, 3089.88, 3085.4,
3309.19, 3149.66, 3278.12, 3084.82, 3315.15, 3166.87, 3038.28,
3114.49, 3151.16, 3242.14, 3211.64, 3269.05, 3262.71, 3070.54,
3085.48, 3289.32, 3260.07, 3170.33, 3214.82, 3314.96, 3196.11,
3147.6, 3284.88, 3132.4, 3156.46, 3165.42, 3194.68, 3480.29,
3310.13, 3274.93, 3190.17, 3362.12, 3334.05, 2993.88, 3315.54,
3259.08, 3351.82, 3544.53, 3422.43, 3394.66, 3352.86, 3200.08,
3205.98, 3269.16, 3280.44, 3239.95, 3355.97, 3463.38, 3249.61,
3405.55, 3410.15, 3322.29, 3319.95, 3466.56, 3244.31, 3126.41,
3348.52, 3367.87, 3206.18, 3475.06, 3296.34, 3499.94, 3411.42,
3273.21)), class = "data.frame", row.names = c(NA, -200L))fit <- lm(log(y) ~ log(x), data = df)
df$p <- exp(fit$fitted.values)
fit##
## Call:
## lm(formula = log(y) ~ log(x), data = df)
##
## Coefficients:
## (Intercept) log(x)
## 6.9418 0.2139係数が0.2139なので、xが1%増加するとyが0.2139％増加するという関係になっている。可視化すると効果が逓減していく様子が見れる。
ggplot() +
geom_point(data = df, aes(x, y), col = 'gray') +
geom_line(data = df, aes(x, p), col = 'royalblue') +
scale_y_continuous(label = comma, breaks = seq(2000, 3500, 100)) +
scale_x_continuous(breaks = seq(00, 300, 10)) +
theme_bw()
このモデルを使って予測値を計算すると、x=100だとy=2771.068、x=200だとy=3214.011となる。
exp(predict(fit, newdata = data.frame(x = 100)))## 1
## 2771.068exp(predict(fit, newdata = data.frame(x = 200)))## 1
## 3214.011下記はメモ。xが1%増加するとyは0.2139%増加するので、どの時点1％増加させても弾力性は常に約0.21%となる。
f <- function(a, b) {
b1 <- fit$coefficients[[1]]
b2 <- fit$coefficients[[2]]
a <- exp(b1 + b2 * log(a))
b <- exp(b1 + b2 * log(b))
diff <- b - a
e <- (diff / a) * 100
return(list(
pre = a,
post = b,
difference = diff,
elasticity = e
))
}# xが1%増加するとyは0.2106%増加するので、どの時点でも弾力性elasticityは常に約0.21%増加する
f(100, 100 * 1.01) # 101## $pre
## [1] 2771.068
##
## $post
## [1] 2776.973
##
## $difference
## [1] 5.905055
##
## $elasticity
## [1] 0.2130967f(150, 150 * 1.01) # 151.5## $pre
## [1] 3022.17
##
## $post
## [1] 3028.61
##
## $difference
## [1] 6.440144
##
## $elasticity
## [1] 0.2130967f(200, 200 * 1.01) # 202## $pre
## [1] 3214.011
##
## $post
## [1] 3220.86
##
## $difference
## [1] 6.848951
##
## $elasticity
## [1] 0.2130967f(400, 400 * 1.01) # 404## $pre
## [1] 3727.756
##
## $post
## [1] 3735.7
##
## $difference
## [1] 7.943724
##
## $elasticity
## [1] 0.2130967f(1000, 1000 * 1.01) # 1010## $pre
## [1] 4535.026
##
## $post
## [1] 4544.69
##
## $difference
## [1] 9.66399
##
## $elasticity
## [1] 0.2130967可視化している通り逓減しているので、実スケールで各時点で1づつ増やすと、増加量はどの時点でも1で同じだが、diffがどんどん小さくなる様子が見れる。1%ではなく1なのでelasticityは21％ではない。
f(100, 101)## $pre
## [1] 2771.068
##
## $post
## [1] 2776.973
##
## $difference
## [1] 5.905055
##
## $elasticity
## [1] 0.2130967f(150, 151)## $pre
## [1] 3022.17
##
## $post
## [1] 3026.469
##
## $difference
## [1] 4.299021
##
## $elasticity
## [1] 0.1422495f(200, 201)## $pre
## [1] 3214.011
##
## $post
## [1] 3217.442
##
## $difference
## [1] 3.431172
##
## $elasticity
## [1] 0.1067567f(400, 401)## $pre
## [1] 3727.756
##
## $post
## [1] 3729.748
##
## $difference
## [1] 1.991765
##
## $elasticity
## [1] 0.05343066f(1000, 1001)## $pre
## [1] 4535.026
##
## $post
## [1] 4535.996
##
## $difference
## [1] 0.9698085
##
## $elasticity
## [1] 0.02138485dput(df)という、オブジェクトを再作成するための関数があるのか。しらなんだ。
dput(df)## structure(list(x = c(37.02, 39.51, 47.73, 48.28, 51.27, 56.39,
## 59.64, 59.71, 60.61, 61.84, 62.14, 62.84, 66.06, 68.67, 69.09,
## 69.18, 70.08, 71.96, 72.86, 73.48, 74.13, 74.3, 74.67, 75.79,
## 75.8, 76.34, 76.5, 77.51, 77.61, 79.91, 80.38, 80.39, 80.43,
## 80.86, 82.51, 86.72, 87.08, 87.36, 89.47, 90.57, 91.97, 92.01,
## 92.13, 92.48, 94.92, 97.7, 100.76, 101.71, 102.57, 103.2, 104.17,
## 107.82, 108.04, 109.72, 109.93, 110.8, 111.77, 111.95, 112.29,
## 113.68, 113.93, 114.16, 115.11, 115.45, 116.13, 116.32, 116.5,
## 116.93, 119.16, 119.71, 120.03, 120.69, 121.02, 121.67, 122.39,
## 122.96, 123.59, 124.63, 125.2, 126.13, 128.2, 129.28, 129.71,
## 130.96, 131.35, 132.23, 132.63, 132.66, 134.41, 135.29, 135.3,
## 136.01, 136.38, 136.93, 137.24, 138.18, 138.32, 140.77, 144.03,
## 144.96, 145.15, 149.93, 151.32, 151.67, 152.52, 157.65, 160.26,
## 160.52, 161.44, 161.58, 161.96, 162.31, 163.3, 164.91, 166.4,
## 166.43, 166.86, 168.35, 170.67, 170.89, 174.2, 174.6, 174.73,
## 174.89, 174.93, 176.41, 176.43, 176.53, 177.04, 178.25, 180.31,
## 181.7, 183.17, 183.3, 183.32, 184.26, 184.35, 185.07, 186.58,
## 189, 189.17, 189.26, 190.96, 191.45, 192.03, 192.04, 193.1, 193.29,
## 194.84, 196.03, 196.11, 197.46, 197.75, 198.5, 200.14, 200.21,
## 201, 202.36, 203.71, 204.43, 205.18, 208.11, 208.93, 209.25,
## 210.83, 211.69, 212.1, 212.18, 213.36, 214.17, 217.43, 219.82,
## 220.69, 221.61, 222.18, 224.6, 226.2, 227.49, 227.89, 228.36,
## 229.47, 230.13, 231.9, 235.72, 236.5, 244.46, 244.51, 244.81,
## 246.81, 247.07, 249.9, 251.02, 253.3, 254.73, 254.86, 255.57,
## 255.75, 256.87, 257.86, 258.54), y = c(2222.32, 2202.38, 2590.29,
## 2281.37, 2370.36, 2348.79, 2463.7, 2511.71, 2385.07, 2339.25,
## 2555.6, 2547.69, 2438.93, 2546.07, 2785.52, 2665.16, 2494.81,
## 2689.2, 2603.23, 2410.28, 2467.54, 2420.38, 2549.48, 2737.52,
## 2539.91, 2349.46, 2803.46, 2636.5, 2528.86, 2617.78, 2788.7,
## 2564.32, 2825.09, 2623.81, 2637.92, 2727.83, 2910.96, 2753.95,
## 2724.13, 3048.19, 2674.23, 2685.19, 2709.63, 2653.29, 2683.19,
## 2819.6, 2908.35, 2773.82, 2844.26, 2645.23, 2873.45, 2685.92,
## 2927.15, 2671.7, 2585.38, 2594.85, 2983.36, 2992.01, 3179.27,
## 2894.2, 2807.12, 3058.08, 2836.37, 2850.06, 2909.13, 3015.87,
## 2888.5, 2873.56, 2975.6, 2883.41, 2852.07, 3064.73, 2835.07,
## 2859.48, 2977.81, 2850.68, 2948.12, 2847.23, 2892.27, 2828.08,
## 2897.52, 2977.36, 2914.79, 2856.47, 2990.32, 2966.5, 3111.06,
## 2606.25, 2781.62, 3001.49, 3145.2, 2821.24, 2742.66, 3097.45,
## 3000.8, 3115.1, 3058.95, 2796.17, 2927.21, 3153.79, 3033.11,
## 3192.74, 2916.45, 3104.92, 2927.72, 3022.24, 3004.02, 3137.29,
## 3096.55, 2965.18, 3157.04, 3285.46, 3068.87, 3104.39, 3199.4,
## 3114.7, 3214.57, 2849.37, 3363.34, 2847.73, 3074.9, 3018.01,
## 3095.36, 3060.31, 3103.62, 3131.26, 2949.09, 3278.77, 2969.07,
## 3233.33, 3073.76, 3126.97, 3236.47, 3355.7, 3089.88, 3085.4,
## 3309.19, 3149.66, 3278.12, 3084.82, 3315.15, 3166.87, 3038.28,
## 3114.49, 3151.16, 3242.14, 3211.64, 3269.05, 3262.71, 3070.54,
## 3085.48, 3289.32, 3260.07, 3170.33, 3214.82, 3314.96, 3196.11,
## 3147.6, 3284.88, 3132.4, 3156.46, 3165.42, 3194.68, 3480.29,
## 3310.13, 3274.93, 3190.17, 3362.12, 3334.05, 2993.88, 3315.54,
## 3259.08, 3351.82, 3544.53, 3422.43, 3394.66, 3352.86, 3200.08,
## 3205.98, 3269.16, 3280.44, 3239.95, 3355.97, 3463.38, 3249.61,
## 3405.55, 3410.15, 3322.29, 3319.95, 3466.56, 3244.31, 3126.41,
## 3348.52, 3367.87, 3206.18, 3475.06, 3296.34, 3499.94, 3411.42,
## 3273.21), p = c(2240.38037034626, 2271.79826319243, 2365.53995494588,
## 2371.34519847714, 2402.02528404025, 2451.44008857735, 2481.00390389683,
## 2481.62658321071, 2489.58178530971, 2500.30506895369, 2502.89504137596,
## 2508.90030415978, 2535.86571510063, 2556.97458284405, 2560.31225777172,
## 2561.02539785959, 2568.11698556226, 2582.70259431961, 2589.57926647347,
## 2594.27778528765, 2599.17030839184, 2600.44432593224, 2603.20928669568,
## 2611.5137716685, 2611.58748307999, 2615.55659544518, 2616.72838960812,
## 2624.08120412066, 2624.80510026217, 2641.25579707851, 2644.57155201724,
## 2644.64193423058, 2644.92339428163, 2647.94216197064, 2659.40997133342,
## 2687.87426165191, 2690.25746769546, 2692.10572377664, 2705.88591424077,
## 2712.96884805259, 2721.88634413486, 2722.13955736939, 2722.89867824478,
## 2725.10835024488, 2740.33294454458, 2757.30855673834, 2775.56047472046,
## 2781.13823246295, 2786.15236461504, 2789.80457706815, 2795.3937076951,
## 2816.06516360143, 2817.29343654178, 2826.60871210033, 2827.76522258412,
## 2832.53805924224, 2837.82489020986, 2838.80198110709, 2840.64423284819,
## 2848.13049123159, 2849.46929711287, 2850.69896049768, 2855.75746436268,
## 2857.55990279063, 2861.15229817574, 2862.15309956176, 2863.1000425472,
## 2865.35753751298, 2876.96143125546, 2879.79710760388, 2881.44224614497,
## 2884.82448568095, 2886.51015611952, 2889.8198764639, 2893.46984294128,
## 2896.34744829406, 2899.51578557881, 2904.71838937826, 2907.55535375076,
## 2912.16236490308, 2922.32160422183, 2927.57097696193, 2929.65141310157,
## 2935.66856093716, 2937.53666935151, 2941.73592378911, 2943.63741702189,
## 2943.77984720367, 2952.04479788754, 2956.16896477605, 2956.21570898072,
## 2959.52762755909, 2961.24817348758, 2963.79897327419, 2965.23314968994,
## 2969.56643097173, 2970.20982919384, 2981.38728119133, 2996.0254176882,
## 3000.15354407544, 3000.99436231981, 3021.86828973059, 3027.84004439716,
## 3029.33692455211, 3032.96093279544, 3054.50206551599, 3065.25077554119,
## 3066.31397392861, 3070.06522437753, 3070.63459276612, 3072.17806887516,
## 3073.59717522137, 3077.59823551478, 3084.06448039375, 3090.00470469955,
## 3090.12387640007, 3091.8301513232, 3097.71598073267, 3106.7994937185,
## 3107.65581283041, 3120.4360871494, 3121.96757085315, 3122.46470922133,
## 3123.07617285583, 3123.22897005338, 3128.86325465055, 3128.93913876362,
## 3129.31845793251, 3131.25036256436, 3135.81645972724, 3143.53440199803,
## 3148.70306489989, 3154.13549864914, 3154.61426660795, 3154.68789952485,
## 3158.14154516067, 3158.47148633165, 3161.10646757032, 3166.60654367825,
## 3175.34870167677, 3175.95950571351, 3176.28269789289, 3182.36486452688,
## 3184.11005540226, 3186.17125968937, 3186.20675475471, 3189.96101916427,
## 3190.63224034402, 3196.08870648765, 3200.25475683069, 3200.53411359082,
## 3205.23480248856, 3206.24128216261, 3208.83887695685, 3214.49218242884,
## 3214.73267056752, 3217.44217929448, 3222.08710324257, 3226.67366759006,
## 3229.11007310954, 3231.64083337856, 3241.45851113718, 3244.18665139009,
## 3245.2490085926, 3250.47575567431, 3253.30776490436, 3254.65472632634,
## 3254.9173093856, 3258.78139892482, 3261.42415477007, 3271.98164075601,
## 3279.6428626893, 3282.41542259136, 3285.3379941895, 3287.14393622623,
## 3294.77095459092, 3299.77821018781, 3303.79507723186, 3305.03698028773,
## 3306.49402876384, 3309.92580872184, 3311.96014494403, 3317.39331286529,
## 3329.00894799259, 3331.36251132133, 3355.0386685954, 3355.18546034688,
## 3356.065715685, 3361.912514717, 3362.66986030468, 3370.87302077027,
## 3374.0993401029, 3380.63239668374, 3384.70633223007, 3385.07579781532,
## 3387.09103843239, 3387.60124512419, 3390.76953779014, 3393.56105340668,
## 3395.47358129925)), row.names = c(NA, -200L), class = "data.frame")