[連載]フリーソフトによるデータ解析・マイニング 9

Rと確率分布


 

 

1.確率と確率分布

1つのサイコロを投げることを考えよう。通常では、各目がでる可能性は1/6であることを前提としている。この1/6が1つのサイコロを投げたとき、1つの目のでる確率である。しかし、サイコロを6回投げたとき各々の目がちょうど1回ずつでることは少ない。確率は、ある試行のくり返しの回数が十分多いとき、事象Aが起こった相対頻度である。つまり、試行回数が十分多いとき、事象Aが起こる相対頻度

 

 

の極限値を事象Aの確率と言う。通常事象Aの確率をで表す。サイコロを6回投げたときには各々の目がでる相対頻度は1/6にならないが、60回、600回、6000回、60000回と投げる回数を増やすと各々の目のでる相対頻度は1/6に近づく。1つのサイコロを投げたとき各々の目のでる確率を表1に示す。

 

1. 1つのサイコロの目と確率

の値

1

2

3

4

5

6

合計

確率

1

 

1では目の値123456で表している。このように変数の値にはそれに対応する確率が存在する場合、その変数を確率変数と言う。表1のような確率変数とその変数の取り得る確率を一緒にして確率分布と言う。表1のように確率変数が離散の場合は離散型確率変数、連続の場合は連続型確率変数と呼んでいる。

 

2.離散型確率分布

2.1 離散型一様分布

試行をくり返し行う際、取り得る結果が表1のように全て同じである離散型確率分布を離散型一様分布(uniform distribution )と呼ぶ。

 1つのサイコロを投げたとき各々の目がでる確率が1/6であることを確かめるためには、サイコロをくり返し投げることが考えられる。しかし、このような方法はあまりにも原始的で、時間と労力を無駄にする。

 プログラム言語やデータ解析ソフトにはこのような試行をシミュレーションするための乱数を発生させる機能がある。乱数(random numbers)とは,意図的にどの数値が多くでるようにと手を加えることなく、偶然にまかせた数値組である。Rでの離散型一様分布の乱数を発生させる関数はrunifである。関数runifは以下のように定義されている。

 

runif(n, min, max)

 

引数nは試行の回数、minは出力する乱数の最小値、maxは出力する乱数の最大値である。ただし、出力される乱数には引数の最大値は含まない。試行の回数を6、最小値を1、最大値を7にしたときに出力された乱数を以下に示す。乱数であるので、同じのコマンドをくり返し実行してもまったく同じの乱数が返される確率は非常に小さい。よって、読者が同じのコマンドを実行した場合、ここに示している結果と異なる結果が返されるのは自然である。

 

> temp<-runif(6,1,7)

> temp

[1] 4.577059 3.561613 6.876066 1.520917 4.693214 4.219530

 
このような乱数を用いてサイコロを投げるシミュレーションを行うことができる。その方法としては、次のように関数floorを用いて、出力値の小数点以下を切り捨てることが考えられる。

 

> floor(temp)

[1] 4 3 6 1 4 4

 
 この値を、サイコロを6回投げたとき毎回でた目の値と考えればよい。このような方法で、関数runifの引数n100100010000にした結果のヒストグラムを図1に示す。n100である一様分布のヒストグラムは次のコマンドで作成できる。関数histに用いた引数probability=Tは相対頻度であるか、絶対頻度であるかを指定する引数である。Tの場合は相対頻度、Fの場合あるいは省略した場合は絶対頻度のヒストグラムを作成する。
 
>temp<-floor(runif(100,1,7))
>hist(temp,breaks=c(0,1,2,3,4,5,6),probability=T,col="gray")
 
1で分かるように一様乱数を用いたサイコロを投げるシミュレーションでは、サイコロを100回投げた場合は、各々の目がでる確率は1/6にならないが、回数を100010000回に増やすと各々目のでる確率は1/6に近く。
1 サイコロを投げるシミュレーションの結果
 
 このシミュレーションは母集団(調査対象全体)から標本(調査対象)を抽出する、いわゆるサンプリングのシミュレーションの基本である。

 

