Rで確率を学ぶ

これは、長畑(2009)『Rで学ぶ統計学』共立出版に基づく解説です。

学習項目です

• 事象と確率
• ベイズの定理
• 確率変数
• 確率分布

事象と確率

x <- floor(runif(100,1,7))

length(x[x==1])

sim <- function(n=100) {
  x <- floor(runif(n,1,7))
  for(y in 1:6) {
     cat(y,"=>\t",length(x[x==y]),"\t",length(x[x==y])/n, "\n")
       }
}

ベイズの定理

確率変数と確率分布

• 確率変数(random variable): 実数の値を取る変数Xの値が確率に基づくとき、Xを確率変数という。 そして、実際に取る値のことを実現値という
  例: サイコロを投げて出る目を表す変数は確率変数で、その実現値は１，２，３，４，５，６の六通りある。
• 離散(計数)型確率変数: サイコロの目のように飛び飛びの値をとる確率変数
• 連続(計量)型確率変数: 身長や体重のように連続した値を取る確率変数
• 確率分布: (離散型)確率変数の実現値をx₁, x₂,…,x_n、 それぞれの確率をp₁, p₂,…,p_nとすると、{p_i, i=1,…n}を確率分布という
• 分布関数 F(x): 確率変数Xの実現値がx以下である確率が以下のように書かれるときのF(x)
  
  特にXが離散型の場合 P(X=x_i)=p_i を確率関数(p.f., probability function)、 Xが連続型の場合(分布関数の定義に現れる) f(x)を(確率)密度関数(p.d.f., probability density function)という

[連続型確率変数の実現値と確率について]

『サイコロを振って出た目』を確率変数とみなすと、サイコロを振って出た目は1,2,3,4,5,6という飛び飛びの『実現値』を取る(例えば、3.5という値にはならない)ので、これは離散型確率変数と言われる。 それに対し、人の体重や身長を確率変数とみなすと、体重の場合0kgから300kgの間の値、身長なら0cmから300cmの間の値（連続した値）をとるので連続型確率変数と言って、離散型確率変数とは別の扱いになる。
ここで、離散型確率変数の場合、確率変数をX、その実現値をx_iで表すと、確率関数はp(X=x_i)と書ける。ところが、連続型確率変数の場合、確率変数をX、その実現値をx_iで表すとしても、確率関数はp(X=x_i)とは書けない。このことについて腑に落ちない人のためにここで説明する。

身長を確率変数とみなそう。するとクラスには実現値170cmの者がいるだろう。50人のクラスで10人がそうだとすれば「適当にクラスから人を選んでその身長が170cmである確率は10/50 = 0.20 」と言いたくなるかもしれない。しかしこれは次のような誤解と混同がある。

まず「誰かの身長がちょうど170cm」という値をとる確率はほとんど0である、という事実を指摘しておこう。なぜなら、170cmという値はどこかで『丸め』を行って得られた値だからである。計測器がcm刻みであったか、mm（ミリメートル）刻みであったか、μm（マイクロメートル）刻みか、nn(ナノメートル)刻みか、。。。ということを考えれば、cm刻みなら170cmということはありそうだが、nn刻みの場合、170cmちょうどいう値が得られるのはほとんど可能性がない、ということは容易に理解できよう。

つまり、我々が日常的に使っている「身長が170cmちょうど」というのは実は（例えば)169.5cmから170.5cmの範囲にある」ことの言い換えになっているのである。連続的に値を取る身長に対して確率のようなものを与える関数を確率密度関数(f(x)で表す)といい、連続型の確率変数の場合、ある値の幅の区間で定積分したもの（つまり、確率変数の実現値がある幅の値の範囲に入る可能性）が確率となるのである。ちなみに先の例でいえば、f(x)を区間[169.5, 170.5]の間で定積分した値が0.20 となる。

なお、身長をcm刻みで表示することにすれば、飛び飛びの値を取ることになり、この場合は「離散型確率変数」として扱われる。

