情報科学と人工知能のノート

初等的な知識から最新論文の解説まで色々集めていきます．備忘録兼用．

R 備忘録 2013/05/10-3

統計学 R

中心極限定理をシミュレーション

平均 \mu、分散 \sigma の分布に従う独立した n 個の確率変数 X_1, \cdots, X_n について \bar{X} = \frac{1}{n} \sum_{i=1}^n X_i と定義する。
n が大きいとき、 \bar{X} は平均 \mu、分散 \frac{\sigma^2}{n} の正規分布に従う。

とのことです。つまり \bar{X} をたくさん作ってヒストグラムを作ったら正規分布の密度関数の曲線が現れるわけです。
0-1 一様分布で n=10 でシミュレーションしてみます。

まず、一様乱数の平均値を 10000 個ほど生成します。

> vec <- mapply(function (n) return (mean(runif(10))), 1:10000)
> summary(vec)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.2085 0.4391 0.5017 0.5010 0.5625 0.8141

最小値と最大値を考慮して、わかり易さ・やり易さのために [0:1] 区間を使うことにして、ヒストグラムを作ります。

> hist(vec, breaks=seq(0, 1, 0.01), probability=T, ylim=c(0, 5))

probability=TRUE を指定すると密度関数のように（おそらく、面積が 1 になるように）高さの単位を調節してくれます。
ちなみに ylim=c(0, 5) の 5 はヒストグラムを見て後から自分で調節しました。

0-1 一様分布の平均は 0.5、分散は 1/12 なので、n=10 のとき、\bar{X} は平均 0.5、分散 1/120 の正規分布に従います。

> par(new=T)
> plot(function(x) return (dnorm(x, mean=0.5, sd=sqrt(1/120))), xlim=c(0, 1), ylim=c(0, 5), xlab="", ylab="")

f:id:ti2236:20130510171903p:plain
それらしいグラフが表示されました。
ちなみに par はグラフィックスパラメータを取得／操作するための関数で、属性 new が TRUE のとき、新しい図が前の図に上書きされます。
new は図を描くたびに FALSE に戻るため、上書きするたびに設定してあげる必要があります。

分布の 100p パーセント点を求める

前回出てきた qDIST 系関数を使えば簡単です。

> qnorm(0.95) # 標準正規分布の95 パーセント点
[1] 1.644854
> qnorm(0.05) # 5 パーセント点
[1] -1.644854

R 備忘録 2013/05/10-2

統計学 R

ちなみに、「意味がわかる 統計解析」の流れに沿って、書いてあることが R で実行できるように R の機能を勉強しています。

確率分布

R には確率分布 DIST ごとに（DICT の累積分布関数を F(x)、確率密度関数を f(x) とおくと）

  • dDIST(x)：f(x) の値を返す。
  • pDIST(q)：F(x) の値を返す。
  • qDIST(p)：F(x) の値が p となる x を求める。
  • rDIST(n)：DIST に従って乱数を n 個生成し、長さ n のベクトルとして返す。

という４つの関数があります。
DIST にあたる部分には norm（正規分布）、unif（一様分布）、pois（ポアソン分布）などがあり、パラメータは引数で調節できます。
例えば

> dunif(0.0, -1, 1) # [-1:1] の一様分布
[1] 0.5
> pnorm(0.0)
[1] 0.5
> qnorm(0.5)
[1] 0
> runif(5)
[1] 0.58321736 0.05653621 0.51690395 0.63476821 0.69568012

となります。

乱数種（シード）の設定

set.seed を使います。

> set.seed(1)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286
> rnorm(3)
[1] 1.5952808 0.3295078 -0.8204684
> set.seed(1)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286

ちなみにベクトルを与えたら第一要素がシードになりました。

> set.seed(1:5)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286
> set.seed(5)
> rnorm(3)
[1] -0.8408555 1.3843593 -1.2554919
> set.seed(1)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286

基本統計量

mean, var sd などその名（省略名ですが）の通りの関数を使えば基本統計量を求めることができます。
なお var, sd は不偏分散、不偏標準偏差です。

> vec <- rnorm(1000)
> mean(vec)
[1] -0.02148036
> var(vec)
[1] 1.097021
> sd(vec)
[1] 1.047388
> median(vec)
[1] -0.03616252
> max(vec)
[1] 3.810277
> min(vec)
[1] -3.008049
> quantile(vec)
0% 25% 50% 75% 100%

  • 3.00804860 -0.73133329 -0.03616252 0.69697293 3.81027668

やっかいなのは度数分布の基本統計量です。

> vec_freq <- table(cut(vec, breaks=seq(-5, 5, 0.1)))

今は個表データから度数分布を生成しているので本当は元の個表データから基本統計量を計算した方が単純なのですが、そういうのは無視して、上記のような度数分布があった場合、

