2015年2月25日水曜日

Rでシミュレーション(正規分布に従う確率変数の和は正規分布に従う)

「シミュレーション」なんていうと、難しそうな気がしてしまいますが、今回のは「乱数を使って試してみる」くらいのものです。

例えば、10万人の学生がいます。彼らの身長は正規分布に従っているとします。また、彼らのテストの点数も正規分布に従っているとします。

このとき、「身長+点数」という確率変数を考えると、これは正規分布に従っているでしょうか?

これを数学的に証明するのではなく、乱数を使ったシミュレーションで、「身長+点数もきっと正規分布に従うんじゃないか」ということを確認しようというわけです。

まずは↓こんな感じでデータを生成しておきます。

# 身長
#   平均    :170cm
#   標準偏差:  5cm
heights <- rnorm(100000, 170,  5)
 
# 点数
#   平均    : 50
#   標準偏差: 10
scores  <- rnorm(100000,  50, 10)


↓身長と点数のそれぞれでヒストグラムを描いてみましょう。

brk <- seq(0,300,1) # ヒストグラムの刻み
hist(heights, breaks=brk, ylim=c(0, 10000))

正規分布に従う「身長」のヒストグラム

hist(scores,  breaks=brk, ylim=c(0, 10000))

正規分布に従う「点数」のヒストグラム

↓両者を足してヒストグラムを描いてみると、

hist(heights + scores,  breaks=brk, ylim=c(0, 10000))

「身長+点数」は正規分布に従うか?

どうやら正規分布になっているっぽいです。

では、どんな正規分布(平均、分散)になっているのでしょうか。

> mean(heights)
[1] 169.9721

> mean(scores)
[1] 50.04487

> mean(heights + scores)
[1] 220.017

「身長の平均」と「点数の平均」を足したものが、「身長+体重の平均」になっているようです。(これは当たり前の感じがしますが)

> var(heights)
[1] 24.88623

> var(scores)
[1] 99.94487

> var(heights + scores)
[1] 124.4847

分散も同様に、「身長の分散」と「点数の分散」を足したものが、「身長+体重の分散」になっているように見えます。

証明したことにはなりませんが、確率統計の教科書に出てくる

 E[X+Y] = E[X] + E[Y]
 V[X+Y] = V[X] + V[Y]

を、体感できた、というところでしょうか。

では、Rのコード例として、3つのヒストグラムをまとめて描画するものをあげておきます。

hist(heights         , breaks=brk, ylim=c(0, 10000), border="red",
     main="点数、 身長、 点数+身長 のヒストグラム", xlab="")
hist(scores          , breaks=brk, ylim=c(0, 10000), border="blue"  , add=T)
hist(heights + scores, breaks=brk, ylim=c(0, 10000), border="purple", add=T)

青(点数)と赤(身長)を足した変数のヒストグラムが紫となる

青と赤を足すと、紫になるという具合です。

ヒストグラムを見るのに慣れていないと、「青と赤を足したのに、紫はどうして両者よりも低くなっているんだ?」なんて思ってしまうかもしれません。

私も最初は混乱しましたが、ヒストグラム(=度数)を足し合わせるというのと、変数を足した上でヒストグラムを描く(度数を求める)というのは、違いますからね。

分散は足されて大きくなる
→ばらつきが大きくなる
→最頻値の階級への集中度は小さくなる
→山は小さくなる

みたいに考えると納得できるのではないでしょうか。

また、それぞれの総度数(10万人)は同じなので、赤と青と紫の面積はすべて等しいということになります。




2015年2月19日木曜日

Rでヒストグラムの上に度数の数字を表示する

例えば、↓こんなヒストグラムがあったとしましょう。

普通にhist関数で描画したもの。 30~40、40~50の階級の度数はゼロ?

ぱっと見、30~40と40~50の階級の度数はゼロなんだな、と見えちゃいますが、そうではありません。

↓実はこんなスクリプトでサンプルデータを作ったので、

d <- c(rep(5,10000), rep(15, 1000), rep(25, 100), rep(35,10), rep(45,1))
hist(d, breaks=seq(0,50,10))


30~40の階級には10の度数、40~50の階級には1の度数があるんですよね。でも、0~10の階級の10000の度数が多すぎるために、つぶれて見えなくなってしまったという状態。

こんなときは、棒グラフの上のところに度数の数字を表示させてやれば、少ない度数の階級についても確認ができます。

hist(d, breaks=seq(0,50,10), labels=T)

hist関数にlabels=Tオプションを指定し、度数を表示させたもの

