R(数値計算)

2023年5月31日 (水)

[R]回帰モデルの信頼区間と予測区間

lm関数により求めた回帰モデルの信頼区間と予測区間はpredict関数で求めることができる。以下の計算により、赤実線が求めた回帰直線。桃破線がその95%信頼区間、橙点線がその95%予測区間。

> # 説明変数xと目的変数y
> x <- c(1, 2, 3, 5, 6, 7, 8, 10, 11, 12)
> y <- c(0, 1, 1, 2, 4, 4, 6, 5, 8, 8)
> # 回帰直線を求める
> r <- lm(y ~ x)
> print(r)
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
-0.8567 0.7318
> # yの推定値
> yest <- fitted(r)
> # 信頼区間と予測区間の計算(後続の手順のためdata.frame化)
> dtf <- data.frame(x)
> rcon <- predict(r, dtf, interval = "confidence", level = 0.95)
> rpre <- predict(r, dtf, interval = "prediction", level = 0.95)
> rcon <- data.frame(rcon)
> rpre <- data.frame(rpre)
> # 図の作成(黒丸:観測値,赤実線:回帰直線,桃破線:信頼区間,橙点線:予測区間)
> plot(x, y, type = "n", asp = 1.0)
> lines(x, rpre$lwr, lty = "dotted", col = "orange", lwd = 2.0)
> lines(x, rpre$upr, lty = "dotted", col = "orange", lwd = 2.0)
> lines(x, rcon$lwr, lty = "dashed", col = "pink", lwd = 2.0)
> lines(x, rcon$upr, lty = "dashed", col = "pink", lwd = 2.0)
> lines(x, yest, col = "red", lwd = 2.0)
> points(x, y, pch = 20, col = "black")

R_conpre

> # 以下はpredict関数と手計算の計算結果の比較
> ah <- 0.95
> n <- length(x)
> k <- 2
> mxxx <- matrix(c(rep(1, n), x), ncol = 2)
> mxy <- matrix(y, ncol = 1)
> # 最小二乗推定量
> mxb <- solve(t(mxxx) %*% mxxx) %*% (t(mxxx) %*% mxy)
> print(mxb)
[,1]
[1,] -0.8567050
[2,] 0.7318008
> # yの推定値
> yest <- mxxx %*% mxb
> # 残差分散
> s <- sqrt(sum((y - yest) ^ 2) / (n - k))
> #
> s0 <- sh0 <- double(n)
> for (i in 1:n) {
+ mxx0 <- matrix(c(1, x[i]), ncol = 1)
+ s0[i] <- s * sqrt(t(mxx0) %*% solve(t(mxxx) %*% mxxx) %*% mxx0)
+ sh0[i] <- s * sqrt(1 + t(mxx0) %*% solve(t(mxxx) %*% mxxx) %*% mxx0)
+ }
> # 95%信頼区間
> ycon1 <- yest - s0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> ycon2 <- yest + s0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> # 95%予測区間
> ypre1 <- yest - sh0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> ypre2 <- yest + sh0 * qt((1. - ah) / 2, n - k, lower.tail = FALSE)
> # predict関数による信頼区間と予測区間
> print(data.frame(rcon$lwr, rcon$upr, rpre$lwr, rpre$upr))
rcon.lwr rcon.upr rpre.lwr rpre.upr
1 -1.1763909 0.9265825 -2.2315171 1.981709
2 -0.3152118 1.5290049 -1.4382140 2.652007
3 0.5349489 2.1424457 -0.6558465 3.333241
4 2.1772621 3.4273356 0.8728263 4.731771
5 2.9513451 4.1168542 1.6179064 5.450293
6 3.6831458 4.8486549 2.3497072 6.182094
7 4.3726644 5.6227379 3.0682286 6.927174
8 5.6575543 7.2650511 4.4667589 8.455846
9 6.2709951 8.1152118 5.1479929 9.238214
10 6.8734175 8.9763909 5.8182913 10.031517
> # 手計算による信頼区間と予測区間
> print(data.frame(ycon1, ycon2, ypre1, ypre2))
ycon1 ycon2 ypre1 ypre2
1 -1.1763909 0.9265825 -2.2315171 1.981709
2 -0.3152118 1.5290049 -1.4382140 2.652007
3 0.5349489 2.1424457 -0.6558465 3.333241
4 2.1772621 3.4273356 0.8728263 4.731771
5 2.9513451 4.1168542 1.6179064 5.450293
6 3.6831458 4.8486549 2.3497072 6.182094
7 4.3726644 5.6227379 3.0682286 6.927174
8 5.6575543 7.2650511 4.4667589 8.455846
9 6.2709951 8.1152118 5.1479929 9.238214
10 6.8734175 8.9763909 5.8182913 10.031517

