R入門

数値計算

最終更新：2024年12月07日 12:45

r-intro

- view

管理者のみ編集可

目次
丸め誤差を考慮して数値の比較を行う
関数の最大値・最小値を得る
関数の最小値（最大値）を検索する
複数の引数を持つ関数の値の最小値（最大値）を求める
サンプルデータを確認する
パスカルの三角形を求める
回帰分析における対数尤度を簡単に求める
回帰分析におけるAICを簡単に求める
回帰モデルの信頼区間と予測区間

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

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

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

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

関数の最大値・最小値を得る

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

> 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を得ることが明らかである。

> 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)が格納されている。

関数の最小値（最大値）を検索する

optimize関数を使う。以下は関数 f(x) = (x-2)^2 - 10 の最小値を検索した例。xの最小値は式と下図より10であり、その時の説明変数 x=2 である。戻り値はリストで、objectiveに最小値が、その時の説明変数の値がminimumに格納されている。

> f <- function(x) return((x - 2) ^ 2 + 10)
> plot(-5:5, f(-5:5), type = "o")
> optimize(f, c(-5, 5))
$minimum
[1] 2
$objective
[1] 10

以下は関数 f(x) = -(x-2)^2 - 10 の最大値を検索した例。xの最大値は式と下図より20であり、その時の説明変数 x=2 である。最大値を検索する場合はmaximumオプションをTRUEにする。

> f <- function(x) return(-(x - 2) ^ 2 + 20)
> plot(-5:5, f(-5:5), type = "o")
> optimize(f, c(-5, 5), maximum = TRUE)
$maximum
[1] 2
$objective
[1] 20

少しふざけて、上図の範囲（x > -5、x < 5）で最小値を検索してみる。2番目の引数に説明変数の検索範囲をベクトル（最小値、最大値）で与える。図のとおりに、説明変数がおおよそ-5のときにおおよそ-29であると求まる（グラフの左端のこと）。

> optimize(f, c(-5, 5))
$minimum
[1] -4.999944
$objective
[1] -28.99922

optimize関数の3番目以降の引数には、その検索を行う関数の引数を指定することができる。以下の例では検索する関数 f(x, a) であり、最小値検索のための説明変数はその関数の1番目にして、2番目以降は任意の値を設定することができる。

> f <- function(x, a) return((x - a) ^ 2 + 10)
> optimize(f, c(-5, 5), a = 2)
$minimum
[1] 2
$objective
[1] 10
> optimize(f, c(-5, 5), a = 3)
$minimum
[1] 3
$objective
[1] 10

複数の引数を持つ関数の値の最小値（最大値）を求める

optim関数を使う。使用時には適切な初期値を与える必要がある。数値的に求めるため、必ずしも期待通りの結果が得られるわけではないので注意。

以下はf1とf2というそれぞれ2つずつ引数を持つ関数の最小値と最大値を求めた例。f1の値が最小となるのはf1(5, 7)、f2の値が最大値となるのはf2(5, 7)であることは明らか。optim関数はデフォルトでは関数の値が最小値となる引数を求めるが、最大値となる引数を求めたい場合はcontrolオプションにlist(fnscale = -1)を指定すること。

> f1 <- function(x) {(x[1] - 5) ^ 2 + (x[2] - 7) ^ 2}
> f2 <- function(x) {-((x[1] - 5) ^ 2 + (x[2] - 7) ^ 2)}
> r <- optim(c(0, 0), f1, method = "BFGS")
> print(r)
$par
[1] 5 7
$value
[1] 1.31369e-22
$counts
function gradient 
      12        3 
$convergence
[1] 0
$message
NULL

f2関数の最小値は負の無限大であり、fnscaleに-1を指定しない状態では求めた引数（以下の$par）がまともに求まっていないことが明らか。fnscaleに-1を指定することで、期待通りの結果が得られる。

> r <- optim(c(0, 0), f2, method = "BFGS")
> print(r)
$par
[1] -4.758163e+13 -1.858368e+13
$value
[1] -2.609365e+27
$counts
function gradient 
      28       28 
$convergence
[1] 0
$message
NULL
> r <- optim(c(0, 0), f2, method = "BFGS", control = list(fnscale = -1))
> print(r)
$par
[1] 5 7
$value
[1] -1.31369e-22
$counts
function gradient 
      12        3 
$convergence
[1] 0
$message
NULL

サンプルデータを確認する

Rには、数値計算に使えるサンプルデータが標準で収録されている。標準で読み込まれるパッケージ（datasets）に多数含まれている。含まれているデータを確認するには、data関数を使う。

> data()

これらのデータはそのまますぐに使うことができる。

> AirPassengers
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
（以下、表示省略）

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

以下のスクリプトを実行する。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

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

回帰分析を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

回帰分析における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

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

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

> # 以下は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

「数値計算」をウィキ内検索

添付ファイル

人気記事ランキング

最近更新されたページ

新規Wikiランキング

最近作成されたWikiのアクセスランキングです。見るだけでなく加筆してみよう！

人気Wikiランキング

atwikiでよく見られているWikiのランキングです。新しい情報を発見してみよう！

全体ページランキング

最近アクセスの多かったページランキングです。話題のページを見に行こう！