labels=T オプションを指定するだけなので、とても簡単です。

が、実は私、このオプションの存在を知らず、いままで↓下記のような方法で書いてました。(このへんの紆余曲折が、この記事を書こうと思い立った理由だったりします)

h <- hist(d, breaks=seq(0,50,10))
text(h$mids, h$counts, h$counts)

hist関数の戻り値を、text関数を使ってプロット

hist関数の戻り値(変数hに格納してます)には、階級値や度数の情報が入っているので、それをtext関数で表示させるというもの。

文字と棒グラフのborderが重なって見にくいと思う場合は、↓こうやればOK。

h <- hist(d, breaks=seq(0,50,10), ylim=c(0,10100))
text(h$mids, h$counts, h$counts, pos=3)

文字が棒グラフの線と重ならないようにしたもの

text関数のオプションに pos=3 をつけてやれば、テキストの描画位置が指定座標のやや上側になります。さりげなく、y軸の範囲を広げている(ylim=c(0,10100)のところ)のは、こうしないと数字の上部が切れてしまうから。

labels=T 指定と同じようなヒストグラムが描けたものの、いろいろと面倒くさいですね。

ただ、この方法にも良い点があって、それは表示する数字を加工することができるというものです。手作りだから当たり前なのですが。

↓こんな感じ。

h <- hist(d, breaks=seq(0,50,10), ylim=c(0,10100))
text(h$mids, h$counts, paste(h$counts/1000, "千"), pos=3)

棒グラフの上に表示する数字に加工を施したもの

度数を1000で割るという演算を行い、単位として「千」という文字を追加しました。

表示する数字をいろいろとカスタマイズしたい場合は、この方法でやるとよさそうです。




Rのduplicated関数を使ってデータフレームの重複チェック(改)

duplicated関数の引数にデータフレームを指定すると、すべての項目(列)を組み合わせた上での重複チェックができますよ、という話を以前書きました↓

Rのデータフレームで複数項目を組み合わせた重複のチェック - Rプログラミングの小ネタ

また、上記の記事では、すべてではなく一部の項目のみを組み合わせ対象にしたい場合は、interaction関数を使って新たにチェック用の列を作ってはどう?、というやり方も紹介しました。

が、duplicated関数がすべての列を対象にする動きなら、対象にしたい列だけで構成された新たなデータフレームを生成して、duplicated関数に渡した方がシンプルではないか、ということに気づきました。

例えば↓こんなデータがあって、

> d
   名字 名前 年齢
1  伊東 二郎   52
2  佐藤 五郎   52
3  遠藤 太郎   53
4  後藤 太郎   53
5  近藤 五郎   28
6  遠藤 二郎   33
7  後藤 四郎   52
8  須藤 次郎   28
9  伊東 二郎   35
10 遠藤 太郎   23
11 遠藤 史郎   38
12 工藤 吾郎   53
13 須藤 太郎   22
14 武藤 二郎   28
15 工藤 吾郎   53
16 伊東 二郎   22
17 遠藤 四郎   33
18 後藤 史郎   38
19 加藤 一郎   28
20 加藤 二郎   22


名字と名前だけを対象として重複チェックをしたい場合は、

> sum(duplicated( data.frame(d$名字, d$名前) ))
[1] 4

とか、

> d[duplicated( data.frame(d$名字, d$名前) ), ]
   名字 名前 年齢
9  伊東 二郎   35
10 遠藤 太郎   23
15 工藤 吾郎   53
16 伊東 二郎   22


という感じですね。interaction関数を使うよりも、シンプルだと思います。

が、duplicated関数の仕様って、2つ目以降が現れたときにTRUEとなる、というものなんですよね。つまり、伊東二郎さんが3回現れる場合、1つ目はFALSE、2つ目はTRUE、3つ目はTRUEという結果になります。なので、sum関数で数えた場合の結果は「2」となるし。TRUEのものを表示させようとすると上記の例のように、1つ目の伊東二郎さんは表示されないことになってしまいます。

重複の有無だけを確認したい場合は上記でもよさそうですが、重複したものをすべて確認したいという場合は、interaction関数を使ってチェック用の列を作った方がいいかもしれません。↓こんな感じ。(分かりやすいように、いちいち表示させてます)

> # チェック用の「氏名」列を作る
> d$氏名 <- interaction(d$名字, d$名前, drop=T)
> d

   名字 名前 年齢      氏名