> weight <- seq(-5, 4.9, 0.1) + 0.05
> weighted.mean(weight, vec_freq)
[1] -0.0232

というように階級値のベクトルを作ってあげれば、weighted.mean で度数分布の平均を求めることができます。
が、そこまでするなら自分で計算してあげた方が早い気もします。

> sum(weight * vec_freq) / sum(vec_freq)
[1] -0.0232
> weight %*% vec_freq / sum(vec_freq)
[,1]
[1,] -0.0232

なお、分散等を求めるための組み込み関数はありません。
重みを考慮しながら自分で計算するか、データ数がそれほど大きくない場合には扱いを簡単にするために仮想的に個票データを作ってしまうのが良いと思います。

> (vir_vec <- do.call(c, mapply(function(f, w) return (rep(w, times=f)), as.vector(vec_freq), weight)))
[1] -3.05 -2.95 -2.95 -2.85 -2.55 -2.55 -2.55 -2.55 -2.45 -2.45 -2.45 -2.45
[13] -2.45 -2.45 -2.35 -2.35 -2.35 -2.35 -2.25 -2.25 -2.25 -2.25 -2.25 -2.25
[25] -2.15 -2.15 -2.15 -2.05 -2.05 -2.05 -2.05 -2.05 -2.05 -2.05 -1.95 -1.95
[37] -1.95 -1.95 -1.95 -1.95 -1.95 -1.95 -1.85 -1.85 -1.85 -1.85 -1.85 -1.85
[49] -1.75 -1.75 -1.75 -1.75 -1.75 -1.75 -1.75 -1.65 -1.65 -1.65 -1.65 -1.65
# 長いので省略
> mean(vir_vec)
[1] -0.0232
> var(vir_vec)
[1] 1.097319

do.call はベクトルを関数の引数に展開するための関数で、例えば do.call(f, c(1, 2, 3)) は f(1, 2, 3) と同じ意味になります。

基本統計量の要約

summary 関数を使えば、基本統計量をいくつかまとめて出力することができます。

> summary(vir_vec)
Min. 1st Qu. Median Mean 3rd Qu. Max.

  • 3.0500 -0.7500 -0.0500 -0.0232 0.6750 3.8500

summary も引数によって挙動がいろいろと変わるので総称関数っぽいですね。

R 備忘録 2013/05/10

統計学 R

データを視覚化する方法について、軽くまとめ。

とりあえずベクトルの値を並べる

plot 関数を使うとベクトルの要素をインデクス順にプロットすることができます。

> plot(c(1:5, 5:1))

ちなみにうちの環境では出力はこのようになりました。
f:id:ti2236:20130510145200p:plain
グラフを画像に保存するには png, jpeg などといったファイル形式と同名の関数を使用します。

> png(filename="r_plot.png") # デバイスを開く
> plot(c(1:5, 5:1))
> dev.off() # デバイスを閉じる

引数 type に "p" を入れると点でプロットされ、"l" を代入すると折れ線が描かれます。
"o" を入れると点と線の両方で描画されます。
他にもありますが、そんなに使わないのではないかと。

散布図を描く

２次元散布図を書くときにも plot 関数を使います。
table のように横軸成分を表すベクトルと縦軸成分を表すベクトルをそれぞれ与えてもよいですし、２列の行列（第一列が横軸成分で第二列が縦軸成分）を与えてもよいです。
以下の二行は（ラベル等を除いて）どちらも同じ結果になるはずです。

> plot(c(1:5, 5:1), c(1:10), type="o")
> plot(matrix(c(1:5, 5:1, 1:10), 10, 2), type="o")

f:id:ti2236:20130510150226p:plain

関数を描く

plot に関数を与えれば、関数をグラフ化することができます。

> plot(function(n) return (n))

横軸の下限・上限はそれぞれ第二引数、第三引数に与えます。

plot 関数は引数の構造が不思議ですね。
総称関数というやつでしょうか。

ヒストグラムを書く

hist 関数でデータの度数分布のヒストグラムを出力できます。

> # 人間の性別と身長と体重が並んだ個表データを作る。
> # 値は適当です。
> sex <- c("F", "F", "M", "F", "M", "M", "F", "M", "F", "F")
> height <- c(154, 145, 170, 167, 180, 165, 160, 177, 148, 162)
> weight <- c(48, 40, 50, 52, 80, 75, 48, 66, 58, 51)
> people <- data.frame(SEX=sex, HEIGHT=height, WEIGHT=weight)
> hist(people$HEIGHT)

データの区切り幅を変えたいときには breaks にベクトルを与えます。

# ヒストグラムの階級幅を 10 にする。
hist(people$HEIGHT, breaks=seq(130, 190, by=10))
# 身長 170 cm 以下の人と 170 cm より大きい人を数えるヒストグラムを作る。
hist(people$HEIGHT, breaks=c(130, 170, 190))