2.1 二項分布

試行の取り得る結果が2通りしかない場合、その確率分布を二項分布(binomial distribution )と呼び、で表す。二項分布は次の式により定義される。

 

 

*個から個を取りだす組み合わせの数は、Rでは関数chooseで計算できる。例えば、10個の中から3個を取り出す組み合わせの数choose(10,3)で求ることができる。

二項分布の確率計算は定義式を用いて計算してもよいが、Rには二項分布の確率を求める関数dbinomがある。の二項分布の確率変数の確率は次のように求めることができる。

 

> dbinom(20,50,0.3)

[1] 0.03703876

 

 確率変数123、・・・、40をとる2つの二項分布を図2に示す。分布の中心位置などは確率の値に依存することが図でわかる。

 

> x<-0:40

> plot(x,dbinom(x,50,prob=0.3),type="h",lw=5,col="gray")

> plot(x,dbinom(x,50,prob=0.5),type="h",lw=5,col="gray")

 

2 二項分布の確率分布

 

 

2.連続型の確率分布

2.1 ヒストグラムと密度曲線

本節では連続の確率変数の例としてRの中にあるirisデータを用いることにする。irisのデータのsetosaと言う品種の花弁(Petal.Length)の長さのデータのヒストグラムを次のコマンドで作成したものを図3に示す。

 

>data(iris)

>hist(iris[1:50,3],probability=T,col="gray")

 

図で分かるように、setosaと言う品種の花弁の長さのデータの範囲は11.9で、1.31.5がもっとも多い。

  

3 setosaの花弁の長さのヒストグラム  図4 setosaの花弁の長さの密度曲線

 

相対度数のヒストグラムの総面積は1である。ヒストグラムの面積は横軸と次のような曲線で囲まれた面積の近似値である。この曲線をデータの密度(density)曲線と言う。ヒストグラムの密度は関数densityで求め、関数linesで曲線を描くことができる。

 

> dens<-density(iris[1:50,3])

>xlim<-range(dens$x)

>ylim<-range(dens$y)

>hist(iris[1:50,3],probability=T,col="gray",xlim=xlim,ylim=ylim)

> lines(dens,lwd=3,col="blue")

 

もしサンプル数が十分多くデータの範囲を区切る階級の間隔を十分小さくすると、そのヒストグラムの頂点をつないだ曲線は密度の曲線と重なる。密度曲線を分布曲線、あるいは分布密度とも呼ぶ。次に頻繁に用いられているいくつの連続型分布とRによる分布図の作成について紹介する。

 

2.2 正規分布

 正規分布とは図5のような左右対称で、平均と分散に依存する確率分布である。図5は次のコマンドで作成することができる。線の太さと色は好みに合わせて引数lwdcolを用いればよい。

 

> x<-seq(-4,4,0.1)

>curve(dnorm(x,0,0.5),from=-4,to=4,type="l")

>curve(dnorm(x,0,1),add=T)

>curve(dnorm(x,0,1.5),add=T)

>curve(dnorm(x,-2,1),add=T)

 

コマンドの中のseqは連続(sequence)に数値を生成する関数である。例えば、seq(-4,4,0.1)-4から4まで、0.1間隔を置いた数値セットを生成する。dnormは正規分布の密度を求める関数である。引数type="l"”l”lineの頭文字である。

 

5 正規分布

 

正規分布はもっとも多く使用されている確率分布で、次の式により定義されている。連続型確率変数の場合、このような関数を確率密度関数(略して密度関数)と呼ぶ。

 式の中のはそれぞれ分布の平均と分散である。平均が、分散がの正規分布をで表す。の正規分布を標準正規分布と呼ぶ。

 

2.3 分布

お互いに独立(ほかの確率変数の確率に影響されないこと)であり、同じの標準正規分布に従う確率変数があるとする。このとき、この確率変数の2乗の和