1  伊東 二郎   52 伊東.二郎
2  佐藤 五郎   52 佐藤.五郎
3  遠藤 太郎   53 遠藤.太郎
4  後藤 太郎   53 後藤.太郎
5  近藤 五郎   28 近藤.五郎
6  遠藤 二郎   33 遠藤.二郎
7  後藤 四郎   52 後藤.四郎
8  須藤 次郎   28 須藤.次郎
9  伊東 二郎   35 伊東.二郎
10 遠藤 太郎   23 遠藤.太郎
11 遠藤 史郎   38 遠藤.史郎
12 工藤 吾郎   53 工藤.吾郎
13 須藤 太郎   22 須藤.太郎
14 武藤 二郎   28 武藤.二郎
15 工藤 吾郎   53 工藤.吾郎
16 伊東 二郎   22 伊東.二郎
17 遠藤 四郎   33 遠藤.四郎
18 後藤 史郎   38 後藤.史郎
19 加藤 一郎   28 加藤.一郎
20 加藤 二郎   22 加藤.二郎


> # 重複している行を抽出
> dup <- d[duplicated(d$氏名), ]
> dup

   名字 名前 年齢      氏名
9  伊東 二郎   35 伊東.二郎
10 遠藤 太郎   23 遠藤.太郎
15 工藤 吾郎   53 工藤.吾郎
16 伊東 二郎   22 伊東.二郎


> # 2回以上登場(=元のデータに3回以上登場)への対策
> dup_u <- unique(dup$氏名)
> dup_u

[1] 伊東.二郎 遠藤.太郎 工藤.吾郎
16 Levels: 加藤.一郎 近藤.五郎 佐藤.五郎 工藤.吾郎 遠藤.史郎 後藤.史郎 ... 武藤.二郎


> # ダブりのあったものを一覧として表示
> for(i in 1:length(dup_u)){
+   print( d[which(d$氏名 == dup_u[i]), ] )
+   cat("---------------------------\n")
+ }

   名字 名前 年齢      氏名
1  伊東 二郎   52 伊東.二郎
9  伊東 二郎   35 伊東.二郎
16 伊東 二郎   22 伊東.二郎
---------------------------
   名字 名前 年齢      氏名
3  遠藤 太郎   53 遠藤.太郎
10 遠藤 太郎   23 遠藤.太郎
---------------------------
   名字 名前 年齢      氏名
12 工藤 吾郎   53 工藤.吾郎
15 工藤 吾郎   53 工藤.吾郎
---------------------------


これで、重複のあった対象の1つ目のものも確認できるようになりました。




2015年2月10日火曜日

Rで複数のヒストグラムを比較する方法あれこれ

2つ以上のヒストグラムを1つのグラフ上に描画して、比較したいようなときってありますよね。

例えば、集団aのテストの点数と、集団bのテストの点数が、それぞれ1000件ずつあるとしましょう。

# サンプルデータを作っておく
a <- rnorm(1000, 40, 10)
b <- rnorm(1000, 60, 10)


borderオプションを使って、ヒストグラムの線に色をつけるってのは1つのやり方ですね↓

# 線に色をつける
hist(a, breaks=seq(0,100,5), border="red", xlim=c(0,100), main="", xlab="")
hist(b, breaks=seq(0,100,5), border="blue", add=T)

borderオプションでヒストグラムの線に色をつける

なんとなく見えてはいるのですが、重なった縦線は後から描画した方で上書きされてしまい、そのへんがちょっと見にくいですね。

colオプションで面に色をつけることもできますが・・・

# 色で塗りつぶす
hist(a, breaks=seq(0,100,5), col="red" , xlim=c(0,100), main="", xlab="")
hist(b, breaks=seq(0,100,5), col="blue", add=T)

colオプションでヒストグラムを塗りつぶす

↑普通の色で塗っちゃうと重なったところが見えなくなってしまいます。

そういうときは、透過色を使えばいいです↓

# 透過色を使う
hist(a, breaks=seq(0,100,5), col="#FF00007F", xlim=c(0,100), main="", xlab="" )
hist(b, breaks=seq(0,100,5), col="#0000FF7F", add=T)
colオプションに透過色を指定しヒストグラムを塗りつぶす

colオプションへの指定は、「##rrggbbaa」という感じで2桁ずつ、「赤・緑・青・不透明度」を16進で指定します。上記の例では、赤の50%透過と青の50%透過を使っているような感じですので、重なったところが紫になります。

赤と青と紫の3つのヒストグラムに見えるというあなたは、ちょっぴり偏屈かもしれません。

3つ以上重なってくると、さすがに上記のようなやり方は無理があるので、いっそのこと一本の線にしてしまうのがいいように思います。

