計量政治分析で便利なRの操作
Masahiko
December 10, 2018
1. Rの基本操作
・Rを動かすためには、コマンドを入力する必要がある。
・ここでは基本的な四則演算等の使い方を紹介する。
四則演算
1 + 1 ## [1] 210 - 2## [1] 8333 * 3## [1] 999333 / 3## [1] 111べき乗
5 ^ 2## [1] 25平方根
sqrt(25)## [1] 5連続した整数を表示（1から10まで）
1:10## [1] 1 2 3 4 5 6 7 8 9 10連続した整数の平均値を表示（1から10までの平均値）
mean(1:10)## [1] 5.5(1:10)・・・引数（引数）
引数とは:
・例えば、スターバックスに行ったとき「キャラメルマキアートひとつ！」と注文して「お金を払う」と、実際にキャラメルマキアートがひとつ渡される。
・この場合「キャラメルマキアートひとつ！」と注文することと「お金を払うこと」が引数。
・「キャラメルマキアートを渡されること」が結果の表示。
・「引数」の定義：「数学における関数やコンピュータプログラムにおける手続きにおいて、外部と値をやりとりするための変数もしくはその値」
・引数は複数指定することもできる。
seq(from = 0, to = 100, by = 20)## [1] 0 20 40 60 80 1002. パイプ演算子：%>%
・R でデータの操作をするときに非常に便利なツール
・RStudio では、 Command(Control) + Shift + m キーで入力できる
・magrittr というパッケージに含まれる演算子
・dplyr を読み込むと magrittr も読み込まれる
library("dplyr")・演算子の左側の処理結果を、演算子右側の関数の第 1 引数として利用する
(40 - 4) %>% sqrt()## [1] 6・まず 40 - 4 を計算し、その値 36 を引数として、右側にあるsqrt() 関数 の第 1 引数として利用する => つまり sqrt(36) = 6 を計算
40 - 4 %>% sqrt()## [1] 38・4 を右側にある sqrt()関数の第 1 引数として計算 => つまり sqrt(4) = 2 を計算
・40 から 2 を引く = 38
パイプ演算子の本領は、複数のデータ操作処理を連続して行うときに発揮される
まず、衆院選挙データ (HR)をダウンロードし HR というデータフレーム名をつける。
download.file(url = "https://git.io/fAnI2", destfile = "hr96-17.csv")library("readr") HR <- read_csv("hr96-17.csv", na = ".")HRは1996年から2017年までに行われた8回の衆院選挙結果データである。
summary()を使って、データの変数の様子を表示してみる。
summary(HR) year ku kun status
Min. :1996 Length:8803 Min. : 1.000 Min. :0.0000
1st Qu.:2000 Class :character 1st Qu.: 2.000 1st Qu.:0.0000
Median :2005 Mode :character Median : 4.000 Median :0.0000
Mean :2007 Mean : 5.738 Mean :0.4857
3rd Qu.:2012 3rd Qu.: 8.000 3rd Qu.:1.0000
Max. :2017 Max. :25.000 Max. :2.0000
name party party_code previous
Length:8803 Length:8803 Min. : 1.00 Min. : 0.000
Class :character Class :character 1st Qu.: 1.00 1st Qu.: 0.000
Mode :character Mode :character Median : 3.00 Median : 0.000
Mean : 12.04 Mean : 1.712
3rd Qu.: 8.00 3rd Qu.: 3.000
Max. :100.00 Max. :20.000
wl voteshare age nocand
Min. :0.0000 Min. : 0.10 Min. :25.0 Min. :2.000
1st Qu.:0.0000 1st Qu.: 8.90 1st Qu.:43.0 1st Qu.:3.000
Median :0.0000 Median :25.76 Median :51.0 Median :4.000
Mean :0.4659 Mean :27.08 Mean :50.9 Mean :3.956
3rd Qu.:1.0000 3rd Qu.:42.90 3rd Qu.:59.0 3rd Qu.:5.000
Max. :2.0000 Max. :95.30 Max. :94.0 Max. :9.000
NA's :5
rank vote eligible turnout
Min. :1.000 Min. : 177 Min. :115013 Min. :48.90
1st Qu.:1.000 1st Qu.: 18240 1st Qu.:269294 1st Qu.:57.80
Median :2.000 Median : 49021 Median :330188 Median :62.70
Mean :2.477 Mean : 54911 Mean :325666 Mean :62.96
3rd Qu.:3.000 3rd Qu.: 86494 3rd Qu.:390637 3rd Qu.:67.53
Max. :9.000 Max. :201461 Max. :495212 Max. :83.80
NA's :959 NA's :1895
exp
Min. : 535
1st Qu.: 2803566
Median : 6541589
Mean : 7542700
3rd Qu.:11044485
Max. :27462362
NA's :2043 ここで例えば「2009 年の選挙だけ抜き出す」ためには filter()関数を使えばよい。
HR2009 <- filter(HR, year == 2009)summary(HR2009$year) Min. 1st Qu. Median Mean 3rd Qu. Max.
2009 2009 2009 2009 2009 2009 2009年のデータだけが抜き出されていることがわかる。
さらに、次のことも同時に行いたいとする。
- 小選挙区での得票率順に並べ替え(arrange)、第 1 位から最下位まで 順位 order をつける(mutate)
- 年齢が 40 歳以下の候補を除外する(filter)
- year, ku、kun、party、name、order だけ残す(select)
このような時、次のようにパイプ演算子を使うと、間違いが少なく効率的に実行できる。
HR09 <- HR %>%
filter(year == 2009) %>% #2009のデータだけ抜き出す
arrange(desc(voteshare)) %>% #得票率が大きい順に並べる
mutate(order = 1:n()) %>% #候補者に順位を付ける
arrange(ku, kun) %>% #選挙区順に並べかえる
filter(age > 40) %>% #候補者の年齢を40歳以上に指定
select(year, ku, kun, party, name, voteshare, order) #変数を指定 head(HR09)# A tibble: 6 x 7
year ku kun party name voteshare order
<int> <chr> <int> <chr> <chr> <dbl> <int>
1 2009 aichi 1 DPJ SATO, YUKO 54.4 130
2 2009 aichi 1 JCP KIMURA, EMI 6.4 744
3 2009 aichi 1 SDP HIRAYAMA, RYOHEI 2.7 838
4 2009 aichi 1 kofuku KAWADA, SEIJI 1.5 952
5 2009 aichi 2 DPJ FURUKAWA, MOTOHISA 66.6 9
6 2009 aichi 2 LDP MIYAHARA, MISAKO 23.9 601tail(HR09)# A tibble: 6 x 7
year ku kun party name voteshare order
<int> <chr> <int> <chr> <chr> <dbl> <int>
1 2009 yamanashi 2 independent NAGASAKI, KOTARO 32.1 549
2 2009 yamanashi 2 LDP HORIUCHI, MITSUO 29.6 572
3 2009 yamanashi 2 kofuku MIYAMATSU, HIROYUKI 0.7 1117
4 2009 yamanashi 3 DPJ GOTO, HITOSHI 62.7 20
5 2009 yamanashi 3 LDP ONO, JIRO 35.3 502
6 2009 yamanashi 3 kofuku SAKURADA, DAISUKE 2 900・選挙区ごと、得票率が大きな順に並べ替えられていることがわかる。
・order の番号を見れば、それぞれの候補者が2009年の総選挙に出馬した全ての候補者の中で得票率のランクがわかる。
3. 複数の中から無作為に選ぶ
連続変数の中から数字を無作為に選ぶ
・1から100までの番号から無作為に数字を抽出したければ、次のようにする。
sample(1:100, 3)[1] 88 60 15アルファベットを無作為に選ぶ
LETTERS [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q"
[18] "R" "S" "T" "U" "V" "W" "X" "Y" "Z"複数のアルファベットを抽出したいときは、番号をベクタとして入力します。
LETTERS[c(1, 2, 3)][1] "A" "B" "C"・三つのアルファベットを無作為に抽出したければ、次のようにする。
LETTERS[sample(1:26, 3)][1] "J" "U" "Z"・10名のゼミ学生の中から一人を選ぶとする。
・まず、ゼミ学生の名簿 members を作成する。
members <- c("吉川", "名村", "近藤", "斉藤", "菅原", "江馬", "大川", "尾崎", "滝瀬", "藤枝")・sample()関数を使って、上で作成した members から一人だけ選びたければ、次のようにする。
sample(members, size = 1)[1] "尾崎"・一気に5人選びたければ、次のコマンドを使う。
sample(members, size = 5)[1] "滝瀬" "吉川" "江馬" "菅原" "尾崎"・ここでは、10人から「別々の学生を5人」選んでいる。
・次のような選び方で5人を選びたければ、replace = TRUE を指定すればよい。
一回目：10人から一人選ぶ
二回目：10人から一人選ぶ
三回目：10人から一人選ぶ
四回目：10人から一人選ぶ
五回目：10人から一人選ぶ
sample(members, size = 5, replace = TRUE)[1] "斉藤" "菅原" "菅原" "江馬" "江馬"4. 衆院選挙データを使った無作為抽出例
・2016年から2017年までの衆院選データから2017年の衆院選データだけを抜き出す。
hr2017 <- HR %>%
filter(year == 2017) %>% #2009のデータだけ抜き出す
arrange(ku, kun, rank) %>% #選挙区、当選ランクごとに並べ替える
mutate(order = 1:n()) %>% #候補者にシリアル番号を付ける
select(order, ku, kun, party, age, name, status, wl) #必要な変数を指定summary(hr2017) order ku kun party
Min. : 1.0 Length:936 Min. : 1.000 Length:936
1st Qu.:234.8 Class :character 1st Qu.: 2.000 Class :character
Median :468.5 Mode :character Median : 4.000 Mode :character
Mean :468.5 Mean : 5.755
3rd Qu.:702.2 3rd Qu.: 8.000
Max. :936.0 Max. :25.000
age name status wl
Min. :25.00 Length:936 Min. :0.0000 Min. :0.0000
1st Qu.:44.00 Class :character 1st Qu.:0.0000 1st Qu.:0.0000
Median :53.00 Mode :character Median :1.0000 Median :0.0000
Mean :52.77 Mean :0.5929 Mean :0.5556
3rd Qu.:62.00 3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :80.00 Max. :2.0000 Max. :2.0000 head(hr2017)# A tibble: 6 x 8
order ku kun party age name status wl
<int> <chr> <int> <chr> <int> <chr> <int> <int>
1 1 aichi 1 LDP 53 KUMADA, HIROMICHI 1 1
2 2 aichi 1 CDP 42 YOSHIDA, TSUNEHIKO 2 2
3 3 aichi 1 POH 54 SATO, YUKO 2 0
4 4 aichi 2 POH 51 FURUKAWA, MOTOHISA 1 1
5 5 aichi 2 LDP 45 TABATA, TSUYOSHI 1 0
6 6 aichi 2 JCP 31 SAKAI, KENTARO 0 0・例えば、データフレーム (hr2017)に納められている 936名の候補者から無作為に10人を選びたければ、次のようにする。
library("dplyr")sample_n(tbl = hr2017, size = 10)# A tibble: 10 x 8
order ku kun party age name status wl
<int> <chr> <int> <chr> <int> <chr> <int> <int>
1 703 saitama 11 independent 42 KONNO, TOMOHIRO 0 0
2 211 gunma 5 JCP 34 ITO, TATSUYA 0 0
3 502 miyazaki 2 others 57 KONO, ICHIRO 1 0
4 647 osaka 15 NIK 44 URANO, YASUTO 0 2
5 145 fukuoka 5 LDP 73 HARADA, YOSHIAKI 1 1
6 517 nagano 3 others 57 OIKAWA, YUKIHISA 1 0
7 131 fukui 2 POH 43 SAIKI, TAKESHI 2 2
8 274 hyogo 1 JCP 39 RIKISHIGE, TOMOYUKI 0 0
9 330 ibaraki 7 independent 68 NAKAMURA, KISHIRO 1 1
10 226 hiroshima 4 LDP 42 SHINTANI, MASAYOSHI 1 120代から40代を younger、50代以上を senior と分類し、 それぞれ2つのグループから5人ずつ選ぶとすれば、次のようにする。
younger <- hr2017 %>%
filter(age < 50) %>% #20代、30代、40代
sample_n(size = 5)
younger# A tibble: 5 x 8
order ku kun party age name status wl
<int> <chr> <int> <chr> <int> <chr> <int> <int>
1 274 hyogo 1 JCP 39 RIKISHIGE, TOMOYUKI 0 0
2 745 shizuoka 3 LDP 42 MIYAZAWA, HIROYUKI 1 1
3 164 fukuoka 11 POH 48 MURAKAMI, TOMONOBU 0 0
4 314 ibaraki 1 others 30 KAWABE, KENICHI 0 0
5 139 fukuoka 2 JCP 45 MATSUO, RITSUKO 0 0senior <- hr2017 %>%
filter(age > 49) %>% #50代以上
sample_n(size = 5)
senior# A tibble: 5 x 8
order ku kun party age name status wl
<int> <chr> <int> <chr> <int> <chr> <int> <int>
1 215 hiroshima 2 LDP 69 HIRAGUCHI, HIROSHI 1 1
2 558 niigata 5 independent 61 ODAIRA, ETSUKO 0 0
3 880 tokyo 25 CDP 58 YAMASHITA, YOKO 1 0
4 225 hiroshima 3 independent 72 NISHIMOTO, AKIHIKO 0 0
5 317 ibaraki 2 JCP 69 HOSHINO, FUMIO 0 05. データを作って図を描く
・-10 から 10 まで、1 区切りで整数の値を並べてベクトルを作り、x という変数名をつける。
x = seq(from = -10, to = 10, by = 1) # ベクトルに x という変数名をつける・y を x の累乗に指定する。
y = x^2・二つの変数をプロットする。
・type = “l” を加えて、直線で表す
plot(x, y, col = "tomato", type = "l")
・三次関数は次の様にして表示できる。
x = seq(from = -10, to = 10, by = 1) # ベクトルに x という変数名をつける
y = x^3
plot(x, y, col = "blue", type = "l")
6. 算術演算子の優先順位
優先順位：累乗＞乗算＞加算
1 + 2 * 2 # 2 * 2 = 4 → 1 * 4[1] 5(1 + 2) * 2 # (1 + 2) = 3 → 3 * 2[1] 6(1 + 2) ^ 2 # (1 + 2) = 3 → 3 ^ 2[1] 9((1 + 2) * 2) ^ 2 #(1 + 2) = 3 → 3 * 2 = 6 → 6 ^ 6[1] 367. csvファイルの読み込みと factor
library("readr")
hr = read.csv("hr_sample.csv")・データを読み込んだ hr はデータフレーム型。
head(hr) year ku kun party name age status wl vote
1 2014 tokyo 10 LDP KOIKE, YURIKO 62 incumbent 1 93610
2 2014 tokyo 10 DPJ EBATA, TAKAKO 54 moto 0 44123
3 2014 tokyo 10 JCP KON, HIDEKO 66 challenger 0 28453
4 2014 tokyo 10 seikatsu TAGAYA, RYO 46 challenger 0 9663
5 2014 tokyo 10 NG KAMITANI, TIZUKO 61 challenger 0 8688
voteshare
1 50.7
2 23.9
3 15.4
4 5.2
5 4.7・hr の変数 (year, ku, kun,…など)はベクトル型、もしくはファクター型。
・数値以外の変数は全てファクター型に変化される。
・party の型をチェック。
hr$party[1] LDP DPJ JCP seikatsu NG
Levels: DPJ JCP LDP NG seikatsu・5つのレベルのファクター型に変換されている。
・status をチェック。
hr$status[1] incumbent moto challenger challenger challenger
Levels: challenger incumbent moto・status は 3 つのレベルのファクター型に変換されている。
・status の水準 (level) をチェック。
as.numeric(hr$status)[1] 2 3 1 1 1・status の level は、デフォルトでアルファベット順: challenger = 1, incumbent = 2, moto = 3
・順番を自在に指定できる。
・例えば、: incumbent = 1, challenger = 2, moto = 3 と指定したければ次の様にする。
hr$status = factor(hr$status,
level = c("incumbent", "challenger", "moto"))
hr$status[1] incumbent moto challenger challenger challenger
Levels: incumbent challenger motoas.numeric(hr$status)[1] 1 3 2 2 2ファクター型をベクトルとして使いたい場合
・name の型をチェック。
hr$name[1] KOIKE, YURIKO EBATA, TAKAKO KON, HIDEKO TAGAYA, RYO
[5] KAMITANI, TIZUKO
5 Levels: EBATA, TAKAKO KAMITANI, TIZUKO KOIKE, YURIKO ... TAGAYA, RYO・nameも5つのレベルのファクター型に変換されている。
・name が重複することは希なので、nameをファクター型として使うメリットはない。
・as.vector()関数を使って、name を元のベクトルに変換する。
hr$name = as.vector(hr$name)
hr$name[1] "KOIKE, YURIKO" "EBATA, TAKAKO" "KON, HIDEKO"
[4] "TAGAYA, RYO" "KAMITANI, TIZUKO"数値情報をファクター型として使いたい場合
・wl をチェック。
hr$wl[1] 1 0 0 0 0・数値ベクトルとして読み込まれている。
・faactor()関数を使って、ファクターに変換。
hr$wl = factor(hr$wl)
hr$wl[1] 1 0 0 0 0
Levels: 0 1・wl が2つのレベルのファクターに変換されたことがわかる。
8. 架空コインを作る
・人工的にコインを作成し、500回コイントスをして、「表」が出る割合をシミュレーションしてみる。
N = 500
probHeads = 0.5 # 偏りのないコインを想定：「表」が出る確率 = 0.5
flip = sample(x = c(0,1),
prob = c(1-probHeads, probHeads),
size = N,
replace = TRUE)
r = cumsum(flip) # 累積和関数を使って、トスの回数を重ねる度に割合を update
n = 1:N
runProp = r / n
plot(n, runProp, type = "o", log = "x", col = "blue",
xlim = c(1, N), ylim = c(0.0, 1.0), cex.axis = 1,
xlab = "コイントスの回数", ylab = "「表」が出る割合", cex.lab = 1,
main = "コイントスで「表」が出る割合のシミュレーション", cex.main = 1)
abline(h = probHeads, lty = "dotted", col = "tomato") # 「表」が出る割合 = 0.5 に点線を引く
flipLetters = paste(c("裏", "表")[flip[1:10]+1], collapse = "")
displayString = paste0("コイン投げの結果 = ", flipLetters, "...")
text(N, .9, displayString, adj = c(1, 0.5), cex = 1)
text(N, .8, paste("最終的に「表」が出る割合 = ", runProp[N]), adj = c(1, 0.5), cex = 1)
・コイントスを試す度に結果が異なる。
・cumsum()関数を使って、累積和を計算している。
・1回目に「表」が出て、2回目に「表」が出た場合・・・「表」が出る割合 = 1
・1回目に「表」が出て、2回目に「裏」が出た場合・・・「表」が出る割合 = 0.5
・1回目に「裏」が出て、2回目に「裏」が出た場合・・・「表」が出る割合 = 0
・結果を見やすくするため、x 軸の「コイントスの回数」は log で表示。
・コイントスの回数が増えるにつれて、累積割合が 0.5 に近づいている。
9. 架空のサイコロを作る
架空のサイコロを作り die と名前を付ける。
die <- 1:6
die## [1] 1 2 3 4 5 6次に、作ったサイコロを三回振ってみる。
sample(die, size = 1)[1] 5sample(die, size = 1)[1] 1sample(die, size = 1)[1] 1サイコロを振る度に異なる値が出ることわかる。
次に、架空のサイコロ die を二回振って出た目をdice と名前を付ける。
dice <- sample(die, size =2)
dice[1] 6 4dice と続けて二回入力しても、同じサイコロの二つ目が出る。
dice[1] 6 4dice[1] 6 4サイコロを振る度に異なる結果を出すためには、独自関数 function を使い、次の様に入力する。サイコロを二回振って出た目の合計を計算する。
roll <- function(){
die <- 1:6
dice <- sample(die, size = 2, replace = TRUE)
sum(dice)
}続けて roll() と二回入力すると、異なるサイコロの目の総和が出る（はず）。
roll()[1] 5roll()[1] 9roll() から() を除いて roll だけ入力すると、関数内に格納されているコード（＝関数の本体）だけを表示する。
rollfunction(){
die <- 1:6
dice <- sample(die, size = 2, replace = TRUE)
sum(dice)
}
<bytecode: 0x7fcd53a7a058>ここで二つのサイコロを同時に10回振って、出た目の総和を確かめてみる。
rolls <- replicate(10, roll())
rolls [1] 8 8 4 2 7 5 5 9 9 7この結果をヒストグラムで表すと次の様になる。
library(ggplot2)qplot(rolls, binwidth = 1)
二つのサイコロを同時に振る回数を10000回に増やして、出た目の総和をヒストグラムにしてみると下図のようなきれいな山の形になる。
rolls <- replicate(10000, roll())
qplot(rolls, binwidth = 1)
10. ウェイトを付けたサイコロを作る
ここで、二個のサイコロの出目がそれぞれ1/6の当確率ではなく、下図のように、2 から 6 までの目が 1/9 の確率、1 だけが 4/9 の確率で出るように細工をしたサイコロを人工的に作る。
roll <- function(){
die <- 1:6
dice <- sample(die, size = 2, replace = TRUE,
prob = c(4/9, 1/9, 1/9, 1/9, 1/9, 1/9))
sum(dice)
}この二個のサイコロを 10000 回投げてその総和をヒストグラムにすると、ウェイトをかけた 1 の出目だけが突出していることがわかる。
rolls <- replicate(10000, roll())
qplot(rolls, binwidth = 1)
11. 単純な条件文
・オブジェクト x （ここでは数字の 9）が「 10 未満の数かどうか」を確かめるプログラムは次のように書くことができる。
x <- 9
if(x < 10) print("yes") else print("no")[1] "yes"つまり、数字の 9 は 10 未満なので “yes” と判定される。
print というコマンドを省略して次のようにも書くことができる。
x <- 9
if(x < 10) {
"yes"
} else {
"no"
}[1] "yes"・\(y = x^2 + a\) という二次関数をプログラムしてみる。
xsqplusa = function(x, a = 5){
y = x^2 + a
return(y)
}・xsqplusa と名付けられた関数は、x を二乗し、a を加え、その後 y として結果を出力。
・x と a にそれぞれ値を指定して出力してみる。
xsqplusa(x = 2, a = 5)[1] 9・関数の呼び出しの際、明示的に引数を指定する場合は、順番はどちらでもいい。
xsqplusa(a = 5, x = 2)[1] 9・引数 x はデフォルトの値をもたないが、引数 a のデフォルトは 1。
xsqplusa(x = 3) # a を指定しなければ、引数 a = 1 のデフォルト値が使われる [1] 1412. 複数の判定条件を並べる場合
また、次のように、判定条件を複数並べることもできる。
ここでは統計の仮説検定において 0.06 というp-値（有意確率）が何%の有意水準で棄却されるかプログラムすると。
p <- 0.06
if (p < 0.01){
"1% significance"
} else if (p < 0.05){
"5% significance"
} else if (p < 0.1){
"10% significance"
} else if (p >= 0.1){
"not significant"
}[1] "10% significance"2.5. ベクトル全体に条件式をあてはめる場合
ベクトル全体に条件式をあてはめて実行する場合には ifelse( ) 関数を使う。
例えば、1〜5 の整数で構成されたベクトル x 一つ一つに「3 より大きいかどうか」という条件式を当てはめる場合は次のように書く。
x <- 1:5
ifelse(x > 3, "yes", "no")[1] "no" "no" "no" "yes" "yes"x <- 1:5
y <- ifelse (x > 3,
{cat("yes"); print(x / 10) },
{cat( "no"); print(x * 10) }
)yes[1] 0.1 0.2 0.3 0.4 0.5
no[1] 10 20 30 40 502.6. switch( )関数
条件として指定された値ごとに異なる値を返す関数。
選択肢が複数ある場合に使う。else if を書くこともできるが煩雑。
switch は第1引数としてテストしたい値（この場合は x）をとる。
univ <- function(x)
{
switch(x,
"t" = "takushoku",
"w" = "waseda",
"k" = "kobe",
"other")
}“t”, “w”, “k” という第1引数を入力してテストしてみると、指定した大学名が表示され、これら以外の文字（ここでは“a”）に対しては “other” が表示される。
univ("t")[1] "takushoku"univ("w")[1] "waseda"univ("k")[1] "kobe"univ("a")[1] "other"ここで第1引数に文字でなく、" " を外して数値を入力するとswitch()関数で指定した引数名とは無関係に、引数の順番どおりに表示される。
univ(1)[1] "takushoku"univ(2)[1] "waseda"univ(3)[1] "kobe"univ(4)[1] "other"13. ループ (for関数)
・繰り返しとは「ループ」とも呼ばれ、ある条件が満たされる限り、同じ処理を繰り返すこと。
・for( ) 関数は「ベクトル X 内の全ての値 (i) について expression 1 から expression N まで実行せよ」というもので、次のように書く。
for(i in X){
expression 1
expression 2
expression 3
・・・
expression N
}・ベクトル X は数値や文字列を集めたもの。
・for ループは、ベクトル X 内の個々の値を一つずつ取り出し、初めから終わりまで、一度ずつ { } (Braces) 内のコードを実行する。
・i はシンボルなので何を使ってもいいが、通常、i が多く使われる。
・“I”, “Have”, “a”, “Dream”という四つの文字列の要素を順番に取り出してみる。
for(i in c("I", "Have", "a", "Dream")){
print(i)
}[1] "I"
[1] "Have"
[1] "a"
[1] "Dream"ループ内の i は関数に対する引数（ひきすう）のように機能する。
for ループは i という名前のオブジェクトを作り、ループを実行する度に ベクトル X 内
（この場合は“I”,“Have”,“a”,“Dream”) から個々の値を一つずつ取り出し i に新しい値を割り当てる。
まず最初の値である “I” を i に割り当てる。
次に二番目の値である “Have” を i に割り当てる。
全ての値が i に割り当てられるまでこのループ作業を繰り返す。
三つの値から構成されるベクトル (values) を 指定し、10 倍するループを作る
values <- c(1, 2, 3) # 1, 2, 3から構成されるベクトル (values) を指定
n <- length(values) # ベクトル (values) の値の数（この場合 3) を n と指定
results <- rep(NA, n) # 結果を置く空のコンテナ (NA) を作り、results と名前を付ける
for (i in 1:n) { # 初めから終わりまで (1:n) 、values の中から値を一つずつ取り出す
results[i] <- values[i] * 10 # 取り出した値に 10 を掛けた結果を results に置く
cat(values[i], "times 10 equal to", results[i], "\n") # 三つを順番に表示 (cat) させる
# "\n" で改行する
results # for ループは自動的に結果を返さないので、結果を表示するには cat() や print() を使う
}1 times 10 equal to 10
2 times 10 equal to 20
3 times 10 equal to 30 三つのベクトル a, b, c から構成されるデータフレームを作り、中央値を求めるループを作る
- 三つのベクトル a, b, c から構成されるデータフレームを作り、data と名前を付け、表示させる
data <- data.frame("a" = 1:5, # 1 から 5 までの連続変数を生成する
"b" = 4:8,
"c" = 7:11)
data a b c
1 1 4 7
2 2 5 8
3 3 6 9
4 4 7 10
5 5 8 11- data から一つずつ値を取り出して求めた「行の」中央値の結果を results に置いて表示させる
n <- length(data) # ベクトル (data) の数（この場合 3) を n と指定
results <- rep(NA, n) # 結果を置く空のコンテナ (NA) を作り、results と名前を付ける
for (i in 1:n) { # 初めから終わりまで (1:n) 、data の中から値を一つずつ取り出す
results[i] <- median(data[, i]) # 計算した三つの結果（行の中央値）を results に置く
cat("iteration", i, "'s median is", results[i], "\n") # 三つを順番に表示 (cat) させる
}iteration 1 's median is 3
iteration 2 's median is 6
iteration 3 's median is 9 ・cat( ) を使わず、結果だけ出力したければ
results[1] 3 6 91 から 100 までの整数の合計を求める方法
①sum( ) を使う方法
sum(1:100)[1] 5050②ループを使う方法
・手順は次のとおり
STEP1: 1から100までの整数を作る
x <- 1:100STEP2:i という空の（つまり０の）オブジェクトを用意する
i <- 0STEP3: 1 から 100 までの整数 x を最初から順番に一つずつ取り出し j に代入
for (j in x){
i <- i + j # 即ち i = i + j = 0 + 1 + 2 + ... + 100 = 5050
}STEP4: 結果を表示する
i[1] 5050STEP5: i = 1 に指定してみる
即ち i = i + j = 1 + 1 + 2 + … + 100 = 5051
x <- 1:100
i <- 1
for (j in x){
i <- i + j
}
i[1] 5051ループとベクトル操作の処理時間の違い
・1 から 100 までの整数の平方根を計算し、結果をベクトルに格納させてみる。
sqrt(1:100) [1] 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751
[8] 2.828427 3.000000 3.162278 3.316625 3.464102 3.605551 3.741657
[15] 3.872983 4.000000 4.123106 4.242641 4.358899 4.472136 4.582576
[22] 4.690416 4.795832 4.898979 5.000000 5.099020 5.196152 5.291503
[29] 5.385165 5.477226 5.567764 5.656854 5.744563 5.830952 5.916080
[36] 6.000000 6.082763 6.164414 6.244998 6.324555 6.403124 6.480741
[43] 6.557439 6.633250 6.708204 6.782330 6.855655 6.928203 7.000000
[50] 7.071068 7.141428 7.211103 7.280110 7.348469 7.416198 7.483315
[57] 7.549834 7.615773 7.681146 7.745967 7.810250 7.874008 7.937254
[64] 8.000000 8.062258 8.124038 8.185353 8.246211 8.306624 8.366600
[71] 8.426150 8.485281 8.544004 8.602325 8.660254 8.717798 8.774964
[78] 8.831761 8.888194 8.944272 9.000000 9.055385 9.110434 9.165151
[85] 9.219544 9.273618 9.327379 9.380832 9.433981 9.486833 9.539392
[92] 9.591663 9.643651 9.695360 9.746794 9.797959 9.848858 9.899495
[99] 9.949874 10.000000・この程度ならほぼ瞬時で計算可能。
・proc.time()関数を使って、プログラムの効率性（＝所要時間）を計測。
・1 から 1,000,000 までの整数の平方根を計算し、結果をベクトルに格納させてみる。
startTime = proc.time()
y = vector(mode = "numeric", length = 1000000)
for(i in 1:1000000){y[i] = sqrt(i)}
stopTime = proc.time()
elapsed = stopTime = startTime
show(elapsed) user system elapsed
3.341 0.297 3.824 ・研究室にある iMac(3.5GHz Intel Core i7)だと、約1.793秒かかる。
startTime = proc.time()
y = sqrt(1:1000000)
stopTime = proc.time()
elapsed = stopTime - startTime
show(elapsed) user system elapsed
0.007 0.003 0.009 ・ループを使わないベクトル操作だと、わずか0.053秒で計算が完了。
・ループを使った計算と比較すると、ベクトル操作だと約32倍の速度で計算できる。
14. ループ (while関数)
for( )関数では指定されたベクトルから要素を一つずつ取り出す回数が引数によって指定され、ベクトルの要素数によって決まる。
while( )関数では、繰り返しの回数が引数で指定されない。
i <- 0
while(i < 100){
i <- i + 1
}
print(i)[1] 100i <- 0
while(i < 100){
i <- i + 1
if (i %% 10 != 0) next
if (i %% 90 == 0){
cat("\n")
break
}
cat(i, "\t")
}10 20 30 40 50 60 70 80 15. ループの応用
サイコロを100回ふって、出た目の合計を計算する。
この作業を1000回実施し、その平均を計算する。
このシミュレーションの過程を一つずつ確認してみよう。
dice という名前を付けた人工のサイコロを作る。
dice <- 1:61000個の値を入れる場所を準備し、その場所に res という名前を付ける。
res <- numeric(1000) res の最初の6つのベクトルは次のとおり。場所を確保しただけで、まだ何もいれていないので数値は全て 0 である。
こういう場所が 1000 個できたことになる。
head(res)[1] 0 0 0 0 0 0サイコロを100回ふってその結果を hundred に収める。
hundred <- sample(dice, 100, replace = TRUE) サイコロを100回ふった結果は次の通り。
hundred [1] 2 2 3 3 5 4 5 1 1 6 5 1 4 1 1 5 4 2 3 5 1 1 3 3 4 2 3 2 2 3 2 3 6 3 5
[36] 6 2 3 4 2 4 4 4 4 4 6 2 5 3 5 4 4 3 6 2 4 5 6 4 3 2 1 4 3 6 1 4 1 1 6
[71] 6 4 1 1 4 3 6 2 5 1 2 3 2 4 2 4 3 5 3 4 2 3 2 2 6 5 3 3 1 5当然だが、正常なサイコロであれば、それぞれの目が1/6の確率で出るわけだから、この結果は妥当だといえる。
さらに、この結果をヒストグラムにすると次の様になる。
hist(hundred)
サイコロを100回ふる作業を1000回実施し、その平均を計算し図示するコマンドは次のとおり。
dice <- 1:6 #dice という名前を付けた人工のサイコロを作る
res <- numeric(1000) #1000個の値を入れる場所を準備、res という名前を付ける
hundred <- sample(dice, 100, replace = TRUE) #サイコロを100回ふってその結果を hundred に収める
hist(hundred) #結果をヒストグラムで表示させる16. Rを使ったプログラミング応用例
16.1 Monty Hall Problem
ここでは R 上でシミュレーションすることでモンティホール問題（ウィキペディア）を考えてみる。
【ゲームのルール】
- 3つのドア (1, 2, 3) に（高級車、ヤギ、ヤギ）がランダムに入っている。
- プレイヤーはドアを1つ選ぶ。（例えば 1 を選んだとする）
- モンティは残りのドアのうち 1 つを開ける。（例えば 2 を開けたとする）
- そのドアにはヤギが入っている。（2 にはヤギが入っている）
- プレイヤーはドアを選びなおすことができる。（つまり、3 を選び直すことができる）
Question: プレイヤーは最初に選んだドアに止まるべきか、それともドアを選び直すべきか？
意見 A : ドアを選び直しても高級車を得られる確率は同じだから、どちらでも良い。
意見 B : 高級車を得られる確率が上がるから、ドアを選び直すべき。この立場をとる説明の一例はこちら。
【R を使ったシミュレーションの前提】
- 高級車はドア 1 にある。
- プレイヤーは三つのドア（1，2，３）のうち無作為に一つのドアを選ぶ。
- （モンティがヤギの入っているドアを 1 つ開ける）
- プレイヤーは次の二つの選択ができる：
① 最初に選んだドアにとどまる
② 残ったドアに選択を変える
プレイヤーは 200 回無作為にドアを選び、①と②それぞれの場合に高級車が当たる確率を計算する。
シミュレーションの全体のスクリプトは ここからダウンロードできます。
# シミュレーションを行う為に 200個の入れものを準備。
set.seed(1838-3-11) # 乱数の設定。大隈重信の誕生日。
trials <- 200
# 「ドア 1」「ドア 2」「ドア 3」を無作為に 200 回選ぶ
prize <- sample(1:3, trials, replace = TRUE)
prize [1] 1 3 2 1 2 2 2 2 3 2 3 3 1 1 1 3 3 2 1 1 3 1 2 3 2 2 3 2 2 3 3 2 3 1 3
[36] 1 3 3 1 2 1 1 1 1 2 2 3 3 2 2 2 3 3 3 2 2 2 1 1 2 3 1 3 3 1 1 1 3 2 2
[71] 1 1 1 3 3 3 2 3 1 1 2 3 2 1 3 2 1 1 1 3 2 3 2 1 3 2 1 1 1 3 1 3 2 1 2
[106] 1 1 3 1 3 1 2 2 2 3 1 2 2 3 2 3 1 1 2 1 1 2 2 3 2 2 1 3 2 2 1 2 1 3 1
[141] 1 2 1 1 1 3 3 1 1 3 3 1 2 1 2 1 2 2 2 2 1 2 1 3 3 2 2 1 2 2 3 3 2 3 2
[176] 3 1 1 1 3 1 3 1 1 2 3 2 2 2 1 1 3 3 2 3 1 3 3 2 2「ドア 1」に高級車があると指定。
stay <- prize == 1
head(stay)[1] TRUE FALSE FALSE TRUE FALSE FALSEstay では 1 が TRUE、2 と 3 が FALSE と表示。
「ドア 2」と「ドア 3」に高級車はないと指定。
move <- prize != 1
head(move)[1] FALSE TRUE TRUE FALSE TRUE TRUEmove では 1 が FALSE、2 と 3 が TRUE と表示。
- Monty Hall の提案を受けず「ドア 1」に止まった時に高級車が当たる割合
prob.stayと「ドア 1」から動いた時に当たる割合prob.moveの累積を 200 回分計算する。
prob.stay <- prob.move <- rep(NA, trials)
for (i in 1:trials) {
prob.stay[i] <- sum(stay[1:i]) / i # 「ドア 1」に止まった時に高級車が当たる割合
prob.move[i] <- sum(move[1:i]) / i # 「ドア 1」から動いた時に高級車が当たる割合シミュレーション結果をプロットする。
plot(1:trials, prob.move, ylim = c(0, 1),
type = "l", col = "red",
xlab = "number of trials",
ylab = "prob. of getting a car",
main = "Simulation of Monty Hall Problem")
lines(1:trials, prob.stay, type = "l", col = "blue")
abline(h = c(0.33, 0.66), lty = "dotted") # 勝つ確率が 33% (1/3)と66% (2/3)に点線を引く
legend("topright", lty = 1, col = c("red", "blue"),
legend = c("change", "stay"))}
・実行方法
・R 上で「ファイル」=> 「新規文書」
=> 下記のコードをペースト
=> Script 上でコード全てを選択: Command + A (Mac), Ctr + A (Windows)
=> Command を押しながら Return (Enter) を押す (Mac の場合)
=> Cntl を押しながら R を押す(Windows の場合)
注意：RMarkdown 上で実行できないこともないが、動画アニメーションではなく「紙芝居」のような動きをする
# シミュレーションを行う為に 200個の入れものを準備。
set.seed(1838-3-11) # 乱数の設定。大隈重信の誕生日。
trials <- 200
prize <- sample(1:3, trials, replace = TRUE) # 「ドア 1」「ドア 2」「ドア 3」を無作為に 200 回選ぶ
prize
stay <- prize == 1 #「ドア 1」に高級車があると指定
head(stay)
move <- prize != 1 #「ドア 2」と「ドア 3」に高級車はないと指定
head(move)
prob.stay <- prob.move <- rep(NA, trials)
for (i in 1:trials) {
prob.stay[i] <- sum(stay[1:i]) / i # 「ドア 1」に止まった時に高級車が当たる割合
prob.move[i] <- sum(move[1:i]) / i # 「ドア 1」から動いた時に高級車が当たる割合
plot(1:trials, prob.move, ylim = c(0, 1), #シミュレーション結果をプロット
type = "l", col = "red",
xlab = "number of trials",
ylab = "prob. of getting a car",
main = "Simulation of Monty Hall Problem")
lines(1:trials, prob.stay, type = "l", col = "blue")
abline(h = c(0.33, 0.66), lty = "dotted") # 勝つ確率が 33% (1/3)と66% (2/3)に点線を引く
legend("topright", lty = 1, col = c("red", "blue"),
legend = c("change", "stay"))}Monty Hall Problem を JavaScript を使ってシミュレーションしているサイトもある。
16.2 Birthday Paradox
・スリラー作家アダム・ファウアーの小説『数学的にありえない』（上）の中で、大学の教室で授業を受けている58人の学生の中に「同じ誕生日のペアが少なくても1組いる確率」を問うシーンがでてくる(pp.175-188)。
・誕生日は1月1日から12月31日まで365日あるから、異なる2人の誕生日が同じ確率は1/365。
・2人が58人になったとしても、「同じ誕生日のペアが少なくても1組いる確率」がそれほど高まるとは、直感的には思えない。
・これをRを使って分析してみる。
・Rには、pbirthday()とqbirthday()という二つの関数が用意されている。
クラスの中で少なくても同じ誕生日のペアが1組いる確率 pbirthday()
・「クラスの人数」を横軸、「クラスの中で少なくても同じ誕生日のペアが1組いる確率」を縦軸に指定し、両者の関係を図示する。
dob <- 1:365
probs <- sapply(dob, pbirthday)
plot(dob, probs, type = "l", col = "tomato", lwd = 2, xlim = c(0, 100),
xlab = "The number of people in a group",
ylab = "Probability",
main = "Prob. that at least two people have the same DOB：pbirthday()")
abline(h = c(0, 1), col = "gray") ## 確率の0%と100%に横線を引く
abline(v = c(23, 58), col = "gray") ## クラス人数23人と58人に縦線を引く
・23人と58人のところに縦線を引いている。
・それぞれのクラス人数において「クラスの中で少なくても同じ誕生日のペアが1組いる確率」を求めてみると。
pbirthday(23, classes =365, coincident = 2) ## クラス人数が23人の時の確率[1] 0.5072972pbirthday(58, classes = 365, coincident = 2) ## クラス人数が58人の時の確率[1] 0.991665クラスの中で少なくても同じ誕生日のペアが1組いるために必要な人数 ｑbirthday()
・「クラスの中で少なくても同じ誕生日のペアが1組いるために必要な人数」を横軸、「クラスの人数」を横軸に指定し、両者の関係を図示する。
prob <- seq(0, 1, by = 0.01)
n.people <- sapply(prob, qbirthday)
plot(prob, n.people, typ = "l", col = "royalblue", lwd = 2,
xlab = "Prob. that at least two people have the same DOB",
ylab = "The number of people in a group",
main = "The number of people needed to have the same DOB：qbirthday()")
abline(v = c(0, 0.5, 0.99, 1), col = "gray") ## 確率の0%, 50%, 99%, 100%に縦線を引く
・50%と99%のところに縦線を引いている。
・50％、99%それぞれの確率を得るために必要なクラス人数を求める
qbirthday(0.5, classes =365, coincident = 2) ## 50％の確率を得るために必要な人数[1] 23qbirthday(0.99, classes = 365, coincident = 2) ## 99％の確率を得るために必要な人数[1] 57
References
- 浅野正彦, 矢内勇生.『Stataによる計量政治学』オーム社、2013年
- http://www2.kobe-u.ac.jp/~yyanai/classes/rm1/contents/
-
Kosuke Imai, Quantitative Social Science: An Introduction, Princeton
University Press, 2017. - John K. Kruschke, Bayesian Data Analysis, A Tutorial with R, JAGS and Stan, 2nd Edition, ELSEVIER INC., 2015.