一方、度数分布表を折れ線グラフで表示するのは単純ではありません。
ベクトルの要素を先頭から順にプロットしていく plot 関数を使うために、データを加工しなければいけません。
もっと簡単な方法があれば教えてください。
なお、棒グラフを書くための barplot という関数もあります。

# 度数をプロット
> height_freq <- cut(people$HEIGHT, breaks=seq(130, 190, by=10))
> plot(table(height_freq), type="l")
# 累積度数をプロット
> plot(cumsum(table(height_freq)), type="l")
> barplot(cumsum(table(height_freq)))

ちなみに hist や plot の戻り値もオブジェクトです。

> (h <- hist(people$HEIGHT, breaks=seq(130, 190, by=10)))
$breaks
[1] 130 140 150 160 170 180 190

$counts
[1] 0 2 2 4 2 0

$intensities
[1] 0.00 0.02 0.02 0.04 0.02 0.00

$density
[1] 0.00 0.02 0.02 0.04 0.02 0.00

$mids
[1] 135 145 155 165 175 185

$xname
[1] "people$HEIGHT"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"
> mode(h)
[1] "list"

再描画はどうすればいいんでしょうか？

ラベルの付け方

plot, hist, barplot のいずれも、引数 xlab や ylab に文字列を与えることによって軸のラベルをつけることができます。

> plot(1:10, xlab="x", ylab="y")

R のデータ構造

統計学 R

とりあえず現状の理解をまとめます。

属性

R オブジェクトには属性というものを割り当てることができます。
そして、全てのオブジェクトは少なくとも mode と length の二つの本質的属性を持っています。
いくつかの属性の値は専用の関数で調べることができます。
また、それ以外の属性の値は attributes 関数で調べることができます。

> a <- 1:4
> mode(a)
[1] "numeric"
> length(a)
[1] 4
> attributes(a)
NULL
> comment(a) <- "a sample"
> attributes(a)
$comment
[1] "a sample"

> b <- matrix(1:4, 2, 2)
> mode(b)
[1] "numeric"
> length(b)
[1] 4
> attributes(b)
$dim
[1] 2 2

> c <- list(1:4)
> mode(c)
[1] "list"
> length(c)
[1] 1
> attributes(c)
NULL

> d <- data.frame(R1=1:4)
> mode(d)
[1] "list"
> length(d)
[1] 1
> attributes(d)
$names
[1] "R1"

$row.names
[1] 1 2 3 4

$class
[1] "data.frame"

> f <- function() return (1)
> mode(f)
[1] "function"
> length(f)
[1] 1
> attributes(f)
$srcref
function() return (1)

ちなみに、属性は「$属性名」で表示されていますが、ドル記号でアクセスできるわけではないようです。

属性の値を変化させたり、属性を追加したりすると、オブジェクトの振る舞いが変わります。

> is(a)
[1] "integer" "numeric" "vector"
[4] "data.frameRowLabels"
> length(a) <- 6
> a
[1] 1 2 3 4 NA NA
> mode(a) <- "complex"
> a
[1] 1+0i 2+0i 3+0i 4+0i NA NA
> is(a)
[1] "complex" "vector"
> attr(a, "dim") <- c(3, 2)
> a
[,1] [,2]
[1,] 1+0i 4+0i
[2,] 2+0i NA
[3,] 3+0i NA
> is(a)
[1] "matrix" "array" "structure" "vector"
> class(a)
[1] "matrix"
> attributes(a)
$dim
[1] 3 2

ちなみに、自分で勝手な属性を与えると、多少挙動が変わりますが本質的に大変なことは特に何もおこりません。

> attr(a, "myattribute") <- c(1, 2, 3)
> a
[,1] [,2]
[1,] 1+0i 4+0i
[2,] 2+0i NA
[3,] 3+0i NA
attr(,"myattribute")
[1] 1 2 3
> attributes(a)
$dim
[1] 3 2

$myattribute
[1] 1 2 3

以上、R オブジェクトの属性についてでした。
内部でどういう処理が行われているのか全く分かりませんが、様々なプリミティブ関数は結局この属性操作をラップしたもの、と見なしてよいのでしょうか？

ベクトルと行列と配列の関係

前の記事でベクトルや行列は配列の特殊化だと書きましたが、意味的にはともかく R の仕様的にはこれが誤りであることが上の記述から分かります。
正しくは、行列や配列はベクトルの特殊化で、dim 属性を持っているベクトルが行列・配列、です。

リストの二重カッコアクセス

[[]] があるのも mode が原因です。

> mode(c[1])
[1] "list"
> mode(c[[1]])
[1] "numeric"