分布関数と確率関数(密度関数)をグラフに書くと:

上段: 確率関数、下段:分布関数	上段: 密度関数, 下段: 分布関数

• 分布関数について次が成り立つ
  

  (1) F(x)は単調非減少

  (2)

  (3) F(x)は右連続な関数:

• 確率関数(密度関数)について次が成り立つ
  

  (1) pxi≥ 0 　　　　 (f(x) ≥ 0)

  (2)

演習問題3-1

裏表のあるコイン2枚を投げたときに、表の出る枚数の確率関数と分布関数を求めよ。ただし表か裏のどちらかになるものとする。 [ヒント]

裏表のあるコイン2枚を投げたとき、表が出る枚数を確率変数Xとする。 Xの実現値は0, 1, 2 の3通りであり、離散型の確率変数である。 確率関数 P(X=x)において、x≠0,1,2の場合は P(X=x=0 は明らか。 そこで、P(X=0), P(X=1), P(X=2)の値を求めれば良い。

例えばP(X=0)の値は、二枚とも裏の場合の確率なので、 1/2 * 1/2 = 1/4 = 0.25 と求められる。

分布関数は、P(X < 0)、P(X ≤ 0)、P(X < 1)、 P(X ≤ 1), P(X < 2)、P(X ≤ 2), P (X < ∞) を求めれば良い。明らかに P(X < 0) = 0, P(X < 2) = P(X ≤ 2), P(X < ∞) = 1.0である。

演習問題3-2

関数f(x )を、 0 < x　または π < x のときf(x ) = 0
0 ≤ x ≤ π のとき f(x) = (1/2)*sin(x) とする。
問1. f(x )が確率関数の性質(1)と(2)を満たすことを示せ。
問2. f(x )を確率関数としたときの分布関数を表示せよ。
[ヒント]

確率関数の性質(1)とは、0 ≤ f(x )
性質(2)は、区間[-∞ , +∞] におけるf(x )の(定)積分が 1 になることである。 どちらも示すのは容易であろう。
分布関数F(x )は、0 < x のときF(x)=0、 0 ≤ x ≤ π のときは、区間 [0, x] におけるf(x)の定積分と等しく、 π < x のときはF(x)=1である。これを表示すれば良い。

なお、この表示には、ベクトルに関数を適用した結果をベクトルにするsapply関数を使った： (分布関数をF(x)とすると、 x <- seq(-2,4,0.01); y <- sapply(x,F) ; plot(x,y,type="l")

確率変数Xが分布F(x)に従う確率変数であることを、X ～ F(x)と表す

• 期待値E(X): 確率変数 X の取る値とその確率の積の総和。(算術)平均、もしくは加重平均ともいう （注: 「平均」にはこれ以外に幾何平均や調和平均などがあるが、これらは期待値(=算術平均)とは別のもの)式で書くと以下：
  
  
• 分散V(X): 確率変数Xの取る値と期待 値との差の二乗の総和の期待値。式で表すと E((X - E(X))²)
  
  

期待値の性質:	a, bを定数とすると、 E(aX + b) = aE(X)+b
分散の性質:	V(X) = E(X2) - {E(X)}2
	V(aX + b) = a2 V(X)

• Xの原点周りのk次モーメント:
• Xの平均μ周りのk次モーメント:
• Xの積率母関数: φ(θ) = E[e^θX]

主要な(1変数の)確率分布

連続型分布

正規分布

正規分布(normal distribution): 密度関数f(x)が以下の式で与えられる分布を、 平均μ、分散σ²の正規分布（またはガウス分布)といい、N (μ,σ²)で表す:
特にμ=0, σ=1の正規分布 (N (0, 1)) を標準正規分布という

参考: 区間[-∞, ∞]におけるf(x)の積分の値は1、 区間[μ-σ, μ+σ]におけるf(x)の積分の値は約0.683、 区間[μ-2σ, μ+2σ]におけるf(x)の積分の値は約0.954、 区間[μ-3σ, μ+3σ]におけるf(x)の積分の値は約0.997