2023年1月 3日 (火)

[R]回帰分析におけるAICを簡単に求める

回帰分析をlm関数で行い、その戻り値を用いてAIC関数を使う。以下は「情報量規準」(朝倉書店)に掲載の計算例(pp.59-61)を再現した結果。

観測値はpp.59-60に掲載されており、以下のとおり。これをtable_p059.csvと保存する。

 i,    x,     y
1, 0.00, 0.854
2, 0.05, 0.786
3, 0.10, 0.706
4, 0.15, 0.763
5, 0.20, 0.772
6, 0.25, 0.693
7, 0.30, 0.805
8, 0.35, 0.739
9, 0.40, 0.760
10, 0.45, 0.764
11, 0.50, 0.810
12, 0.55, 0.791
13, 0.60, 0.798
14, 0.65, 0.841
15, 0.70, 0.882
16, 0.75, 0.879
17, 0.80, 0.863
18, 0.85, 0.934
19, 0.90, 0.971
20, 0.95, 0.985

以下、計算結果。当該書籍では、上記の観測値を用いた多項式回帰モデルの残差分散とAICの一覧を表にまとめて掲載している(p.61)。

> dtf <- read.csv("table_p059.csv", header = TRUE)
> n <- nrow(dtf)
> for (k in 1:9) {
+ r <- lm(y ~ poly(x, k), data = dtf)
+ rss <- sum(r$re ^ 2)
+ cat(sprintf("次数:%d, 残差分散:%8.6f, AIC:%5.2f\n", k, rss / n, AIC(r)))
+ }
次数:1, 残差分散:0.002587, AIC:-56.38
次数:2, 残差分散:0.000922, AIC:-75.03
次数:3, 残差分散:0.000833, AIC:-75.04
次数:4, 残差分散:0.000737, AIC:-75.50
次数:5, 残差分散:0.000688, AIC:-74.89
次数:6, 残差分散:0.000650, AIC:-74.00
次数:7, 残差分散:0.000622, AIC:-72.89
次数:8, 残差分散:0.000607, AIC:-71.38
次数:9, 残差分散:0.000599, AIC:-69.66

2023年1月 1日 (日)

[R]回帰分析における対数尤度を簡単に求める

回帰分析をlm関数で行い、その戻り値を用いてlogLik関数を使う。以下は「情報量規準」(朝倉書店)に掲載の計算例(pp.59-61)を再現した結果。

観測値はpp.59-60に掲載されており、以下のとおり。これをtable_p059.csvと保存する。

 i,    x,     y
1, 0.00, 0.854
2, 0.05, 0.786
3, 0.10, 0.706
4, 0.15, 0.763
5, 0.20, 0.772
6, 0.25, 0.693
7, 0.30, 0.805
8, 0.35, 0.739
9, 0.40, 0.760
10, 0.45, 0.764
11, 0.50, 0.810
12, 0.55, 0.791
13, 0.60, 0.798
14, 0.65, 0.841
15, 0.70, 0.882
16, 0.75, 0.879
17, 0.80, 0.863
18, 0.85, 0.934
19, 0.90, 0.971
20, 0.95, 0.985

以下、計算結果。当該書籍では、上記の観測値を用いた多項式回帰モデルの残差分散と対数尤度の一覧を表にまとめて掲載している(p.61)。

> dtf <- read.csv("table_p059.csv", header = TRUE)
> n <- nrow(dtf)
> for (k in 1:9) {
+ r <- lm(y ~ poly(x, k), data = dtf)
+ rss <- sum(r$re ^ 2)
+ cat(sprintf("次数:%d, 残差分散:%8.6f, 対数尤度:%5.2f\n", k, rss / n, logLik(r)))
+ }
次数:1, 残差分散:0.002587, 対数尤度:31.19
次数:2, 残差分散:0.000922, 対数尤度:41.51
次数:3, 残差分散:0.000833, 対数尤度:42.52
次数:4, 残差分散:0.000737, 対数尤度:43.75
次数:5, 残差分散:0.000688, 対数尤度:44.44
次数:6, 残差分散:0.000650, 対数尤度:45.00
次数:7, 残差分散:0.000622, 対数尤度:45.44
次数:8, 残差分散:0.000607, 対数尤度:45.69
次数:9, 残差分散:0.000599, 対数尤度:45.83

