［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

［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

［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 

［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"

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

［R］単回帰モデルにおける回帰係数のp値を求める（44の例題で学ぶ計量経済学、オーム社、p.181）

データは以下のとおり（p.97の表3.2）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- nrow(dtf)
> k <- 2
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi), ncol = 2)
> mxb <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> b <- as.vector(mxb)
> yyhi <- mxxx %*% mxb
> sssse <- sum((yyi - yyhi) ^ 2)
> sh2 <- sssse / degf
> mxcc <- sh2 * solve(t(mxxx) %*% mxxx)
> seb <- sqrt(diag(mxcc))
> tb <- b / seb
> pb <- 2 * pt(abs(tb), degf, lower.tail = FALSE)
> rr2 <- sum((yyhi - mean(yyi)) ^ 2) / sum((yyi - mean(yyi)) ^ 2)
> arr2 <- 1 - (sssse / (n - k)) / (sum((yyi - mean(yyi)) ^ 2) / (n - 1))
> cat(sprintf("回帰式 Y^i =  %.1f  +  %.1f Xi\n", b[1], b[2]))
回帰式 Y^i =  0.5  +  1.1 Xi
> cat(sprintf("            (%.3f) (%.3f)\n", tb[1], tb[2]))
            (0.430) (5.185)
> cat(sprintf("αの推定値の t値 tα = %.3f\n", tb[1]))
αの推定値の t値 tα = 0.430
> cat(sprintf("αの推定値の p値 pα = %.3f\n", pb[1]))
αの推定値の p値 pα = 0.709
> cat(sprintf("βの推定値の t値 tβ = %.3f\n", tb[2]))
βの推定値の t値 tβ = 5.185
> cat(sprintf("βの推定値の p値 pβ = %.3f\n", pb[2]))
βの推定値の p値 pβ = 0.035

［R］単回帰モデルにおけるt分布による切片の仮説検定（両側検定）を行う（44の例題で学ぶ計量経済学、オーム社、p.181）

データは以下のとおり（p.97の表3.2）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- nrow(dtf)
> k <- 2
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi), ncol = 2)
> mxb <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> b <- as.vector(mxb)
> yyhi <- mxxx %*% mxb
> ssu2 <- sum((yyi - yyhi) ^ 2)
> sh2 <- ssu2 / (n - k)
> ssxx <- sum((xxi - mean(xxi)) ^ 2)
> sa2 <- sh2 * sum(xxi ^ 2) / (n * ssxx)
> sa <- sqrt(sa2)
> ta <- b[1] / sa
> cat(sprintf("残差分散 σ^2 = %.1f\n", sh2))
残差分散 σ^2 = 0.9
> cat(sprintf("αの推定値の分散 sa^2 = %.2f\n", sa2))
αの推定値の分散 sa^2 = 1.35
> cat(sprintf("αの推定値の標準誤差 sα = %.5f\n", sa))
αの推定値の標準誤差 sα = 1.16190
> cat(sprintf("αの推定値のt値 tα = %.3f\n", ta))
αの推定値のt値 tα = 0.430
> cat(sprintf("両側検定 t0.05(2) = %.3f\n", abs(qt(0.05 / 2, n - k))))
両側検定 t0.05(2) = 4.303

［R］単回帰モデルにおけるt分布による傾きの仮説検定（両側検定、片側検定）を行う（44の例題で学ぶ計量経済学、オーム社、pp.178-180）

データは以下のとおり（p.97）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- nrow(dtf)
> k <- 2
> degf <- n - k
> mxy <- matrix(yyi, ncol = 1)
> mxxx <- matrix(c(rep(1.0, n), xxi), ncol = 2)
> mxb <- solve(t(mxxx) %*% mxxx) %*% t(mxxx) %*% mxy
> b <- as.vector(mxb)
> yyhi <- mxxx %*% mxb
> ssu2 <- sum((yyi - yyhi) ^ 2)
> sh2 <- ssu2 / (n - k)
> ssxx <- sum((xxi - mean(xxi)) ^ 2)
> sb2 <- sh2 / ssxx
> sb <- sqrt(sb2)
> tb <- b[2] / sb
> cat(sprintf("残差分散 σ^2 = %.1f\n", sh2))
残差分散 σ^2 = 0.9
> cat(sprintf("βの推定値の分散 sb^2 = %.3f\n", sb2))
βの推定値の分散 sb^2 = 0.045
> cat(sprintf("βの推定値の標準誤差 sβ = %.5f\n", sb))
βの推定値の標準誤差 sβ = 0.21213
> cat(sprintf("βの推定値のt値 tβ = %.3f\n", tb))
βの推定値のt値 tβ = 5.185
> cat(sprintf("両側検定 t0.05(2) = %.3f\n", abs(qt(0.05 / 2, n - k))))
両側検定 t0.05(2) = 4.303
> cat(sprintf("片側検定 t0.05(2) = %.3f\n", qt(0.05, n - k, lower.tail = FALSE)))
片側検定 t0.05(2) = 2.920

［R］最小二乗法による回帰直線の傾き、切片、決定係数、自由度調整済み決定係数を求める（44の例題で学ぶ計量経済学、オーム社、pp.96-108）

データは以下のとおり（p.97）。xxiは説明変数、yyiは目的変数（被説明変数）。これをメモ帳に貼り付けて「table3_2.csv」というテキストファイルで保存をして、カレントディレクトリに置いておく。

i, xxi, yyi
1,   2,   3
2,   4,   5
3,   6,   6
4,   8,  10

以下、計算。

> dtf <- read.csv("table3_2.csv", header = TRUE)
> xxi <- dtf$xxi
> yyi <- dtf$yyi
> n <- length(xxi)
> ssxy <- sum((xxi - mean(xxi)) * (yyi - mean(yyi)))
> ssxx <- sum((xxi - mean(xxi)) ^ 2)
> b <- ssxy / ssxx
> xxm <- mean(xxi)
> yym <- mean(yyi)
> a <- yym - b * xxm
> yyei <- a + b * xxi
> ui <- yyi - yyei
> ssu2 <- sum(ui ^ 2)
> ssyy <- sum((yyi - yym) ^ 2)
> ssyeye <- sum((yyei - yym) ^ 2)
> rr2 <- ssyeye / ssyy
> sig2 <- ssu2 / (n - 2)
> sy2 <- ssyy / (n - 1)
> arr2 <- 1 - sig2 / sy2
> cat(sprintf("傾き β=%.2f\n", b))
傾き β=1.10
> cat(sprintf("切片 α=%.2f\n", a))
切片 α=0.50
> cat(sprintf("決定係数 R^2=%.3f\n", rr2))
決定係数 R^2=0.931
> cat(sprintf("自由度調整済み決定係数 adj.R^2=%.3f\n", arr2))
自由度調整済み決定係数 adj.R^2=0.896

［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

［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)