正規分布の性質：　X ～N (μ,σ²) のとき、E(X) = μ, V(X) = σ²

• 中心極限定理:X₁,…, X_nが互いに独立で、平均がいずれもμ、分散がσ²のとき(ここで、確率変数X_iは正規分布に従う必要はない)、データ数が増えてくるとその算術平均の分布は平均μ、分散σ²/nの正規分布で近似される
• X～N(μ,σ²⁾のとき、Y=aX+b(ここでa, bは定数)とすると、 演習問題4-1,および4-2から
• X₁～N(μ₁,σ₁²)、 X₂～N(μ₂,σ₂²)、かつX₁とX₂が独立のとき、次が成り立つ: a₁X₁+a₂X₂ ～ N(a₁μ₁+a₂μ₂, (a₁σ₁)²+ (a₂σ₂)²)

参考: 正規分布の密度関数の描画:

curve(dnorm,-4,4,main="正規分布の分布関数",xlab="",ylab="")
abline(h=0,v=0,lty=1)
segments(2,0,2,0.054, lty=3)
segments(-2,0,-2,0.054, lty=3)
segments(-1,0,-1,0.24, lty=3)
segments(1,0,1,0.24, lty=3)

[多変数の正規分布]

下に示される多変数の正規分布(ガウス分布)も重要である(xとμはd次元ベクトル、Σはd×dの行列(共分散行列), dは正整数)。 そして、次のような性質を持つ:

• n次元ベクトルX が N(μ,Σ) に従うとする。また、 A がm×nの行列、b がm次元ベクトルとすると、 Y = AX + b は N(Aμ + b, AΣA^TX) に従う (ここで A^TはA の転置行列)

その応用としては、例えば、ロボットの位置推定がある。 ロボットに例えば1m移動させようとしても、正確に1m移動することは不可能であり、1m移動した後の位置はその地点を平均とした正規分布に従うと考えられる。このような移動を繰り返していくと、位置はどんどん不確定になるが、そうであっても、ロボットの位置は正規分布で表すことができる。ここでセンサー(例えばGPS)によりロボットの位置についての何らかの情報が得られるとしよう。センサーで得られる情報はやはり誤差を含んだものである（車のナビを使っていて、道路を走っていても、ナビでは道路を外れたところを表示していることがしばしばある）。このセンサー情報も正規分布で表すことができる。カルマンフィルタはロボットの移動の情報とセンサーの情報により、確率に基づいて、ロボットの位置推定を行うための優れた手法の一つである。

(連続)一様分布

一様分布(uniform distribution): U(a,b)で表される。密度関数f(x)が次式で与えられる分布で、一様乱数が従う分布。a=1, b=2の場合の密度関数f(x)は右図:

ccurve(dunif(x,1,2), xlim=c(0,3), ylim=c(0,2), main="U(1,2)の密度関数")

指数分布

指数分布(exponential distribution): Exp(λ)で表す。密度関数f(x)が次式で与えられる分布で、 いくつかのλに対しf(x)は右図のようになる。待ち行列でよく利用される分布。

x <- seq(0,6,length=500)
u  <-  c(0.5, 1, 3,5)
y <- sapply(u, function(a) { a*exp(-a*x)})
matplot(x,y,type="l", lty=1, col=1:4, ylab="density", 
		main="指数分布のグラフ(λ=0.5, 1, 3, 5)")
legend(4,5,paste("λ=",u,sep=""),col=1:4, lty=1, bty="n")
abline(h=0, v=0, lty=2)

離散型分布

二項分布

二項分布(binominal distribution): B(n, p) で表される、ここでn, pの意味は例で示す: 「不良率p (0 < p < 1) の工程からランダムにn個の製品をとったとき、そのうちのx個が不良品である確率は次式で与えられ、 この確率変数Xは二項分布B(n,p)に従う(X ～ B(n,p))： P(X = x) = _nC_x p^x(1-p)^n-x 　　(x=0,1,…, n)

x=0:10
p<-seq(0.1, 0.9, 0.1)
y <- sapply(p,function(a){dbnom(x,10,a)})
y <- sapply(p,function(a){dbinom(x,10,a)})
matplot(x,y,type="o",lty=1,ylim=c(0,1),pch=1:9,col=1:9,ylab="density",
     main="二項分布")
legend(2,1.1,paste("p=",p,"(n=10)",sep=""),col=1:9,lty=1,bty="n")
abline(h=0,v=0,lty=2,col=2)

ポアソン分布

ポアソン分布(Poisson distribution): P_o(λ)で表される。ここでλは平均。 「1時間あたりに来店する客の数の平均をλ人であるとき、単位あたりの来客数Xがxである確率は次式で与えられ、Xは平均λのポアソン分布に従う(X ～ P_o(λ))：　P(X=x) = p_x = e^-λλ^x / x! 　　　(x=0,1,…)

x <- 0:15
lam <- c(0.5, 1,2,3,4,5,7)
y <- matrix(0,nrow=16, ncol=7)
for(j in 1:7) { d<-dpois(x,lam[j])
 	y[,j] <- d
}
matplot(x,y,type="o",lty=1,ylim=c(0,0.7),pch=1:7,col=1:7,ylab="density",
     main="ポアソン分布")
legend(6,0.7,paste("λ=",lam,sep=""),col=1:7,lty=1,bty="n")
abline(h=0,v=0,lty=2,col=1)

統計量の分布

カイ二乗(カイ自乗, カイ2乗、χ²)分布

χ²分布(Chi square distribution) χ²_n : X₁,…, X_nが互いに独立に標準正規分布N (0,1²)に従うとき、 は自由度nの カイ二乗分布に従うといい、T～χ²_n と表す。
自然数nに対し、Tの密度関数は次で与えられる:
なお、Γ(x)はガンマ関数で、その定義は:

x <- seq(0,10,length=500)
n <- c(1,2,3,6)
y <- sapply(n/2, function(a) { dgamma(x,shape=a,scale=2)})
matplot(x,y,type="l",ylim=c(0,1),lty=1,col=1:4,ylab="density",
     main="カイ二乗分布")
legend(5,0.8,paste("自由度=",n,sep=""),col=1:4,lty=1,bty="n")
abline(h=0,v=0,lty=3)

(スチューデントの)t分布

t分布(Student t distribution): Xが標準正規分布 N (0,1²)に従い、それと独立なYが自由度nのカイ二乗分布に従うとき、は自由度nの t分布に従うといい、T～t_n と表す。
自然数nに対し、Tの密度関数は次で与えられる:

x <- seq(-6,6,length=500)
n <- c(1,3,100)
y <- sapply(n,function(a) { dt(x,a)})
matplot(x,y,type="l",ylim=c(0,0.6),lty=1,col=1:3,ylab="density",
     main="t分布の密度関数")
legend(1,0.6,paste("自由度=",n,sep=""),col=1:3, lty=1,bty="n")
abline(h=0,v=0,lty=3)

F分布

F分布(F distribution): Xが自由度mのカイ二乗分布に従い、それと独立なYが自由度nのカイ二乗分布に従うとき、は自由度(m,n)のF分布に従うといい、T ～ F_m,nと表す。
自然数m, nに対してTの密度関数は次で与えられる:

x <- seq(0,6,length=500)
y <- matrix(0, 500,4)
n1 <- c(1,2,6,10)
n2 <- c(1,2,8,10)
for(j in 1:4) {
	d <- df(x,n1[j],n2[j])
	y[,j] <- d
}
matplot(x,y,type="l",ylim=c(0,1.2),lty=1,col=1:4, ylab="density",
     main="F分布の密度関数")
legend(3,1,paste("n1=",n1," n2=",n2,sep=""),col=1:4,lty=1,bty="n")
abline(h=0,v=0,lty=3)