実はhist関数には戻り値があって、階級値やら度数やらの情報が格納されています。それを使って、線でプロットする(plotのtype="l"指定とか、lines関数での描画とか)という方法があります↓

#
# ヒストグラムではなく線グラフで表現
#
 
# hist関数を使って階級値ごとの度数を求めておく(描画結果は不要)
ha <- hist(a, breaks=seq(0,100,5))
hb <- hist(b, breaks=seq(0,100,5))
 
# mids(=階級値)とcounts(=度数)の情報を使って描画する
plot(0, 0, type="n", xlim=c(0,100), ylim=c(0, 200), xlab="", ylab="")
lines(ha$mids, ha$counts, col="red")
lines(hb$mids, hb$counts, col="blue")
plot関数のtype="l"オプション または lines関数で描画する

縦軸のスケールが度数ではなく確率密度になっても構わないなら、density関数による密度推定を使うという方法もあります↓

# density関数を使ったカーネル密度推定
plot(0, 0, type="n", xlim=c(0,100), ylim=c(0, 0.04), xlab="", ylab="")
lines(density(a), col="red")
lines(density(b), col="blue")

density関数を使ってカーネル密度推定結果で描画する

度数をそのままプロットするより、滑らかな感じになりました。おおよその分布の傾向を見ることが目的なら、こちらの方法の方が簡単な気がします。

また、度数の差が大きいものどうし(100人の得点データと1000人の得点データなど)の場合でも、確率密度を使うと同じ土俵で比較することが可能になります。




2015年2月7日土曜日

Rのinteraction関数で重複チェックするときは、drop=Tを指定する

前回の記事で、データフレームの重複をチェックするときに、interaction関数を使うという方法を紹介しました↓

Rのデータフレームで複数項目を組み合わせた重複のチェック


例えば、こんな↓データだったら、

> d
   名字 名前 年齢 出身         職業
1  佐藤 一郎   20 埼玉   アルバイト
2  武藤 二郎   22 大阪         学生
3  加藤 三郎   25 愛知       会社員
4  伊東 四郎   77 東京 コメディアン
5  近藤 五郎   32 千葉         無職
6  後藤 太郎   33 兵庫       公務員
7  工藤 次郎   35 福岡     専業主夫
8  須藤 三朗   38 静岡       自営業
9  遠藤 史郎   52 茨城     会社役員
10 斉藤 吾郎   83 広島     定年退職


↓こんなふうにして、

> intr <- interaction(d$名字, d$名前, d$年齢, d$出身, d$職業, drop=T)
> sum(duplicated(intr))

[1] 0

TRUEの数がゼロだから重複なしだな、みたいな。

で、前回の記事でも上記のサンプルでも、interaction関数にさりげなく、drop=Tオプションを指定していますが、これが結構重要です。何も指定していないと、デフォルトはdrop=Fです。

両者で実行結果を比較してみましょう。

↓まずはdrop=Tで、

> interaction(d$名字, d$名前, d$年齢, d$出身, d$職業, drop=T)
 [1] 佐藤.一郎.20.埼玉.アルバイト   武藤.二郎.22.大阪.学生       
 [3] 加藤.三郎.25.愛知.会社員       伊東.四郎.77.東京.コメディアン
 [5] 近藤.五郎.32.千葉.無職         後藤.太郎.33.兵庫.公務員     
 [7] 工藤.次郎.35.福岡.専業主夫     須藤.三朗.38.静岡.自営業     
 [9] 遠藤.史郎.52.茨城.会社役員     斉藤.吾郎.83.広島.定年退職   
10 Levels: 佐藤.一郎.20.埼玉.アルバイト ... 近藤.五郎.32.千葉.無職


↓次はdrop=Fで、

> interaction(d$名字, d$名前, d$年齢, d$出身, d$職業, drop=F)
 [1] 佐藤.一郎.20.埼玉.アルバイト   武藤.二郎.22.大阪.学生       
 [3] 加藤.三郎.25.愛知.会社員       伊東.四郎.77.東京.コメディアン
 [5] 近藤.五郎.32.千葉.無職         後藤.太郎.33.兵庫.公務員     
 [7] 工藤.次郎.35.福岡.専業主夫     須藤.三朗.38.静岡.自営業     
 [9] 遠藤.史郎.52.茨城.会社役員     斉藤.吾郎.83.広島.定年退職   
100000 Levels: 伊東.一郎.20.愛知.アルバイト ... 武藤.二郎.83.兵庫.無職