がしたがう分布を自由度(カイ2乗、chi-square)分布と呼ぶ。カイ2乗分布は自由度に依存する。カイ2乗分布の密度を求める関数はdchisqである。次のコマンドで作成した5つの異なる自由度のカイ2乗分布を図6に示す。

 

>x<-seq(0,20,0.1)

> curve(dchisq(x,2),from=0,to=20)

> curve(dchisq(x,4),add=T)

> curve(dchisq(x,6),add=T)

> curve(dchisq(x,8),add=T)

> curve(dchisq(x,10),add=T)

 

6 分布

 

2.4 分布

確率変数が標準正規分布にしたがい、確率変数が自由度分布にしたがうとき、

がしたがう分布を自由度t分布と呼ぶ。t分布の密度を求める関数はdtである。次のコマンドで作成した3つの異なる自由度のt分布を図7に示す。

 

> curve(dt(x,10),from=-4,to=4)

> curve(dt(x,3),add=T)

>curve(dt(x,1),add=T)

7 分布

 

8に標準正規分布と自由度1t 分布を示す。図から分かるようにt分布は正規分布より両すそが広がっている。自由度が大きくなるとt分布は正分布に近づく。

8 分布と正規分布

 

2.5 分布

 確率変数が自由度分布にしたがい、確率変数が自由度分布にしたがうとき、

 

がしたがう分布を自由度(,)分布と言う。分布の密度を求める関数はdfである。次のコマンドで作成した3つ異なる自由度の分布を図9に示す。

 

>x<-seq(0,3,0.1)

> curve(df(x,10,10),from=0,to=3)

> curve(df(x,5,10),add=T)

> curve(df(x,5,2),add=T)

 

9 分布

 

2.6 確率密度と乱数発生の関数

2にデータ解析でよく用いられている分布の確率密度と乱数発生の関数を示す。


 

2 Rにおける確率密度と乱数発生の関数

分布の総称

確率密度関数

乱数発生関数

一様(Uniform)分布

dunif(x, min=0, max=1,・・・)

runif(n, min=0, max=1)

二項(Binomial)分布

dbinom(x, size, prob,・・・)
rbinom(n, size, prob)

ポアソン(Poisson)分布

dpois(x, lambda,・・・)

rpois(n, lambda)

正規(Normal)分布

dnorm(x, mean=0, sd=1,・・・)

rnorm(n, mean=0, sd=1)

カイ2(Chi-square )分布

dchisq(x, df, ncp=0,・・・)

rchisq(n, df, ncp=0)

t分布

dt(x, df,・・・)

rt(n, df)

F分布

df(x, df1, df2,・・・)

rf(n, df1, df2)

ガンマ(Gamma)分布

dgamma(x, shape,・・・)

rgamma(n, shape)

ベータ(Beta)分布

dbeta(x, shape1, shape2,・・・)

rbeta(n, shape1, shape2)

対数正規(Lognormal )分布

dlnorm(x, meanlog = 0, sdlog = 1,・・・)

rlnorm(n, meanlog = 0, sdlog = 1)

ロジスティック(Logistic)分布

dlogis(x, ・・・)

rlogis(n)

指数(Exponential)分布

dexp(x, rate = 1, ・・・ )
rexp(n, rate = 1)

負二項(Negbinomail)分布

dnbinom(x, size, prob, mu,・・・ )

rnbinom(n, size, prob, mu)

多項(Multinomial)分布

dmultinom(x, prob, ・・・ )

rmultinom(n, size, prob)

幾何(Geometric)分布

dgeom(x, prob, ・・・ )
rgeom(n, prob)

超幾何(Hypergeometric)分布

dhyper(x, m, n, k, ・・・ )

rhyper(nn, m, n, k)

コーシー(Cauchy)分布

dcauchy(x,location=0,scale= 1,・・・ )
rcauchy(n,location=0,scale = 1)

ワイブル(Weibull)分布

dweibull(x,shape,scale=1,・・・  )
rweibull(n, shape, scale = 1)