確率を理解するための基礎となる順列（permutation）について勉強します。“イラスト・図解 確率・統計のしくみがわかる本 長谷川勝也（著）”の2.2節 順列の例題をRで解きます。
例題
以下のような日本の野球選手、アメリカの野球選手、ドミニカの野球選手が9人います。
日本: イチロー、松井、高橋
アメリカ: マグワイアー、グリフィー、ピアッツァ、グウィン
ドミニカ: ソーサ、ラミレス
（1）この9人で何通りの打順が作れますか？
（2）ソーサ、マグワイヤー、グリフィーは必ずクリーンナップにするとしたら、何通りの打順が作れますか？
（3）日本の選手をかたまって並べるものとすると、何通りの打順が作れますか？
重複なく対象物を並べたものを順列といいます。\(n\)個の対象物から\(r\)個を並べる順列の数は\({}_n \mathrm{P} _r\)通りあります。
\[ {}_n \mathrm{P} _r = \frac{n!}{(n-r)!} \]
9人の野球選手を1番から9番まで並べたものは順列です。当たり前ですが、野球選手に重複はないからです。9人を1番から9番まで並べた順列の数は\({}_9 \mathrm{P} _9\)通りあります。
\[ {}_9 \mathrm{P} _9 = \frac{9!}{0!} = 9! \]
Rでは階乗をfactorial関数で計算できます。
factorial(9)## [1] 362880よって、問（1）「9人で何通りの打順が作れますか？」の答えは362880通りです1。
問の解答は終わりましたが、これから打順の順列をRで作ります。
日本、アメリカ、ドミニカの9人の野球選手の名前のデータをRのオブジェクトに入力します。
ja <- c("イチロー", "松井", "高橋")
us <- c("マグワイアー", "グリフィー", "ピアッツァ", "グウィン")
do <- c("ソーサ", "ラミレス")
c(ja, us, do)## [1] "イチロー" "松井" "高橋" "マグワイアー"
## [5] "グリフィー" "ピアッツァ" "グウィン" "ソーサ"
## [9] "ラミレス"順列を作る前に、順列データのサイズ（bytes）を見積もります。\(データのサイズ = 順列の数\timesオブジェクトのサイズ\)
size <- factorial(9) * object.size(c(ja, us, do))
print(size, units="Mb") # Mb = Mbytes = 1024^2 bytes## 224.3 Mb私がRStudioをインストールしているホストのRAMは12 Gbytesですので、この例題の順列データのサイズより十分大きく、順列を作っても問題ありません2。
Rの標準パッケージに順列を作る関数はありませんので、整数\(n\)と整数\(r\)を引数にして、整数の順列の行列を返す関数myperm(n, r)を作ります3。
# Output the matrix of permutation nPn
permutations <- function (n) {
if (n == 1) {
return(matrix(1))
} else {
SP <- permutations(n-1)
p <- nrow(SP)
A <- matrix(nrow=n*p, ncol=n)
for(i in 1:n){
A[(i-1)*p+1:p,] <- cbind(i, SP+(SP>=i))
}
return(A)
}
}
# Output the matrix of permutation nPr
myperm <- function (n, r) {
if (n < r) {
stop("Permutation is intended for the case n >= r.")
}
P <- permutations(n)
if (n == r) {
return(P)
} else {
A <- unique(P[,1:r])
return(A)
}
}このmyperm関数は再帰的に順列を作ります。使い方は、例えば\(\{1,2,3\}\)から2つの数字を選び並べる順列を作るには次のようにします。
# {1,2,3}から2つを選び並べる順列
myperm(3, 2)## [,1] [,2]
## [1,] 1 2
## [2,] 1 3
## [3,] 2 1
## [4,] 2 3
## [5,] 3 1
## [6,] 3 2myperm関数と選手の名前を入力したオブジェクトc(ja, us, do)を組み合わせて、野球選手9人で1番から9番までの打順の順列を作るには次のようにします。
# Make the permutation of batting orders choice by 9 baseball players
P <- matrix(c(ja, us, do)[myperm(9, 9)], ncol=9)
# The number of rows of P
nrow(P)## [1] 362880# The number of columns of P
ncol(P)## [1] 9# The data size of P
print(object.size(P), units="Mb")## 24.9 Mb順列は、行数362880、列数9の行列の形をしています。もちろん、行数は\({}_9 \mathrm{P} _9 = 9!\)と一致しています。
ここまで計算して気づいたのですが、データのサイズが見積もりのサイズより実際のサイズが小さいのは、データサイズが大きい文字列（野球選手の名前）はRのオブジェクトのラベル名に使われていて、順列の中身は1から9の整数だからです。9の階乗くらいでデータサイズを気にする必要はありませんでした。
作成した順列は大きくてすべてを表示することはできませんので、sample関数でランダムに5つ表示します。
# 行番号が表示されるなど見やすくするためにデータフレームにします
df <- data.frame(P)
# 9番までの打順が1行に表示されるように設定します
options(width = 500)
# sample関数で1から362880までの整数をランダムに5つ選んで、選んだ整数の行番号のデータを表示します
df[sample(nrow(df), 5),]## X1 X2 X3 X4 X5 X6 X7 X8 X9
## 61566 松井 ピアッツァ 高橋 グウィン イチロー マグワイアー ラミレス ソーサ グリフィー
## 201341 グリフィー ラミレス ソーサ マグワイアー グウィン イチロー ピアッツァ 松井 高橋
## 114344 高橋 ソーサ ピアッツァ グウィン ラミレス 松井 イチロー グリフィー マグワイアー
## 245764 グウィン イチロー ソーサ マグワイアー 松井 高橋 ピアッツァ ラミレス グリフィー
## 202461 ピアッツァ イチロー 高橋 マグワイアー 松井 ラミレス グウィン グリフィー ソーサ順列が作れていることが確認できます。
ソーサ、マグワイアー、グリフィーの3人から3番、4番、5番を選び、残りの6人から1番、2番、6番から9番を選ぶ順列の数を求めます。
複数の試行（確率でいう実験）を同時に行う場合の起こり方は、それぞれの試行の起こり方の積です。試行1の起こり方が\(n_1\)通り、試行2の起こり方が\(n_2\)通り、・・・、試行rの起こり方が\(n_r\)通りあると、試行1から試行rを同時に行う場合の試行の起こり方は、\(n_1 \times n_2 \times ... \times n_r\)通りになります。
この問では、
という2つの試行が同時に起こると考えられます。
試行1の起こり方は、\({}_3 \mathrm{P} _3 = 3!\)通りあり、試行2の起こり方は\({}_6 \mathrm{P} _6 = 6!\)通りあるので、試行1と試行2が同時に起こるのは\(3! \times 6!\)通りです。
factorial(3) * factorial(6)## [1] 4320よって、問（2）「ソーサ、マグワイヤー、グリフィーは必ずクリーンナップにするとしたら、何通りの打順が作れますか？」の答えは4320通りです。
順列の数が求まったので、次は順列を求めます。
先に求めた9人を1番から9番まで並べた順列dfから3番、4番、5番がソーサdo[1]、マグワイヤーus[1]、グリフィーus[2]で構成されている順列を抽出するのが簡単そうです。
# dfの3列、4列、5列がソーサ、マグワイヤー、グリフィーの場合のサブセット
df2 <- subset(df,
(df[,3]==do[1] | df[,4]==do[1] | df[,5]==do[1]) &
(df[,3]==us[1] | df[,4]==us[1] | df[,5]==us[1]) &
(df[,3]==us[2] | df[,4]==us[2] | df[,5]==us[2])
)
nrow(df2)## [1] 4320df2[sample(nrow(df2), 5),]## X1 X2 X3 X4 X5 X6 X7 X8 X9
## 205618 ピアッツァ イチロー ソーサ グリフィー マグワイアー 高橋 グウィン ラミレス 松井
## 24191 イチロー ピアッツァ ソーサ グリフィー マグワイアー ラミレス グウィン 松井 高橋
## 63169 松井 ピアッツァ グリフィー ソーサ マグワイアー イチロー 高橋 グウィン ラミレス
## 41244 松井 イチロー マグワイアー グリフィー ソーサ ピアッツァ ラミレス グウィン 高橋
## 89441 高橋 松井 ソーサ マグワイアー グリフィー グウィン ラミレス イチロー ピアッツァ作成した順列の行数は\(3! \times 6!\)と等しいです。
次に、3番から5番の順列（試行1の結果）P_345と、3番から5番以外の順列（試行2の結果）P_126789をそれぞれ作り、組み合わせて2つの試行が同時に起こる順列df2bを作ります。
# ソーサ、マグワイアー、グリフィーの名前データ
smg <- c(do[1], us[1], us[2])
# ソーサ、マグワイアー、グリフィー以外の名前データ
others <- c(ja, us, do)[!(c(ja, us, do) %in% smg)]
# ソーサ、マグワイアー、グリフィー3人の順列
P_345 = matrix(smg[myperm(3, 3)], ncol = 3)
# ソーサ、マグワイアー、グリフィー以外6人の順列
P_126789 = matrix(others[myperm(6, 6)], ncol = 6)
# 2つの順列を組み合わせるためにデータをリピートしてデータフレームにする
tmp <- NULL
for (i in 1:nrow(P_345)) {
tmp <- rbind(tmp,
cbind(P_126789[,1:2],
P_345[i,1], P_345[i,2], P_345[i,3],
P_126789[,3:6])
)
}
df2b <- data.frame(tmp)順列の行列を作るときに注意する点は、扱うデータの塊が行か列か区別することです。cbind関数で列ごとに操作するか、rbind関数で行ごとに操作するか注意します。
9人の順列から抽出して作った順列df2と3人の順列と6人の順列を組み合わせて作った順列df2bが等しいかcompareパッケージを使い調べます4。
# Install compare
if (!require("compare")) {
install.packages("compare", repos="https://cloud.r-project.org/")
library(compare)
}## Loading required package: compare## Warning in library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, : there is no package called 'compare'## Installing package into 'C:/Users/maeda/Documents/R/win-library/3.4'
## (as 'lib' is unspecified)## package 'compare' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\maeda\AppData\Local\Temp\RtmpmmPT08\downloaded_packages##
## Attaching package: 'compare'## The following object is masked from 'package:base':
##
## isTRUEcompare(df2, df2b, allowAll=T)## TRUE
## sorted
## [X1] dropped [unused] levels
## [X2] dropped [unused] levels
## [X3] dropped [unused] levels
## [X4] dropped [unused] levels
## [X5] dropped [unused] levels
## [X6] dropped [unused] levels
## [X7] dropped [unused] levels
## [X8] dropped [unused] levels
## [X9] dropped [unused] levels
## renamed rows
## dropped row namescompare関数の結果がTRUEなのでdf2とdf2bは等しいです。
日本の選手だけかたまって並べるということは、「日本の選手」と残りの6人を1番から7番まで並べるという試行と、日本の選手3人から1番から3番までを並べるという試行を同時に行うと考えられます。
試行1の起こり方は\({}_7 \mathrm{P} _7 = 7!\)通りあり、試行2の起こり方は\({}_3 \mathrm{P} _3 = 3!\)通りあるので、試行1と試行2が同時に起こるのは\(7! \times 3!\)通りです。
factorial(7) * factorial(3)## [1] 30240よって、問（3）「日本の選手をかたまって並べるものとすると、何通りの打順が作れますか？」の答えは30240通りです。
順列の数が求まったので、次は順列を求めます。
問（2）と同様に、まずは、はじめに作った\(9!\)通りの順列データdfから日本の選手jaが連続して並んでいるデータを抽出して30240通りの順列df3を作ります。
tmp <- NULL
df3 <- NULL
for (x in 1:7) {
tmp <- subset(df,
(df[,x]==ja[1] | df[,x]==ja[2] | df[,x]==ja[3]) &
(df[,x+1]==ja[1] | df[,x+1]==ja[2] | df[,x+1]==ja[3]) &
(df[,x+2]==ja[1] | df[,x+2]==ja[2] | df[,x+2]==ja[3])
)
df3 <- rbind(df3, tmp)
}
nrow(df3)## [1] 30240df3[sample(nrow(df3), 5),]## X1 X2 X3 X4 X5 X6 X7 X8 X9
## 131777 マグワイアー 高橋 松井 イチロー グリフィー ソーサ ラミレス ピアッツァ グウィン
## 80991 高橋 イチロー 松井 ピアッツァ ラミレス グウィン グリフィー マグワイアー ソーサ
## 140839 マグワイアー グリフィー ラミレス ピアッツァ グウィン 松井 イチロー 高橋 ソーサ
## 318984 ソーサ ラミレス 高橋 イチロー 松井 グウィン ピアッツァ グリフィー マグワイアー
## 185425 グリフィー ピアッツァ ソーサ マグワイアー 松井 イチロー 高橋 グウィン ラミレス順列の数は30240通りで、サンプルは日本の選手3人がかたまって並べられています。
次に、「日本の選手」と残りの6人から作る順列（試行1の結果）P_1234567と、日本の選手の3人から作る順列（試行2の結果）P_x123を組み合わせて順列df3bを作ります。
P_1234567 <- matrix(c("ja", us, do)[myperm(7,7)], ncol=7)
P_x123 <- matrix(ja[myperm(3,3)], ncol=3)
tmp <- NULL
st <- NULL
for (i in 1:nrow(P_x123)) {
for (j in 1:nrow(P_1234567)) {
x <- grep("ja", P_1234567[j,]) # get the index of "ja"
if (x == 1) {
tmp <- c(P_x123[i,], P_1234567[j,2:7])
} else if (x == 7) {
tmp <- c(P_1234567[j,1:6], P_x123[i,])
} else {
tmp <- c(P_1234567[j,1:(x-1)], P_x123[i,], P_1234567[j,(x+1):7])
}
st <- rbind(st, tmp)
}
}
df3b <- data.frame(matrix(st, ncol=9)) # matrix() works to avoid the warning of duplicated row.names\(9!\)通りの順列dfから抽出して作ったdf3と、2つの順列を組み合わせて作った順列df3bをcompare関数で比べると両者が等しいことがわかります。
compare(df3, df3b, allowAll=T)## TRUE
## sorted
## renamed rows
## dropped row namesこのポストはR Markdownで書いてHTMLに自動変換して作りました（knitrがR MarkdownをMarkdownに変換し、pandocがMarkdownをHTMLに変換します）。R Markdownのインラインコードで` r factorial(9)` と書くと、変換したHTMLページではgetOption("digits")で設定している桁数（多分、デフォルトで7桁）で浮動小数点表示されます。つまり、R Markdownの` r factorial(9)`は、HTMLページで3.62880 10^{5}と表示されます。これだとひと目でわかりにくい。is.numeric(factorial(9)) = TRUEより、factorial関数の戻り値は整数型のオブジェクトなのにknitrはおかしな挙動をするものです。そこで、factorial(9)をインラインコードで整数表示するために、` r sprintf("%i", factorial(9))`としています。この方法は、Knitr with R MarkdownのRoundingセクションを読んで知りました。このポストの著者はsprintf関数で表示がうまくいかない場合があるから、オリジナルの関数myround()を使っているとのことです。また、knitrを設定して表示桁数をカスタマイズすることもできるそうです。Fixing knitr: Formatting statistical output to 2 digits in R - Jason A. FrenchやControl output - Yihui Xieにカスタマイズの方法が書いてあります。↩
私のアーカイブRで重複のない組み合わせをexpand.grid()とsubset()で作るに書いたように、expand.grid関数で\(9^9\)通りの組み合わせを作り、subset関数で\(9!\)通りの順列を抽出しようと考えましたが、データのサイズが200 Gbytes超と大きすぎるのでやめました。↩
permutations関数はGenerating all distinct permutations of a list in R - Stack Overflowを転用しました。↩
r - Sample random rows in dataframe - Stack Overflowを参考にして、compareパッケージを使いました。↩