2022年12月19日 (月)

[R]パスカルの三角形を求める

以下のスクリプトを実行する。1~12のパスカルの三角形を求めている。

for (i in 1:12) {
cat(sprintf("%2d: ", i))
for (j in 0:i) {
cat(sprintf("%3d ", choose(i, j)))
}
cat("\n")
}

出力

 1:   1   1 
2: 1 2 1
3: 1 3 3 1
4: 1 4 6 4 1
5: 1 5 10 10 5 1
6: 1 6 15 20 15 6 1
7: 1 7 21 35 35 21 7 1
8: 1 8 28 56 70 56 28 8 1
9: 1 9 36 84 126 126 84 36 9 1
10: 1 10 45 120 210 252 210 120 45 10 1
11: 1 11 55 165 330 462 462 330 165 55 11 1
12: 1 12 66 220 495 792 924 792 495 220 66 12 1

2022年9月 8日 (木)

[R]丸め誤差を考慮して数値の比較を行う

コンピューターは内部では浮動小数点演算を行っている都合上、丸め誤差は避けられない。それは実数同士の演算でよく見られる。

> (1 + 2) == 3
[1] TRUE
> (0.1 + 0.2) == 0.3
[1] FALSE
> print(sprintf("%.20f", 0.1 + 0.2))
[1] "0.30000000000000004441"
> print(sprintf("%.20f", 0.3))
[1] "0.29999999999999998890"

Rには標準でVisual BasicのDecimal型に相当するベクトルの型は存在せず、計算時に工夫する方法がある。

簡単に回避する方法としてsignif関数を使う方法がある。この関数は指定した有効数字を指定した桁まで丸めることができるため(デフォルトは6桁)、比較時にこの関数を使用すればよい。

> signif(0.1 + 0.2) == 0.3
[1] TRUE
> signif(0.1 + 0.2, digits = 16) == 0.3
[1] TRUE
> signif(0.1 + 0.2, digits = 17) == 0.3
[1] FALSE
> print(sprintf("%.20f", signif(0.1 + 0.2, digits = c(15, 16, 17, 18))))
[1] "0.29999999999999998890" "0.29999999999999998890" "0.30000000000000004441" "0.30000000000000004441"
> print(sprintf("%.20f", signif(0.3, digits = c(15, 16, 17, 18))))
[1] "0.29999999999999998890" "0.29999999999999998890" "0.29999999999999998890" "0.29999999999999998890"

ちなみに、round関数でも同様のことができるが、round関数とsignif関数ではdigitsオプションの意味が違うため注意。この手の処理でround関数を使うことは推奨しない。

> round(0.1 + 0.2, digits = 15) == 0.3
[1] TRUE
> round(0.1 + 0.2, digits = 16) == 0.3
[1] FALSE
> print(sprintf("%.20f", round(0.1 + 0.2, digits = 14:17)))
[1] "0.29999999999999998890" "0.29999999999999998890" "0.30000000000000004441" "0.30000000000000004441"
> print(sprintf("%.20f", round(0.3, digits = 14:17)))
[1] "0.29999999999999998890" "0.29999999999999998890" "0.29999999999999998890" "0.29999999999999998890"

丸め誤差を完全に自動でうまく処理できる方法はなく、その計算で求められる精度を基に、その都度工夫をする必要がある。

2022年5月15日 (日)

[R]単回帰分析(入門 統計解析 医学・自然科学編、東京書籍、pp.236-237)

本計算に使用するデータは、以下のとおり(p.236)。

 i,  x,  y
1, 7, 6
2, 7, 9
3, 12, 10
4, 11, 13
5, 10, 13
6, 5, 7
7, 9, 11
8, 11, 14
9, 11, 15
10, 12, 7
11, 11, 13
12, 15, 14
13, 14, 10
14, 16, 16
15, 8, 8
16, 2, 8
17, 14, 8
18, 8, 12
19, 12, 16
20, 16, 12

