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

Rと確率分布

1.確率と確率分布

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

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

表1 1つのサイコロの目と確率
X の値 1 2 3 4 5 6 合計
確率
1

  表1では目の値1、2、3、4、5、6 を X で表している。このように変数 X の値にはそれに対応する確率が存在する場合、その変数を確率変数と言う。表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] 5.267019  1.262095  5.904341  3.844168  5.441433  2.633377

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

> floor(temp)
[1] 5  1  5  3  5  2

  この値を、サイコロを6回投げたとき毎回でた目の値と考えればよい。このような方法で、関数 runif の引数 n を100、1000、10000にした結果のヒストグラムを図1に示す。n が100である一様分布のヒストグラムは次のコマンドで作成できる。関数 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にならないが、回数を1000、10000回に増やすと各々目のでる確率は1/6に近く。

図1 サイコロを投げるシミュレーションの結果

  このシミュレーションは母集団(調査対象全体)から標本(調査対象)を抽出する、いわゆるサンプリングのシミュレーションの基本である。

2.1 二項分布

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

 n 個から k 個を取りだす組み合わせの数 nCk は、Rでは関数 choose で計算できる。例えば、10個の中から3個を取り出す組み合わせの数は 10C3 は choose(10,3)で求ることができる。
 二項分布の確率計算は定義式を用いて計算してもよいが、Rには二項分布の確率を求める関数 dbinom がある。 n = 50、 p= 0.3 の二項分布 B (50,0.3) の確率変数 k = 20 の確率は次のように求めることができる。

> dbinom(20,50,0.3)
[1] 0.03703876

 確率変数 k が1、2、3、…、40をとる B (50,0.3)、 B (50,0.5) の2つの二項分布を図2に示す。分布の中心位置などは確率 p の値に依存することが図でわかる。

> x<-0:40
> plot(x,dbinom(x,50,prob=0.3),type="h",lwd=5,col="gray")
> plot(x,dbinom(x,50,prob=0.5),type="h",lwd=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 と言う品種の花弁の長さのデータの範囲は1〜1.9で、1.3〜1.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は次のコマンドで作成することができる。線の太さと色は好みに合わせて引数 lwd と col を用いればよい。

> 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 はそれぞれ分布の平均と分散である。平均が μ 、分散が σ2 の正規分布を N (μ,σ2) で表す。μ = 0 、σ2 = 1 の正規分布 N(0,1) を標準正規分布と呼ぶ。

2.3 χ2 分布

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

χ2 = X12 + X22 + … + Xn2

がしたがう分布を自由度 n の χ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 分布

2.4 t 分布

 確率変数 X が標準正規分布 N (0,1) にしたがい、確率変数 Y が自由度 n の χ2 分布にしたがうとき、

がしたがう分布を自由度 nt 分布と呼ぶ。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 t 分布

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

図8 t 分布と正規分布

2.5  F 分布

 確率変数 X が自由度 m の χ2 分布にしたがい、確率変数 Y が自由度 n の χ2 分布にしたがうとき、

がしたがう分布を自由度 (m,n) の F 分布と言う。F 分布の密度を求める関数は df である。次のコマンドで作成した3つ異なる自由度の分布 F を図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 F 分布

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)