« 2023年10月 | トップページ | 2023年12月 »
その環境下に搭載されているタイムゾーンを識別するID（タイムゾーンID）を一覧形式で得るには、OlsonNames関数を使う。
> OlsonNames()
[1] "Africa/Abidjan" "Africa/Accra"
[3] "Africa/Addis_Ababa" "Africa/Algiers"
[5] "Africa/Asmara" "Africa/Asmera"
（以下、表示省略）
> grep("Tokyo", OlsonNames())
[1] 322
> grep("Tokyo", OlsonNames(), value = TRUE)
[1] "Asia/Tokyo"
Sys.timezone関数を使う。以下は、Windowsの日本語版バージョン4.3.1で実行した例。
> Sys.timezone()
[1] "Asia/Tokyo"
以下の表5.4（p.232）の値をCSV形式で入力したtable5_4.csvを、カレントディレクトリに置いておく。
i, sex, xxi, yyi
1, 女性, 255.5, 15.8
2, 女性, 259.7, 16.8
3, 女性, 385.3, 16.4
4, 女性, 286.5, 15.7
5, 女性, 302.1, 16.3
6, 女性, 317.8, 17.4
7, 女性, 321.5, 18.3
8, 女性, 393.0, 17.6
9, 女性, 393.7, 17.4
10, 女性, 289.2, 15.3
11, 男性, 352.2, 9.7
12, 男性, 367.1, 5.2
13, 男性, 276.9, 7.1
14, 男性, 348.4, 9.6
15, 男性, 277.8, 6.9
16, 男性, 386.5, 8.2
17, 男性, 337.7, 7.3
18, 男性, 219.2, 6.6
19, 男性, 330.2, 8.7
20, 男性, 395.8, 10.1
計算する。
> dtf <- read.csv("table5_4.csv")
> r <- lm(yyi ~ 1 + xxi, dtf, trimws(sex) == "女性")
> print(r$coefficients) # 女性だけ
(Intercept) xxi
13.722088020 0.009293487
> r <- lm(yyi ~ 1 + xxi, dtf, trimws(sex) == "男性")
> print(r$coefficients) # 男性だけ
(Intercept) xxi
3.5590101 0.0133088
> ddi <- ifelse(trimws(dtf$sex) == "女性", 1, 0)
> r <- lm(dtf$yyi ~ 1 + dtf$xxi + ddi)
> print(summary(r)) # 女性と男性
Call:
lm(formula = dtf$yyi ~ 1 + dtf$xxi + ddi)
Residuals:
Min 1Q Median 3Q Max
-3.1720 -0.4944 -0.1475 0.7592 1.5878
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.19004 1.74965 2.395 0.0284 *
dtf$xxi 0.01139 0.00519 2.195 0.0423 *
ddi 8.85968 0.53539 16.548 6.45e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.193 on 17 degrees of freedom
Multiple R-squared: 0.9417, Adjusted R-squared: 0.9348
F-statistic: 137.2 on 2 and 17 DF, p-value: 3.234e-11
> yyi <- c(3, 6, 8, 12, 11)
> xx2i <- c(2, 4, 6, 8, 10)
> xx3i <- c(1, 1, 2, 2, 4)
> n <- length(yyi)
> kk <- 3
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- as.matrix(data.frame(rep(1, n), xx2i, xx3i))
> mxxxtxxi <- solve(t(mxxx) %*% mxxx)
> mxbh <- mxxxtxxi %*% t(mxxx) %*% mxy
> bhk <- as.vector(mxbh)
> mxyh <- mxxx %*% mxbh
> ui <- as.vector(mxy - mxyh)
> sh2 <- sum(ui ^ 2) / (n - kk)
> sk2 <- as.vector(diag(sh2 * mxxxtxxi))
> sk <- sqrt(sk2)
> tk <- bhk / sk
> print(n - kk) # 自由度
[1] 2
> print(sh2) # 残差分散
[1] 0.1818182
> print(sk) # 推定した回帰モデルのパラメーターの標準誤差
[1] 0.4490578 0.1574592 0.4065578
> print(tk) # t値
[1] 2.631773 11.835681 -5.366563
以下の表4.5（p.186）の値をCSV形式で入力したtable4_5.csvを、カレントディレクトリに置いておく。
i, xx2i, xx3i, yyi
1, 2, 26, 3
2, 4, 5, 5
3, 6, 21, 6
4, 8, 8, 10
5, 10, 20, 11
計算する。
> dtf <- read.csv("table4_5.csv", header = TRUE)
> xx2i <- dtf$xx2i
> xx3i <- dtf$xx3i
> yyi <- dtf$yyi
> kk <- 3
> n <- nrow(dtf)
> degf <- n - kk
> qv <- 0.05
> tv <- qt(1 - 0.05 / 2, 2)
> r <- lm(yyi ~ 1 + xx2i + xx3i)
> bhk <- as.vector(r$coefficient)
> sh2 <- sum(r$residuals ^ 2) / degf
> sk2 <- sk <- tk <- double(kk)
> xx2m <- mean(xx2i)
> xx3m <- mean(xx3i)
> ss22 <- sum((xx2i - xx2m) * (xx2i - xx2m))
> ss23 <- sum((xx2i - xx2m) * (xx3i - xx3m))
> ss33 <- sum((xx3i - xx3m) * (xx3i - xx3m))
> d1 <- xx2m ^ 2 * ss33 - 2 * xx2m * xx3m * ss23 + xx3m ^ 2 * ss22
> d2 <- ss22 * ss33 - ss23 ^ 2
> sk2[1] <- sh2 * (1 / n + d1 / d2)
> sk2[2] <- sh2 * ss33 / d2
> sk2[3] <- sh2 * ss22 / d2
> sk <- sqrt(sk2)
> tk <- bhk / sk
> # 残差分散 σ^^2
> print(sh2)
[1] 0.7197232
> # 推定値 β^^k の分散
> print(sk2)
[1] 1.583391003 0.018451538 0.002263992
> # 推定値 β^^k の標準誤差
> print(sk)
[1] 1.25832865 0.13583644 0.04758143
> # 推定値 β^^k のt値
> print(tk)
[1] 1.1219364 7.6037916 -0.7999399
stringrパッケージのstr_replace_all関数を使う。元となる文字列にはベクトルを指定することができる。str_replaceという別の関数もあり、こちらは各文字列で最初にパターンにマッチした箇所しか置換をしない。
> library("stringr")
> s <- c("ABCD", "BCDE", "CDEF")
> str_replace_all(s, "B", "b")
[1] "AbCD" "bCDE" "CDEF"
> str_replace_all(s, "^B", "b")
[1] "ABCD" "bCDE" "CDEF"
> str_replace_all(s, "[E-F]$", "b")
[1] "ABCD" "BCDb" "CDEb"
> str_replace_all(s, "[E-F]", "b")
[1] "ABCD" "BCDb" "CDbb"
> str_replace(s, "[E-F]", "b")
[1] "ABCD" "BCDb" "CDbF"
ただ文字列を置換するのではなく、パターンに含まれる文字列を使用して置き換えたい場合は、キャプチャによるグループ化を利用する。パターン内で（ ）（括弧）に挟まれた範囲はグループ化（キャプチャという）される。このグループには、キャプチャされた順番に1から始まる連番が付けられ、置換後の文字列で\1, \2, …と挿入することで、置換後の文字列に挿入することができる（実際にはエスケープシーケンスを使用する必要があるので、\\1, \\2, …とする）。
以下は7桁の数字（○○○○○○○）を郵便番号の表記（○○○－○○○○）に置き換えた例。
> s <- c("1234567", "9876543", "3210")
> str_replace_all(s, "([0-9]{3})([0-9]{4})", "\\1-\\2")
[1] "123-4567" "987-6543" "3210"
gsub関数を使う。元となる文字列にはベクトルを指定することができる。
> s <- c("ABCD", "BCDE", "CDEF")
> gsub("B", "b", s)
[1] "AbCD" "bCDE" "CDEF"
> gsub("^B", "b", s)
[1] "ABCD" "bCDE" "CDEF"
> gsub("[E-F]", "b", s)
[1] "ABCD" "BCDb" "CDbb"
> gsub("[E-F]$", "b", s)
[1] "ABCD" "BCDb" "CDEb"
ただ文字列を置換するのではなく、パターンに含まれる文字列を使用して置き換えたい場合は、キャプチャによるグループ化を利用する。パターン内で（ ）（括弧）に挟まれた範囲はグループ化（キャプチャという）される。このグループには、キャプチャされた順番に1から始まる連番が付けられ、置換後の文字列で\1, \2, …と挿入することで、置換後の文字列に挿入することができる（実際にはエスケープシーケンスを使用する必要があるので、\\1, \\2, …とする）。
以下は7桁の数字（○○○○○○○）を郵便番号の表記（○○○－○○○○）に置き換えた例。
> s <- c("1234567", "9876543", "3210")
> gsub("([0-9]{3})([0-9]{4})", "\\1-\\2", s)
[1] "123-4567" "987-6543" "3210"
condaコマンドを実行して、すべてのパッケージをアップデートしてみる。
>conda update --all
Collecting package metadata (current_repodata.json): failed
CondaHTTPError: HTTP 000 CONNECTION FAILED for url <https://repo.anaconda.com/pkgs/main/win-64/current_repodata.json>
Elapsed: -
An HTTP error occurred when trying to retrieve this URL.
HTTP errors are often intermittent, and a simple retry will get you on your way.
If your current network has https://www.anaconda.com blocked, please file
a support request with your network engineering team.
'https//repo.anaconda.com/pkgs/main/win-64'
上記のように表示された場合は、接続に失敗している。会社などの環境であれば、プロキシの設定を行うことで解決できる場合があるため、プロキシの設定を行う。
フォルダー「C:\Users\○○○\anaconda3」（○○○はアカウント名）にテキストファイル.condarcを作成（既にある場合はそれを使用）し、以下の3行を追加する。
proxy_servers:
http: http://○○○:×××@proxy.company.co.jp:8080
https: https://○○○:×××@proxy.company.co.jp:8080
それぞれ以下のとおり。なお、「http:」や「https:」の前は必ずスペース（0x20）でインデントすること。タブ文字や全角空白（0x8140）は絶対に使用してはいけない。
○○○： プロキシに接続するためのアカウント
×××： プロキシに接続するためのパスワード
proxy.company.co.jp： プロキシのアドレス
8080： プロキシのポート
設定が反映されて確認をする。以下のコマンドを実行する。
>conda config --show proxy_servers
proxy_servers:
http: http://○○○:×××@proxy.company.co.jp:8080
https: https://○○○:×××@proxy.company.co.jp:8080
設定したとおりに表示されればよい。行頭の空白の数は同じでなくてもかまわない。
再度、condaですべてのパッケージをアップデートしてみる。
>conda update --all
Collecting package metadata (current_repodata.json): done
Solving environment: done
（表示省略）
Proceed ([y]/n)?
最後に「Proceed ([y]/n)?」と表示されてキー入力の状態になれば、接続は成功。ここで「y」を選択すると全パッケージのアップデートが始まる。取りやめる場合は「n」を選択する。
なお、正しく設定しているのにどうしてもうまくいかない場合は、以下のように「http:」も「https:」も「http://～」と設定してみること。
proxy_servers:
http: http://○○○:×××@proxy.company.co.jp:8080
https: http://○○○:×××@proxy.company.co.jp:8080
lm関数を使って計算している。以下に表示している表3.4（p.116）の値をCSV形式で入力したtable3_4.csvをカレントディレクトリに置いておくこと。
i, W, YEAR, DX
1, 209.0, 0.3, 1
2, 237.0, 2.0, 1
3, 301.7, 5.9, 1
4, 371.6, 10.1, 1
5, 455.3, 14.9, 1
6, 526.6, 20.9, 1
7, 579.9, 25.4, 1
8, 587.3, 29.3, 1
9, 477.0, 33.8, 1
10, 201.7, 0.3, 0
11, 228.7, 2.6, 0
12, 271.5, 6.3, 0
13, 308.1, 9.8, 0
14, 342.7, 14.5, 0
15, 389.0, 20.0, 0
16, 394.6, 23.8, 0
17, 399.8, 28.3, 0
18, 317.0, 32.6, 0
実行してみる。
> dtf <- read.csv("table3_4.csv", header = TRUE)
> xx1i <- as.double(dtf$YEAR)
> xx2i <- as.double(dtf$DX)
> yyi <- as.double(dtf$W)
> r <- lm(yyi ~ 1 + log(xx1i) + I(xx1i ^ 2) + I(xx1i ^ 4) + I(xx2i * xx1i))
> print(summary(r))
Call:
lm(formula = yyi ~ 1 + log(xx1i) + I(xx1i^2) + I(xx1i^4) + I(xx2i *
xx1i))
Residuals:
Min 1Q Median 3Q Max
-17.136 -7.683 2.097 7.070 17.824
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.196e+02 4.748e+00 46.26 8.20e-16 ***
log(xx1i) 1.469e+01 3.257e+00 4.51 0.000586 ***
I(xx1i^2) 5.088e-01 3.531e-02 14.41 2.27e-09 ***
I(xx1i^4) -4.446e-04 2.817e-05 -15.78 7.36e-10 ***
I(xx2i * xx1i) 6.523e+00 2.722e-01 23.96 3.84e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11.01 on 13 degrees of freedom
Multiple R-squared: 0.9937, Adjusted R-squared: 0.9917
F-statistic: 511.7 on 4 and 13 DF, p-value: 3.746e-14
> yyhi <- r$fitted.value
> ei <- yyi - yyhi
> ri <- ei / yyi * 100
> ai2 <- ei ^ 2 / sum(ei ^ 2) * 100
> dtf <- data.frame(yyi, xx1i, yyhi, ei, ri, ai2)
> names(dtf) <- c("W", "YEAR", "Wの推定値", "残差", "誤差率", "平方残差率")
> print(dtf)
W YEAR Wの推定値 残差 誤差率 平方残差率
1 209.0 0.3 203.9620 5.0379825 2.4105179 1.609669358
2 237.0 2.0 244.9042 -7.9041701 -3.3350929 3.962195699
3 301.7 5.9 301.3806 0.3194241 0.1058747 0.006470811
4 371.6 10.1 366.7761 4.8238794 1.2981376 1.475761750
5 455.3 14.9 447.5628 7.7372332 1.6993704 3.796599095
6 526.6 20.9 538.0347 -11.4347032 -2.1714210 8.292264754
7 579.9 25.4 576.0246 3.8753715 0.6682827 0.952467251
8 587.3 29.3 569.4762 17.8237961 3.0348708 20.147621237
9 477.0 33.8 492.7775 -15.7775142 -3.3076550 15.787033282
10 201.7 0.3 202.0051 -0.3050563 -0.1512426 0.005901786
11 228.7 2.6 237.1031 -8.4031158 -3.6742964 4.478205967
12 271.5 6.3 266.1784 5.3215546 1.9600569 1.795975490
13 308.1 9.8 297.9372 10.1627947 3.2985377 6.550126198
14 342.7 14.5 346.2444 -3.5444215 -1.0342636 0.796735341
15 389.0 20.0 396.0212 -7.0212080 -1.8049378 3.126418007
16 394.6 23.8 411.7361 -17.1361000 -4.3426508 18.622901255
17 399.8 28.3 391.0281 8.7718944 2.1940706 4.879891365
18 317.0 32.6 309.3476 7.6523586 2.4139932 3.713761354
5つのモデルについて、lm関数を使用して回帰パラメーターを推定したところまで。
> dtf <- read.csv("table2_5.csv", header = TRUE)
> xxi <- as.double(dtf$AGE)
> yyi <- as.double(dtf$BP)
> # (AGE, BP)
> r <- lm(yyi ~ 1 + xxi)
> print(r)
Call:
lm(formula = yyi ~ 1 + xxi)
Coefficients:
(Intercept) xxi
112.3167 0.4451
> # (AGE, log(BP))
> r <- lm(log(yyi) ~ 1 + xxi)
> print(r)
Call:
lm(formula = log(yyi) ~ 1 + xxi)
Coefficients:
(Intercept) xxi
4.732231 0.003365
> # (log(AGE), log(BP))
> r <- lm(log(yyi) ~ 1 + log(xxi))
> print(r)
Call:
lm(formula = log(yyi) ~ 1 + log(xxi))
Coefficients:
(Intercept) log(xxi)
4.411 0.127
> # (AGE^2 / 10000, log(BP))
> r <- lm(log(yyi) ~ 1 + I(xxi ^ 2 / 10000))
> print(r)
Call:
lm(formula = log(yyi) ~ 1 + I(xxi^2/10000))
Coefficients:
(Intercept) I(xxi^2/10000)
4.7942 0.3863
> # (AGE^1.825 / 10000, log(BP))
> r <- lm(log(yyi) ~ 1 + I(xxi ^ 1.825 / 10000))
> print(r)
Call:
lm(formula = log(yyi) ~ 1 + I(xxi^1.825/10000))
Coefficients:
(Intercept) I(xxi^1.825/10000)
4.7883 0.8272
> dtf <- read.csv("table2_1.csv", header = TRUE)
> n <- nrow(dtf)
> xxi <- 1 / dtf$WV
> yyi <- dtf$DC
> k <- 2
> degf <- n - k
> xxm <- mean(xxi)
> yym <- mean(yyi)
> ssx2 <- sum((xxi - xxm) ^ 2)
> ssxy <- sum((xxi - xxm) * (yyi - yym))
> bh <- ssxy / ssx2
> ah <- mean(yyi) - bh * xxm
> sse2 <- sum((yyi - (ah + bh * xxi)) ^ 2)
> s <- sqrt(sse2 / (n - k))
> xx0i <- c(xxi, 1 / 20)
> yy0i <- c(yyi, 2.632)
> yyh0i <- ah + bh * xx0i
> #
> tl2 <- qt(0.025, degf, lower.tail = FALSE)
> yycon <- tl2 * s * sqrt(1 / n + (xx0i - xxm) ^ 2 / ssx2)
> yypre <- tl2 * s * sqrt(1 + 1 / n + (xx0i - xxm) ^ 2 / ssx2)
> dtf <- data.frame(yy0i, yyh0i, yyh0i - yycon, yyh0i + yycon, yyh0i - yypre, yyh0i + yypre)
> names(dtf) <- c("DC", "DCの予測値", "E(DC)予測下限", "E(DC)予測上限", "DC予測下限", "DC予測上限")
> print(dtf)
DC DCの予測値 E(DC)予測下限 E(DC)予測上限 DC予測下限 DC予測上限
1 1.582 1.5919507 1.55297389 1.6309275 1.39328146 1.7906199
2 1.822 1.8231023 1.78198202 1.8642225 1.62400142 2.0222031
3 1.057 0.9392874 0.88252512 0.9960497 0.73637800 1.1421968
4 0.500 0.4105093 0.32701897 0.4939996 0.19856376 0.6224548
5 2.236 2.2854054 2.22839667 2.3424142 2.08242693 2.4883839
6 2.386 2.2639584 2.20790651 2.3200103 2.06124655 2.4666702
7 2.294 2.2527296 2.19717273 2.3082864 2.05015405 2.4553051
8 0.558 0.7052381 0.63727035 0.7732058 0.49891337 0.9115628
9 2.166 2.1279955 2.07762556 2.1783654 1.92678065 2.3292103
10 1.866 1.8603848 1.81847394 1.9022956 1.66111915 2.0596504
11 0.653 0.5876369 0.51361863 0.6616553 0.37924072 0.7960332
12 1.930 1.8868055 1.84426818 1.9293428 1.68740713 2.0862038
13 1.562 1.4713499 1.43146889 1.5112309 1.27250127 1.6701985
14 1.737 1.7832486 1.74284605 1.8236511 1.58429470 1.9822024
15 2.088 2.0417592 1.99457587 2.0889425 1.84131831 2.2422000
16 1.137 1.0525970 1.00068770 1.1045063 0.85099133 1.2542027
17 2.179 2.0954783 2.04635313 2.1446036 1.89457150 2.2963852
18 2.112 2.1908434 2.13793584 2.2437510 1.98897840 2.3927085
19 1.800 1.9882106 1.94280548 2.0336156 1.78818081 2.1882403
20 1.501 1.7064662 1.66705050 1.7458818 1.50771036 1.9052220
21 2.303 2.2168220 2.16281924 2.2708248 2.01466717 2.4189768
22 2.310 2.2990026 2.24137972 2.3566255 2.09585075 2.5021544
23 1.194 1.2875072 1.24378717 1.3312272 1.08785318 1.4871611
24 1.144 1.2232786 1.17762784 1.2689293 1.02319292 1.4233642
25 0.123 0.1484327 0.05037867 0.2464867 -0.06966102 0.3665264
26 2.632 2.6321328 2.55808466 2.7061810 2.42372598 2.8405396
> # 図2.13
> idx <- order(xx0i)
> png("fig2_13.png", width = 512, height = 512)
> plot(xx0i, yy0i, type = "n")
> lines(xx0i[idx], yyh0i[idx] - yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] + yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] - yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx] + yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx], lty = "solid", col = "gray", lwd = 4)
> points(xx0i, yy0i, pch = 20)
> dev.off()
null device
1
> #
> tl2 <- qt(0.005, degf, lower.tail = FALSE)
> yycon <- tl2 * s * sqrt(1 / n + (xx0i - xxm) ^ 2 / ssx2)
> yypre <- tl2 * s * sqrt(1 + 1 / n + (xx0i - xxm) ^ 2 / ssx2)
> dtf <- data.frame(yy0i, yyh0i, yyh0i - yycon, yyh0i + yycon, yyh0i - yypre, yyh0i + yypre)
> names(dtf) <- c("DC", "DCの予測値", "E(DC)予測下限", "E(DC)予測上限", "DC予測下限", "DC予測上限")
> print(dtf)
DC DCの予測値 E(DC)予測下限 E(DC)予測上限 DC予測下限 DC予測上限
1 1.582 1.5919507 1.53905601 1.6448454 1.3223405 1.8615609
2 1.822 1.8231023 1.76729876 1.8789058 1.5529063 2.0932982
3 1.057 0.9392874 0.86225639 1.0163185 0.6639229 1.2146519
4 0.500 0.4105093 0.29720616 0.5238124 0.1228821 0.6981365
5 2.236 2.2854054 2.20803993 2.3627709 2.0099472 2.5608637
6 2.386 2.2639584 2.18789146 2.3400253 1.9888620 2.5390547
7 2.294 2.2527296 2.17733445 2.3281247 1.9778182 2.5276409
8 0.558 0.7052381 0.61300037 0.7974758 0.4252388 0.9852374
9 2.166 2.1279955 2.05963943 2.1963515 1.8549307 2.4010603
10 1.866 1.8603848 1.80350838 1.9172612 1.5899652 2.1308044
11 0.653 0.5876369 0.48718811 0.6880858 0.3048264 0.8704475
12 1.930 1.8868055 1.82907892 1.9445320 1.6162058 2.1574051
13 1.562 1.4713499 1.41722815 1.5254716 1.2014962 1.7412035
14 1.737 1.7832486 1.72841908 1.8380780 1.5132521 2.0532450
15 2.088 2.0417592 1.97772761 2.1057908 1.7697447 2.3137736
16 1.137 1.0525970 0.98215187 1.1230422 0.7790018 1.3261922
17 2.179 2.0954783 2.02881146 2.1621452 1.8228315 2.3681252
18 2.112 2.1908434 2.11904355 2.2626433 1.9168963 2.4647906
19 1.800 1.9882106 1.92659220 2.0498289 1.7167540 2.2596671
20 1.501 1.7064662 1.65297592 1.7599564 1.4367385 1.9761939
21 2.303 2.2168220 2.14353588 2.2901081 1.9424816 2.4911625
22 2.310 2.2990026 2.22080369 2.3772015 2.0233091 2.5746961
23 1.194 1.2875072 1.22817560 1.3468387 1.0165606 1.5584538
24 1.144 1.2232786 1.16132684 1.2852303 0.9517462 1.4948110
25 0.123 0.1484327 0.01536546 0.2814999 -0.1475381 0.4444035
26 2.632 2.6321328 2.53164348 2.7326221 2.3493079 2.9149577
> # 図2.14
> png("fig2_14.png", width = 512, height = 512)
> plot(xx0i, yy0i, type = "n")
> lines(xx0i[idx], yyh0i[idx] - yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] + yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] - yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx] + yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx], lty = "solid", col = "gray", lwd = 4)
> points(xx0i, yy0i, pch = 20)
> dev.off()
null device
1
pp.64-67に掲載されている、p.23の式(2.30)による結果を流れに沿って計算している。表2.1（p.62）の値をCSV形式で入力したtable2_1.csvをカレントディレクトリに置いておくこと。
> dtf <- read.csv("table2_1.csv", header = TRUE)
> xxi <- as.double(1 /dtf$WV)
> yyi <- as.double(dtf$DC)
> n <- nrow(dtf)
> k <- 2
> ssxx <- sum(xxi)
> ssyy <- sum(yyi)
> ssxx2 <- sum(xxi ^ 2)
> ssyy2 <- sum(yyi ^ 2)
> ssxxyy <- sum(xxi * yyi)
> xxm <- mean(xxi)
> yym <- mean(yyi)
> ssx2 <- ssxx2 - ssxx ^ 2 / n
> ssy2 <- ssyy2 - ssyy ^ 2 / n
> ssxy <- ssxxyy - ssxx * ssyy / n
> bh <- ssxy / ssx2
> ah <- yym - bh * xxm
> ssyh2 <- bh * ssxy
> r2 <- ssyh2 / ssy2
> sse2 <- ssy2 - ssyh2
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> sah2 <- s2 * (1 / n + xxm ^ 2 / ssx2)
> sah <- sqrt(sah2)
> ta <- ah / sah
> sbh2 <- s2 / ssx2
> sbh <- sqrt(sbh2)
> tb <- bh / sbh
> yyhi <- ah + bh * xxi
> ei <- yyi - yyhi
> ri <- 100 * ei / yyi
> ai2 <- 100 * (ei ^ 2 / sse2)
> sh2 <- sum(ei ^ 2) / n
> print(ah) # 推定されたパラメーター α^
[1] 2.97886
> print(bh) # 推定されたパラメーター β^
[1] -6.934547
> print(sh2) # 誤差項の分散の最尤推定量 σ^^2
[1] 0.008158786
> dtf <- data.frame(yyi, yyhi, ei, ri, ai2)
> colnames(dtf) <- c("DC", "DC推定値", "残差", "誤差率(%)", "平方残差率(%)")
> print(dtf)
DC DC推定値 残差 誤差率(%) 平方残差率(%)
1 1.582 1.5919507 -0.009950703 -0.62899514 0.048544719
2 1.822 1.8231023 -0.001102282 -0.06049844 0.000595689
3 1.057 0.9392874 0.117712578 11.13647849 6.793290635
4 0.500 0.4105093 0.089490699 17.89813971 3.926361222
5 2.236 2.2854054 -0.049405439 -2.20954556 1.196696381
6 2.386 2.2639584 0.122041615 5.11490423 7.302143226
7 2.294 2.2527296 0.041270439 1.79906010 0.835050285
8 0.558 0.7052381 -0.147238090 -26.38675451 10.628569260
9 2.166 2.1279955 0.038004532 1.75459519 0.708117342
10 1.866 1.8603848 0.005615206 0.30092206 0.015458445
11 0.653 0.5876369 0.065363052 10.00965576 2.094590369
12 1.930 1.8868055 0.043194527 2.23805841 0.914727865
13 1.562 1.4713499 0.090650121 5.80346482 4.028758432
14 1.737 1.7832486 -0.046248561 -2.66255390 1.048650839
15 2.088 2.0417592 0.046240820 2.21459864 1.048299802
16 1.137 1.0525970 0.084402980 7.42330521 3.492609407
17 2.179 2.0954783 0.083521654 3.83302682 3.420051394
18 2.112 2.1908434 -0.078843429 -3.73311692 3.047652653
19 1.800 1.9882106 -0.188210552 -10.45614178 17.366903572
20 1.501 1.7064662 -0.205466164 -13.68861852 20.697366094
21 2.303 2.2168220 0.086177997 3.74198856 3.641055032
22 2.310 2.2990026 0.010997410 0.47607833 0.059294615
23 1.194 1.2875072 -0.093507161 -7.83142048 4.286710893
24 1.144 1.2232786 -0.079278565 -6.92994450 3.081385381
25 0.123 0.1484327 -0.025432682 -20.67697705 0.317116449
> # 図2.8と2.9の作成
> png("fig2_8.png", width = 512, height = 512)
> plot(1 / xxi, yyi, xlim = c(2, 11), ylim = c(0, 2.5), xlab = "WV", ylab = "DC", pch = 20)
> idx <- order(1 / xxi)
> lines(1 / xxi[idx], yyhi[idx])
> dev.off()
null device
1
> png("fig2_9.png", width = 512, height = 512)
> plot(xxi, yyi, xlim = c(0.09, 0.42), ylim = c(0, 2.5), xlab = "1 / WV", ylab = "DC", pch = 20)
> idx <- order(xxi)
> lines(xxi[idx], yyhi[idx])
> dev.off()
null device
1
p.63に掲載されている3つのモデルについて、それぞれlm関数を使って計算している。表2.1（p.62）の値をCSV形式で入力したtable2_1.csvをカレントディレクトリに置いておくこと。
> dtf <- read.csv("table2_1.csv", header = TRUE)
> xxi <- as.double(dtf$WV)
> yyi <- as.double(dtf$DC)
> # 式(2.28)
> r1 <- lm(yyi ~ 1 + xxi)
> print(summary(r1))
Call:
lm(formula = yyi ~ 1 + xxi)
Residuals:
Min 1Q Median 3Q Max
-0.59869 -0.14099 0.06059 0.17262 0.32184
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.13088 0.12599 1.039 0.31
xxi 0.24115 0.01905 12.659 7.55e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2361 on 23 degrees of freedom
Multiple R-squared: 0.8745, Adjusted R-squared: 0.869
F-statistic: 160.3 on 1 and 23 DF, p-value: 7.546e-12
> # 式(2.29)
> r2 <- lm(yyi ~ 1 + log(xxi))
> print(summary(r2))
Call:
lm(formula = yyi ~ 1 + log(xxi))
Residuals:
Min 1Q Median 3Q Max
-0.31619 -0.07685 0.02395 0.11139 0.23029
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.83036 0.11083 -7.493 1.3e-07 ***
log(xxi) 1.41677 0.06234 22.728 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1376 on 23 degrees of freedom
Multiple R-squared: 0.9574, Adjusted R-squared: 0.9555
F-statistic: 516.6 on 1 and 23 DF, p-value: < 2.2e-16
> # 式(2.30)
> r3 <- lm(yyi ~ 1 + I(1 / xxi))
> print(summary(r3))
Call:
lm(formula = yyi ~ 1 + I(1/xxi))
Residuals:
Min 1Q Median 3Q Max
-0.20547 -0.04940 0.01100 0.08352 0.12204
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.9789 0.0449 66.34 <2e-16 ***
I(1/xxi) -6.9345 0.2064 -33.59 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.09417 on 23 degrees of freedom
Multiple R-squared: 0.98, Adjusted R-squared: 0.9792
F-statistic: 1128 on 1 and 23 DF, p-value: < 2.2e-16
表3.7の値をCSV形式で入力したtable3_7.csvを、カレントディレクトリに置いておく。
> dtf <- read.csv("table3_7.csv", header = TRUE)
> xx1i <- 1 / as.double(dtf$RU)
> xx2i <- as.double(dtf$CPIDOT)
> yyi <- as.double(dtf$WDOT)
> n <- length(xx1i)
> k <- 3
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xx1i, xx2i), ncol = 3)
> mxxxtxxi <- solve(t(mxxx) %*% mxxx)
> mxbh <- mxxxtxxi %*% t(mxxx) %*% mxy
> mxyh <- mxxx %*% mxbh
> mxe <- mxy - mxyh
> sse2 <- t(mxe) %*% mxe
> s2 <- as.vector(sse2 / degf)
> s <- sqrt(s2)
> mxvv <- s2 * mxxxtxxi
> mxs <- matrix(sqrt(diag(mxvv)), ncol = 1)
> tj <- as.vector(mxbh / mxs)
> pj <- 2 * pt(abs(tj), degf, lower.tail = FALSE)
> bh <- as.vector(mxbh)
> yyhi <- as.vector(mxyh)
> ei <- yyi - yyhi
> yym <- mean(yyi)
> rr2 <- sum((yyhi - yym) ^ 2) / sum((yyi - yym) ^ 2)
> rrb2 <- 1 - (n - 1) / (n - k) * (1 - rr2)
> print(bh) # 推定した回帰モデルのパラメーター
[1] -4.6884718 16.2320984 0.8915398
> print(tj) # t値
[1] -8.938476 12.501712 12.623158
> print(pj) # p値
[1] 6.096213e-12 5.248937e-17 3.628214e-17
> print(rr2) # 決定係数
[1] 0.9509842
> print(rrb2) # 自由度修正済み決定係数
[1] 0.9490236
> print(s) # 誤差項の分散の不偏推定量の標準誤差
[1] 1.580779