大して変わらないじゃないかと言われそうですが、実行時間が大きく違います。前者は一瞬で結果が出るのですが、後者の実行には私の環境では3秒くらいかかりました。

3秒くらいだったら我慢できるでしょうが、これがあと2,3列増えるだけで、作業するには現実的でないような実行時間に陥ってしまいます。

drop=Fを指定すると、データ上には存在しないものも含めて、全ての組み合わせを水準として保持します。

100000 Levels: 伊東.一郎.20.愛知.アルバイト ... 武藤.二郎.83.兵庫.無職

とあるように、名字が10種類、名前が10種類、年齢が10種類、出身が10種類、職業が10種類あるので、10×10×10×10×10 = 10万 の水準が作られます。これが時間を食う原因になっているようです。

目的が重複のチェックだけだったら、存在しない組み合わせの水準は必要ないので、drop=Fと指定すれば、大抵の場合は一瞬で組み合わせ因子が作れるはずです。




2015年2月1日日曜日

Rのデータフレームで複数項目を組み合わせた重複のチェック

例えば、↓こんな感じのデータフレームがあったとします。

> d
   名字 名前 年齢
1  伊東 二郎   52
2  佐藤 五郎   52
3  遠藤 太郎   53
4  後藤 太郎   53
5  近藤 五郎   28
6  遠藤 二郎   33
7  後藤 四郎   52
8  須藤 次郎   28
伊東 二郎   35
10 遠藤 太郎   23
11 遠藤 史郎   38
12 工藤 吾郎   53
13 須藤 太郎   22
14 武藤 二郎   28
15 工藤 吾郎   53
16 伊東 二郎   22
17 遠藤 四郎   33
18 後藤 史郎   38
19 加藤 一郎   28
20 加藤 二郎   22

仮に、名字と名前と年齢の一致によって個人を特定できるとし、このデータに重複がないかチェックをしたいということにします。

全項目を加味した重複チェックを行うのは実は簡単で、duplicated関数の引数にデータフレームを渡してやればOKです。TRUEがあれば重複があるということになります↓

> duplicated(d) 
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[11] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE


重複の有無を一目で確認したいときは、sum関数を使います。1以上なら重複ありということになります↓

> sum(duplicated(d))
[1] 1

重複したのがどの行だったのか見たい場合は、↓こんな感じでしょうか。

> d[duplicated(d), ]
   名字 名前 年齢
15 工藤 吾郎   53

データフレームの全項目ではなく、一部の項目の組み合わせでの重複チェックをしたい場合は、interaction関数を使ってそれらの項目を組み合わせた因子を作って、それをIDのように扱うという方法もあります。

年齢の項目は無視して、同姓同名を見つけたいような場合は、名字と名前の2つを組み合わせた因子の項目を作ります。

> d$ID <- interaction(d$名字, d$名前, drop=T)
> d
   名字 名前 年齢        ID
伊東 二郎   52 伊東.二郎
2  佐藤 五郎   52 佐藤.五郎
3  遠藤 太郎   53 遠藤.太郎
4  後藤 太郎   53 後藤.太郎
5  近藤 五郎   28 近藤.五郎
6  遠藤 二郎   33 遠藤.二郎
7  後藤 四郎   52 後藤.四郎
8  須藤 次郎   28 須藤.次郎
伊東 二郎   35 伊東.二郎
10 遠藤 太郎   23 遠藤.太郎
11 遠藤 史郎   38 遠藤.史郎
12 工藤 吾郎   53 工藤.吾郎
13 須藤 太郎   22 須藤.太郎
14 武藤 二郎   28 武藤.二郎
15 工藤 吾郎   53 工藤.吾郎
16 伊東 二郎   22 伊東.二郎
17 遠藤 四郎   33 遠藤.四郎
18 後藤 史郎   38 後藤.史郎
19 加藤 一郎   28 加藤.一郎
20 加藤 二郎   22 加藤.二郎

> duplicated(d$ID)
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
[11] FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE


> sum(duplicated(d$ID))
[1] 4

ちなみに、サンプルのデータはsample関数を使って、組み合わせをランダムに生成して作ったのですが、「伊東四郎」は現れませんでした。■

[2015年2月7日追記]
 サンプルコードにつけている「drop=T」は結構重要だという話↓
Rのinteraction関数で重複チェックするときは、drop=Tを指定する - Rプログラミングの小ネタ

■ 

[2015年2月19日追記]
重複の有無だけでなく、どれが重複したのかちゃんと確認するところまで、サンプルを作ってみました↓

Rのduplicated関数を使ってデータフレームの重複チェック(改) - Rプログラミングの小ネタ