［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］指定された数の空白（0x20）からなる文字列を得る

C#でいうところのStrings.Spaceメソッド。そのような機能を持つ関数はRには標準で搭載されていないためその機能を持つ関数を自作する。以下の3行がその自作関数。

space <- function(n) {
  return(paste(rep(" ", n), collapse = ""))
}

試しに使ってみる

> space(5)
[1] "     "
> paste("A", space(1), "B", space(2), "C", sep = "")
[1] "A B  C"

［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］数値を文字列に変換する

簡単に変換するならas.character関数、変換時に書式を指定したいのであればsprintf関数を使う。

> n <- c(1, 2, -4)
> d <- c(1.1, 0.1, -2.2)
> as.character(n)
[1] "1"  "2"  "-4"
> as.character(d)
[1] "1.1"  "0.1"  "-2.2"
> sprintf("%03d", n)
[1] "001" "002" "-04"
> sprintf("%+04d", n)
[1] "+001" "+002" "-004"
> sprintf("%8.4f", d)
[1] "  1.1000" "  0.1000" " -2.2000"

［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

[Octave］t分布におけるパーセント点を求める

tinv関数を使う。以下は自由度5のt分布において、上側5%点（=2.015048…）、両側5%点（=2.570582…）を求めた例。

>> pkg load statistics
>> tinv(1 - 0.05, 5)
ans = 2.0150
>> tinv(1 - 0.05 / 2, 5)
ans = 2.5706

tinv関数は下側パーセント点の値しか得られないため、上側パーセント点の値を得るには1からその確率の値を引くことで得られる。両側パーセント点は下側パーセント点の2倍のため、それを2倍している。

［Octave］t分布における分布関数の値を求める

tcdf関数を使う。以下は自由度5のt分布において、確率変数Xが2における上側（=0.05096974…）、両側（=0.1019395…）の分布関数の値を求めた例。

>> pkg load statistics
>> 1 - tcdf(2, 5)
ans = 0.050970
>> 2 * (1 - tcdf(2, 5))
ans = 0.1019

tcdf関数は下側の値しか得られないため、上側の値を得るには1からその値を引くことで得られる。両側は下側の2倍のため、それを2倍している。

［Octave］代数方程式の根を求める

roots関数を使う。以下は2x3+3x2+8x-5=0の3つの根（0.5,-1+2i,-1-2i）を求めた例。虚部も求めることができる。引数には、次数の高いほうから係数を行列でまとめて与える（この場合は[2 3 8 -5]とする）。

>> roots([2 3 8 -5])
ans =
  -1.0000 + 2.0000i
  -1.0000 - 2.0000i
   0.5000 +      0i

戻り値は複素数型。複素数型から実部だけを取り出すにはreal関数を、虚部だけを取り出すにはimag関数を使う。

>> ri = roots([2 3 8 -5]);
>> iscomplex(ri)
ans = 1
>> real(ri)
ans =
  -1.0000
  -1.0000
   0.5000
>> imag(ri)
ans =
   2
  -2
   0

［Octave］固有値と固有ベクトルを求める

eig関数を使う。固有値と固有ベクトルが行列でそれぞれ戻り値になるため、戻り値は2つ指定する。

>> mxaa = [3, 1, 1; 1, 2, 0; 1, 0, 2]
mxaa =
   3   1   1
   1   2   0
   1   0   2
>> [mxv, mxl] = eig(mxaa)
mxv =
  -5.7735e-01  -1.7174e-16   8.1650e-01
   5.7735e-01  -7.0711e-01   4.0825e-01
   5.7735e-01   7.0711e-01   4.0825e-01
mxl =
Diagonal Matrix
   1   0   0
   0   2   0
   0   0   4

この例では固有値と固有ベクトルは3つあり、例えば3番目の固有値と固有ベクトルを取り出すには、以下のようにする。

>> mxl(3, 3)
ans = 4.0000
>> mxv(:, 3)
ans =
   0.8165
   0.4082
   0.4082

［R］固有値と固有ベクトルを求める

eigen関数を使う。戻り値はリストで、valuesに固有値が、vectorsにその固有値に対応する固有ベクトルが含まれる。

> mxaa <- matrix(c(3, 1, 1, 1, 2, 0, 1, 0, 2), nrow = 3)
> print(mxaa)
     [,1] [,2] [,3]
[1,]    3    1    1
[2,]    1    2    0
[3,]    1    0    2
> lam <- eigen(mxaa)
> print(lam)
eigen() decomposition
$values
[1] 4 2 1
$vectors
          [,1]       [,2]       [,3]
[1,] 0.8164966  0.0000000  0.5773503
[2,] 0.4082483 -0.7071068 -0.5773503
[3,] 0.4082483  0.7071068 -0.5773503

この例では固有値と固有ベクトルは3つあり、例えば3番目の固有値と固有ベクトルを取り出すには、以下のようにする。

> lam$values[3]
[1] 1
> lam$vectors[, 3]
[1]  0.5773503 -0.5773503 -0.5773503

［R］単位行列を作る

diag関数を使う。引数を一つだけ指定した場合は、n次の単位行列を作成する。

> diag(3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
> diag(4)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

［R］コレスキー分解する

chol関数を使う。

> mx0 <- matrix(c(1, 0, 1, 0, 2, 0, 1, 0, 3), 3, byrow = TRUE)
> mx0
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    2    0
[3,]    1    0    3
> mx <- chol(mx0)
> mx
     [,1]     [,2]     [,3]
[1,]    1 0.000000 1.000000
[2,]    0 1.414214 0.000000
[3,]    0 0.000000 1.414214

得られた結果が正しいか、確認する。

> t(mx) %*% mx
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    2    0
[3,]    1    0    3

上記と同じ計算（A^T A）はcrossprod関数を使うと同じ結果を得ることができる。

> crossprod(mx)
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    2    0
[3,]    1    0    3

［R］複素数を使う

complex型ベクトルを使う。complex関数で簡単に作ることができる。第一引数には作成するベクトルの個数、第二引数には実部の値、第三引数には虚部の値を与える。

以下は、2+3iを作成した例。

> ri <- complex(1, 2, 3)
> print(ri)
[1] 2+3i

四則演算や累乗は、実数と同じ演算子を使用して簡単に計算することができる。

> ri1 <- complex(1, 3, -2)
> ri2 <- complex(1, 2, 5)
> print(ri1)
[1] 3-2i
> print(ri2)
[1] 2+5i
> print(ri1 + ri2)
[1] 5+3i
> print(ri1 - ri2)
[1] 1-7i
> print(ri1 * ri2)
[1] 16+11i
> ri <- complex(1, 1, 1)
> print(ri)
[1] 1+1i
> print(ri ^ 4)
[1] -4+0i

戻り値は複素数型。複素数型から実部だけを取り出すにはRe関数を、虚部だけを取り出すにはIm関数を使う。これら関数の戻り値はそれぞれ実数型。

> ri <- c(complex(1, 2, 3), complex(1, 4, 5))
> print(ri)
[1] 2+3i 4+5i
> class(ri)
[1] "complex"
> Re(ri)
[1] 2 4
> Im(ri)
[1] 3 5
> class(Re(ri))
[1] "numeric"
> class(Im(ri))
[1] "numeric"
> Re(ri) + Im(ri)
[1] 5 9

［R］代数方程式の根を求める

polyroot関数を使う。以下は2x³+3x²+8x-5=0の3つの根（0.5,-1+2i,-1-2i）を求めた例。虚部も求めることができる。引数には、次数の低いほうから係数をベクトルでまとめて与える（この場合はc(-5,8,3,2)とする）。

> polyroot(c(-5, 8, 3, 2))
[1]  0.5+0i -1.0+2i -1.0-2i

戻り値は複素数型。複素数型から実部だけを取り出すにはRe関数を、虚部だけを取り出すにはIm関数を使う。これら関数の戻り値はそれぞれ実数型。

> ri <- polyroot(c(-5, 8, 3, 2))
> class(ri)
[1] "complex"
> Re(ri)
[1]  0.5 -1.0 -1.0
> Im(ri)
[1]  3.372483e-16  2.000000e+00 -2.000000e+00
> class(Re(ri))
[1] "numeric"
> class(Im(ri))
[1] "numeric"

［Octave］t分布における確率密度関数の値を求める

tpdf関数を使う。以下は確率変数が-5～5の範囲において、自由度1と5のt分布における確率密度関数の値を求めてグラフにした例。

>> pkg load statistics
>> xx = -5:0.1:5;
>> hold off
>> plot(xx, tpdf(xx, 5), 'b')
>> hold on
>> plot(xx, tpdf(xx, 1), 'k')
>> text(1, 0.3, "n=5")
>> text(0.8, 0.1, "n=1")

［R］行列から非対角成分だけを取り出す

row関数とcol関数を使うと、ベクトル形式で簡単に取り出すことができる。

> mx <- matrix(1:16, 4, 4)
> mx
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> mx[row(mx) != col(mx)]
 [1]  2  3  4  5  7  8  9 10 12 13 14 15

応用として、上の例で2、7、12（対角成分のそれぞれ1行下の成分）を取り除くには、以下のようにする。

> mx[row(mx) != col(mx) + 1]
 [1]  1  3  4  5  6  8  9 10 11 13 14 15 16

［R］ファイル、フォルダー（ディレクトリ）の存在を確認する

ファイルの存在の有無を確認するにはfile.exists関数を、フォルダーの場合はdir.exists関数を使う。Rはフォルダー（ディレクトリ）の区切りを示す記号に「￥」（円マーク）と「／」（スラッシュ）の両方を使うことができる。

> file.exists("C:/Windows/win.ini")
[1] TRUE
> file.exists("C:/Windows/win.inii")
[1] FALSE
> dir.exists("C:/Windows")
[1] TRUE
> dir.exists("C:/Windowss")
[1] FALSE

file.exists関数は、フォルダーを指定した場合、そのフォルダー名の最後に区切り記号を付けないとTRUEを返すので注意。

> file.exists("C:/Windows")
[1] TRUE
> file.exists("C:/Windows/")
[1] FALSE

dir.exists関数はファイルであればFALSE、フォルダーであればTRUEを返す。

> dir.exists("C:/Windows")
[1] TRUE
> dir.exists("C:/Windows/")
[1] TRUE
> dir.exists("C:/Windows/win.ini")
[1] FALSE

［R］文字列の前後の空白を取り除く

trimws関数を使う。オプションを特に指定しなければ、文字列の前後の空白（とタブ（0x09）と復帰（CR,\r,0x0d）と改行（LF,\n,\x0a））を取り除く。第二引数に"left"と"right"をそれぞれ指定すると、前の空白だけ、後ろの空白だけを取り除く。デフォルトは両方（both）となっている。

> s <- "   Hello.  "
> trimws(s)
[1] "Hello."
> trimws(s, "left")
[1] "Hello.  "
> trimws(s, "right")
[1] "   Hello."

なお、いわゆる全角空白（0x8140）はそのままでは取り除くことができない。

> s <- "　あいうえお  "
> trimws(s)
[1] "　あいうえお"

ヘルプによると、whitespaceオプションに取り除く文字を正規表現で指定し、その指定に従って取り除いているとのこと。デフォルトは"[ \t\r\n]"。これに全角空白分を加えた文字列を指定すると、全角空白も同時に取り除くことができる。以下は、シフトJISの環境で行った例。

> trimws(s,  whitespace = "[ \t\r\n]")
[1] "　あいうえお"
> trimws(s,  whitespace = "[ \t\r\n　]")
[1] "あいうえお"
> trimws(s,  whitespace = "[ \t\r\n\x81\x40]")
[1] "あいうえお"