これをメモ帳に貼り付けてファイル「table10-1.csv」で保存し、カレントディレクトリに置いておく。

> dtf <- read.csv("table10-1.csv", header = TRUE)
> xi <- dtf$x
> yi <- dtf$y
> n <- nrow(dtf)
> p <- 2
> xm <- mean(xi)
> ym <- mean(yi)
> ssx2 <- sum((xi - xm) ^ 2)
> ssy2 <- sum((yi - ym) ^ 2)
> ssxy <- sum((xi - xm) * (yi - ym))
> b1 <- ssxy / ssx2
> b0 <- ym - b1 * xm
> yt <- b0 + b1 * xi
> sse <- sum((yi - yt) ^ 2)
> sh2 <- sse / (n - p)
> se <- sqrt(sh2)
> seb1 <- sqrt(sh2 / sum((xi - xm) ^ 2))
> tb1 <- b1 / seb1
> cat(sprintf("回帰式 y~ = %.4fx + %.3f\n", b1, b0))
回帰式 y~ = 0.4468x + 6.387
> cat(sprintf("誤差分散σ2の不偏推定量 σ2^ = %.2f\n", sh2))
誤差分散σ2の不偏推定量 σ2^ = 7.61
> cat(sprintf("回帰値の標準誤差 s.e. = %.2f\n", se))
回帰値の標準誤差 s.e. = 2.76
> cat(sprintf("β1の推定値の標準誤差 s.e.(b1) = %.3f\n", seb1))
β1の推定値の標準誤差 s.e.(b1) = 0.173
> cat(sprintf("β1の推定値のt値 t = %.2f\n", tb1))
β1の推定値のt値 t = 2.59
> cat(sprintf("t0.025(18) = %.2f\n", qt(1 - 0.025, n - p)))
t0.025(18) = 2.10

2022年2月 3日 (木)

[R]毎年の太陽黒点数(Wolfer sunspot number)のデータ

パッケージTSSSに格納されている。データ名はSunspot。1749年から1979年までの毎年の値で、データ数は231個。

> library(TSSS)
> Sunspot
エラー: オブジェクト 'Sunspot' がありません
> data(Sunspot)
> Sunspot
Time Series:
Start = 1749
End = 1979
Frequency = 1
[1] 80.9 83.4 47.7 47.8 30.7 12.2 9.6 10.2 32.4 47.6 54.0 62.9 85.9 61.2 45.1
[16] 36.4 20.9 11.4 37.8 69.8 106.1 100.8 81.6 66.5 34.8 30.6 7.0 19.8 92.5 154.4
[31] 125.9 84.8 68.1 38.5 22.8 10.2 24.1 82.9 132.0 130.9 118.1 89.9 66.6 60.0 46.9
[46] 41.0 21.3 16.0 6.4 4.1 6.8 14.5 34.0 45.0 43.1 47.5 42.2 28.1 10.1 8.1
[61] 2.5 0.1 1.4 5.0 12.2 13.9 35.4 45.8 41.1 30.1 23.9 15.6 6.6 4.0 1.8
[76] 8.5 16.6 36.3 49.6 64.2 67.0 70.9 47.8 27.5 8.5 13.2 56.9 121.5 138.3 103.2
[91] 85.7 64.6 36.7 24.2 10.7 15.0 40.1 61.5 98.5 124.7 96.3 66.6 64.5 54.1 39.0
[106] 20.6 6.7 4.3 22.7 54.8 93.8 95.8 77.2 59.1 44.0 47.0 30.5 16.3 7.3 37.6
[121] 74.0 139.0 111.2 101.6 66.2 44.7 17.0 11.3 12.4 3.4 6.0 32.3 54.3 59.7 63.7
[136] 63.5 52.2 25.4 13.1 6.8 6.3 7.1 35.6 73.0 85.1 78.0 64.0 41.8 26.2 26.7
[151] 12.1 9.5 2.7 5.0 24.4 42.0 63.5 53.8 62.0 48.5 43.9 18.6 5.7 3.6 1.4
[166] 9.6 47.4 57.1 103.9 80.6 63.6 37.6 26.1 14.2 5.8 16.7 44.3 63.9 69.0 77.8
[181] 64.9 35.7 21.2 11.1 5.7 8.7 36.1 79.7 114.4 109.6 88.8 67.8 47.5 30.6 16.3
[196] 9.6 33.2 92.6 151.6 136.3 134.7 83.9 69.4 31.5 13.9 4.4 38.0 141.7 190.2 184.8
[211] 159.0 112.3 53.9 37.5 27.9 10.2 15.1 47.0 93.8 105.9 105.5 104.5 66.6 68.9 38.0
[226] 34.5 15.5 12.6 27.5 92.5 155.4
> plot(Sunspot)


