一様分布
離散型一様分布
質的データ (例:勝ち/負け、左右、赤/青/白、東西南北) で、各値が発生する確率が等しいときの確率分布を離散型一様分布といます。
例えばサイコロの目は、最小値が1で最大値が6の離散型一様分布だと考えられます。
サイコロの一様分布は、各目が出る確率が1/6で等しいと考えられるからです。
また、ジャンケンのグー、チョキ、パーが出る確率も全て1/3で等しい離散型一様分布だと考えられます。
離散型一様分布の確率関数
離散型一様分布は以下の確率関数で表されます。
\( f(x) = \frac{1}{n} \)
\( f(x) = 0 (ただしxが上記N以外の選択肢) \)
ここでnは値の数であり、サイコロならば6、ジャンケンならば3です。
確率分布の総面積は1であることも確認しておきましょう。1/6が6個あるから、全部で1です。また、1/3が3個あるから、全部で1です。
離散型一様分布からランダムなデータを得る
sample()を利用して、離散型一様分布からランダムな値を得ることができます。
例えば、サイコロのように6つの選択肢から、ランダムに3回値を出すには
sample(x=6, size=3)
です。ただし、上記の呪文では同じ目は出なかったでしょう。
この呪文は、1つの数字を一度しか選び出さないデフォルト設定(replace=FALSE)になっているためです。同じ数字を繰り返し選んでもかまわないようにする場合、(replace=TRUE) に設定変更する必要があります。
sample(x=6, size=30, replace=T)
離散型一様分布のヒストグラム
このsampleの数字から、ヒストグラムを描いてみましょう。ちなみに、sampleを実行するたびにランダムに数字が選ばれるので、以下の呪文によって、上と同じ結果にはなるわけではありません。
hist(sample(x=6, size=30, replace=T), breaks=seq(0,6,1))
よっぽど幸運な人でないかぎり、30個のデータでサイコロのような試行をしても、各値がだいたい5個ずつ選ばれる、ということにはなりません。では、もっと沢山の数字を選ばせるとどうなるでしょうか。
hist(sample(x=6, size=30000, replace=T), breaks=seq(0,6,1))
#30000個のデータを取れば、各値がほぼ5000ずつ選ばれる結果になったことと思います。よっぽど悪運な人以外は。とはいえ、6個の目が出た回数が全く同じ数字になることもめったにありません。
データ数を変えてみれば、データ数が少ないと、なかなか等頻度にならないことや、データ数が多いと等頻度に近づくことも分かります。ここではデータ数が20の場合と300000の場合を試してみます。
hist(sample(6, size=300, replace=T), breaks=0:6)
この呪文を何度も繰り返してみると、サイコロのような事象を扱う場合、データ数が300くらいでは「ちょっと運が悪い」だけで、各値が6分の1で選ばれているとは感じられないデータになることがあることも実感できますね。
hist(sample(1:6, size=300000, replace=T), breaks=0:6)
データ数が300000になると、ほぼ等しくなりますが、それでもわずかな凸凹は、なかなかなくなりません。
連続型一様分布
連続変数のデータ (例:長さ、重さ、時間) で、全ての値が発生する確率が等しいときの確率分布を連続型一様分布といます。離散型一様分布と違って、連続型一様分布では、同じ数字が2度起こることはまずありません。これは連続型一様分布だけでなく、正規分布など、連続変数を想定している他の確率分布も同じです。
連続型一様分布の確率密度関数
連続型一様分布は以下の確率関数で表されます。
\( f(x) = \frac{1}{max - min} \)
\( f(x) = 0 (ただしx < minあるいはx > max) \)
ここでmaxは確率分布の最大値、minは最小値です。
離散型との違いが明確になるように、具体例を挙げておきます。
サイコロのように最小値が1、最大値が6だけど、連続型の場合は6-1=5なので、2.12や3.56などの各値となる確率密度は1/5です。X軸の長さが5で、Y軸である高さが1/5だから、やはり確率分布の総面積は1ですね。
確率と確率密度
連続変数を想定している確率分布で、ある値になる「確率」は実質的に0なのです。その理由は、ある値の両隣の値との距離が無限小だからです。不連続変数(整数)や質的データ(赤、白、黄色などをそれぞれ1色とみなせる)では、各値と両隣の値の間には、「1」という距離がありますが、連続変数では、例えば3とその隣の数字(3.000......1, ただし0は無限大に続く)の間には無限小の距離しかありません。けれども最小値よりも小さな値が得られる確率は、本当に0であるのに対して、最小値と最大値の間の数字が得られる確率は、とても小さいだけで、真に0というわけではありません。このような事情から、連続変数の確率分布では、確率ではなく、確率「密度」を求めなければなりません。確率密度とは、本当は無限小の確率に対して、X軸が「1」だけ幅があるように値を変換したものです。単位X軸の値あたりなので「密度」というわけです。
例として、連続型一様分布の各値の確率密度を計算してみましょう。
呪文はdunifです。dは密度densityという意味です。例えば最小値が1で最大値が6の連続型一様分布で、値が4になる確率密度を求めると
dunif(x=4, min=1, max=6)
たしかに0.2と出力されましたね。1以上6以下という範囲以外の数値なら確率密度が0となることも確認しておきましょう。
dunif(x=7, min=1, max=6)
連続型一様分布からランダムなデータを得る
連続型一様分布から、ランダムに数値を選ぶ場合はrunif()を用います。
runifのrはrandomという意味。unifはuniform、つまり「一様」という意味です。
今回は最小値1、最大値6の連続型一様分布から3つの値をランダムに選んでみましょう。
runif(n=3, min=1, max=6)
この呪文を何回繰り返しても、理論上は同じ値は出ないはずです。実際にはRが小数点以下6桁までしか表示しないので、100000分の1で同じ値が出ますが。
連続型一様分布の作図
連続型一様分布を作図してみましょう。
まず、0から7までの範囲で0.01間隔の数列xを設定します。
そしてx の各値における確率密度をもとに、折れ線グラフにして図に記します。
x <- seq(0, 7, 0.01)
plot(x, dunif(x, min=1, max=6))
#plot()を折れ線グラフにするにはtype="l"と設定します。ここで、l はエルである。数字のイチではないので、注意してください。
plot(x, dunif(x, min=1, max=6), type="l")
#本当は、一様分布だから、あいだの数字は必要ありません。つまり、0, 0.99, 1, 6, 6.01, 7くらいの数字だけで、同じグラフになります。
x <- c(0, 0.99, 1, 6, 6.01, 7)
plot(x, dunif(x, min=1, max=6), type="l")
連続型一様分布の最尤推定
連続型一様分布の母集団推定を最尤推定によって試みてみましょう。ここから先は難しい話なので、興味がある人だけどうぞ。
さて、連続型一様分布の母数は最小値と最大値なので、これらを推定します。
その前に、肩慣らしの具体的類題として、最小値0、最大値2の連続型一様分布、つまり0≤x≤2のとき、f(x) = 1/(2-0) =0.5となる分布で、複数の値に関する確率密度を求めてみます。仮に、データは2と1.5、0の3つのデータだとします。まず、各データの確率密度を求めます。
dunif(2, min=0, max=2)
dunif(1.5, min=0, max=2)
dunif(0, min=0, max=2)
#尤度は、各データが得られる確率の総乗なので、これらを全て掛けた値となります。
dunif(2, min=0, max=2) * dunif(1.5, min=0, max=2)*dunif(0, min=0, max=2)
#総乗の関数はprod()です。この関数を使うと、尤度は以下のように求めることができます。
prod(dunif(c(2, 1.5, 0), min=0, max=2))
#最大値が2より大きい場合や、最小値が0より小さい場合は、尤度が0.125よりも小さくなることを確認しておきましょう。
prod(dunif(c(2, 1.5, 0), min=0, max=2.1))
prod(dunif(c(2, 1.5, 0), min= - 0.1, max=2))
#逆に、最大値が2より小さい場合や、最小値が0より大きい場合、データが2や0となる確率は0なので、尤度も0になります。
prod(dunif(c(2, 1.5, 0), min=0.1, max=2))
prod(dunif(c(2, 1.5, 0), min= 0, max=1.9))
では、もう少し難しくやってみましょう。いわゆるシミュレーションです。
まず連続一様分布のデータをランダムにとります。そして、まずは、その平均、最小値、最大値を計算します。
n <- 10
x <- sample(seq(0, 5, 0.01), n)
x
mean(x) #平均である必要はないのですが、minとmaxの間の値ということです。
min(x)
max(x)
#ここから先、求めた尤度を一時的に保管しておく場所が必要なので、それをyudo, yudo2, yudo3などとしておきます。
#これらの中身は、最初は空っぽなので、ただc()とだけ書きます。
yudo <- c()
yudo2 <- c()
yudo3 <- c()
#以下で for (i in 1:1000){ } となっているのは「iに1から1000までの整数を代入して{ } 内の作業を繰り返せ」という意味です。
その{ }内にfor (j in 1:1000){ }が入っているから、iに1を代入した状態で、jに1から1000までの数字を代入して作業を繰り返して、それが終わったらiに2を代入した状態で、同じように繰り返して、、、という作業をiに1000を入れるまで繰り返す、ということになります。
手作業でやっていたら文字通り日が暮れてしまいますが、Rはこの作業をかなり速くやってくれます(10秒くらい?パソコンの能力次第だけど)。
for(i in 1:1000){<br> for(j in 1:1000){<br> prod(dunif(x, min=mean(x)-0.01*i, max=mean(x)+0.01*j)) -> yudo[j]<br>}
max(yudo) -> yudo2[i]
}
saisyo <- mean(x)-0.01*1:1000
plot(saisyo, yudo2)
saisyo[yudo2==max(yudo2)]
for(i in 1:1000){<br> for(j in 1:1000){<br> prod(dunif(x, min=mean(x)-0.01*j, max=mean(x)+0.01*i)) -> yudo[j]<br>}
max(yudo) -> yudo3[i]
}
saidai <- mean(x)+0.01*seq(1,1000,1)
plot(saidai, yudo3)
saidai[yudo3==max(yudo3)]
このように、難しーく最尤推定しましたが、要するに最尤推定すれば「データの最大値が連続型一様分布の最大値の推定値、データの最小値が連続型一様分布の最小値の推定値」となるので、計算するまでもないことです。ただ、このようなシミュレーションに挑戦しておくと、Rでのシミュレーションに慣れることにもなりますし、最尤推定がどのようなことをするのかがよく分かるのではないかと思って紹介しました。