Sunspot

2021年9月11日 (土)

[R]「Rによる時系列モデリング入門」に掲載の時系列解析を行う

パッケージTSSSを使用する。

> install.packages("TSSS")
パッケージを ‘C:/Users/○○/Documents/R/win-library/○○’ 中にインストールします
(‘lib’ が指定されていないため)
URL 'https://○○/TSSS_1.3.1.zip' を試しています
Content type 'application/zip' length 852282 bytes (832 KB)
downloaded 832 KB
package ‘TSSS’ successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\○○\AppData\Local\Temp\○○\downloaded_packages
> library(TSSS)

最新情報はhttps://jasp.ism.ac.jp/ism/TSSS/を参照のこと。

2021年8月 6日 (金)

[R]あるハードウェアの毎月の卸売高(WHARD(Wholesale hardware)データ)

パッケージTSSSに格納されている。パッケージをインストールして読み出す。

> install.packages("TSSS")
(表示省略)
> library(TSSS)
> data(WHARD)
> WHARD
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1967 626 614 689 686 723 778 711 824 793 831 775 689
1968 692 718 757 769 791 809 836 878 856 935 850 763
1969 761 796 830 902 910 932 931 908 934 995 865 822
1970 763 778 841 845 863 952 909 899 952 963 893 831
1971 773 803 918 967 963 1065 1014 1051 1054 1051 1039 960
1972 930 956 1072 1023 1136 1181 1088 1247 1164 1251 1218 1062
1973 1114 1088 1253 1254 1354 1349 1305 1420 1313 1481 1387 1284
1974 1310 1262 1446 1573 1634 1612 1591 1640 1590 1696 1456 1296
1975 1311 1232 1274 1388 1374 1443 1466 1454 1538 1587 1406 1341
1976 1351 1367 1553 1588 1591 1703 1643 1711 1731 1678 1678 1580
1977 1515 1544 1817 1838 1925 2017 1898 2068 1961 2027 1974 1820
1978 1796 1708 2021 2102 2306 2360 2247 2412 2159 2455 2250 2057
1979 2142 1984 2319 2374 2592 2461 2524 2678 2399 2794 2415
> plot(WHARD)

Whard

2021年6月17日 (木)

[R]関数の最大値・最小値を得る

optim関数を使う。最小値を得るための実行例は以下のとおり。関数 f(x) = x ^ 2 - 5 の最小値を求めている。図からf(x)が最小となるのはx = 0でその時のf(x)は5である。簡単な計算であれば、methodオプションはBFGSを指定する。計算には適当な初期値(この例では、f(x)が最小となるxの探索のための最初の値)を与える必要があり、以下の例では10としている。

Optimmin
> f <- function(xf) xf ^ 2 - 5
> x <- seq(-5, 5, by = 0.2)
> y <- f(x)
> optim(10, f, method = "BFGS")
$par
[1] 2.441602e-11
$value
[1] -5
$counts
function gradient
6 3
$convergence
[1] 0
$message
NULL

戻り値のparに最小となる関数の引数が、valueにはその最小となるf(x)値が格納されている。

逆に、関数の最大値を得るためには、controlオプションにリストとしてfnscaleに負の値を与える。以下の例では、関数f(x) = -(x ^ 2) - 5の最大値を求めており、図より、x=0の時に最大値f(x)=-5を得ることが明らかである。

Optimmax
> f <- function(xf) -(xf ^ 2) - 5
> x <- seq(-5, 5, by = 0.2)
> y <- f(x)
> optim(10, f, control = list(fnscale = -1), method = "BFGS")
$par
[1] 2.441602e-11
$value
[1] -5
$counts
function gradient
6 3
$convergence
[1] 0
$message
NULL

戻り値parに関数が最大となる場合の引数の値(=0)、その最大となった関数の値(=-5)が格納されている。

より以前の記事一覧