明治大学新領域創造

統計特論３

東京大学大学院農学生命科学研究科　大森宏

この講義の目的

　統計解析ソフトＲを用いて，統計解析が自在に行えるスキルを身に着ける．また，解析的証明より，シミュレーションによる数値的な証明（厳密ではないが直感的に理解しやすい）を行い，統計学の視覚的で直感的な理解をめざす．

Rの基本的使い方を学ぶ．

**今日のスクリプト**
x <- c(2, 3, 5)		#ベクトルの定義
y <- c(10, 14, -7)
a <- matrix(1:9, nrow=3)		#行列の定義
b <- x %*% t(x)		#行列の積
z <- c(6, -5, 2)
a1 <- cbind(x,y,z)		#列ベクトル－＞行列生成
solve(a1)		#逆行列
eigen(a1)		#固有値
x <- 1:10		#連続した自然数
y <- c(3,7,8,6,9,10,6,5,10,7)
plot(x, y)		#散布図
plot(x, y, type="l")		#直線でつなぐ
x <-seq(-20, 20, by=0.1)		#等間隔点列
y <- x^2		#関数定義
y <- 3x^3 - 2x^2 + 6*x + 5
y <- -2cos(x)+4sin(x)
plot(x, y, type="l")		#関数のグラフ
plot(x, y, type="l", xlim=c(-20,20), ylim=c(-20,20))		#表示部分の指定

０．方程式の根

ニュートン法

　Rでは，ニュートン法により方程式の根を求めることができる．

R による方程式の解法のスクリプト
# 関数の作成
f1 <- function(x){
exp(x)-2
}
#関数のグラフ
x <- seq(-1, 2, by=0.01)
plot(x, f1(x), type="l")
# これでもよい
curve(f1(x), -1,2)
#水平線と垂直線
abline(h=0)
abline(v=0)
#グラフタイトル
title(main="f(x) = exp(x) - 2 のグラフ")
# 0 < x < 1 で根を探す
uniroot(f1, c(0,1))
exp

R による方程式の解法のスクリプト
# 関数の作成 f1 <- function(x){ exp(x)-2 } #関数のグラフ x <- seq(-1, 2, by=0.01) plot(x, f1(x), type="l") # これでもよい curve(f1(x), -1,2) #水平線と垂直線 abline(h=0) abline(v=0) #グラフタイトル title(main="f(x) = exp(x) - 2 のグラフ") # 0 < x < 1 で根を探す uniroot(f1, c(0,1))

　方程式，exp(x) - 2 = 0，の根は，x = log(2) = 0.6931472，で，ニュートン法による近似解は， x = 0.6931457，であり，四捨五入した値で比較すると，小数第４位まで正確であった．

ニュートン法の R の出力
> uniroot(f1, c(0,1))          #根を探すコマンド　
$root
[1] 0.6931457          方程式の根　
$f.root
[1] -2.943424e-06        　
$iter        　
[1] 5      ニュートン法の反復回数　
$estim.prec      根の有効数字　
[1] 6.103516e-05      正確さの桁数は小数点4位まで　
> log(2)          #Rのコマンド　
[1] 0.6931472      正確な根　

多項式方程式の根

R による方程式の解法のスクリプト
# ２次方程式
# (x + 1)(x + 2) = 2 + 3x + x^2 = 0 の根
polyroot(c(2, 3, 1))
# ３次関数のグラフ
# f(x) = x³ - x² - 2x + 2，のグラフ（赤色）
curve(2 - 2*x - 1*x^2 + 1*x^3, -2, 2.5, col="red")
# f(x) = x³ - x² - 2x + 3，のグラフ（黒色）の重ね書き
curve(3 - 2*x - x^2 + x^3, -2, 2.5, add=TRUE)
#水平線と垂直線
abline(h=0)
abline(v=0)
title(main="３次関数のグラフ")
# ３次方程式
# (x^2 - 2)(x - 1) = 2 - 2*x - 1*x^2 + 1*x^3 = 0
# ３つの実数根
polyroot(c(2, -2, -1, 1))
# 3 - 2*x - 1*x^2 + 1*x^3 = 0
# 虚数解もでる．
polyroot(c(3, -2, -1, 1))

3th

多項式方程式の R の出力
> polyroot(c(2, 3, 1))          # 2 + 3x + x^2 = 0 の根　
[1] -1+0i    -2-0i          "0i" なので，実数解，x = -1，-2 である．　
> polyroot(c(2, -2, -1, 1))          # 2 - 2*x - 1*x^2 + 1*x^3 = 0，の根　
[1] 1.000000-0i    -1.414214+0i    1.414214+0i　
         "0i" なので，実数解，x = 1, -1.414214（≒-√2），1.414214，である．
> polyroot(c(3, -2, -1, 1))          # 3 - 2*x - 1*x^2 + 1*x^3 = 0，の根　
[1] 1.273409+0.563821i    -1.546818-0.000000i    1.273409-0.563821i　
         実数解：x = -1.546818，虚数解：x = 1.273409 ± 0.563821i　

１．確率分布

離散確率分布

　ある量の集まり P = { p₁，…，p_n } の中で，

(1. 1)

という性質をもつものを離散（discrete）確率分布（probability distribution）という．
　各 p_i を確率密度（probability density）という．なお，n は可算無限（自然数と対応づけられる）であるならば離散的である．

離散確率変数

　離散的な変数 X = { x₁，…，x_n } のおのおのの値に対し，それが生起する確率 p_i が与えられているとき，X を離散確率変数（discrete random variable）という．これは，

表 1 ：離散確率変数（n = 5）
　変数 X 　x₁　　x₂　　x₃　　x₄　　x₅　

　確率 P 　p₁ 　p₂ 　p₃ 　p₄ 　p₅

表 1 ：離散確率変数（n = 5）
変数 X	x₁	x₂	x₃	x₄	x₅
確率 P	p₁	p₂	p₃	p₄	p₅

と表せる．

連続型確率分布

　関数 f(x) が，

(1. 2)

という性質を持つとき，これを連続型確率分布（continuous distribution）という．また，f(x) を確率密度関数（probability density function）とも呼ぶ．なお，確率密度関数 f(x) が，パラメータ θ を持つとき，f(x；θ) と表記することもある．

連続型確率変数

　連続な変数 X が分布 f(x) をもつとき，連続型確率変数（continuous random variable）もしくは変量（variate）という．本稿では，X ～ f(x)，と表記する．

累積分布関数

　Pr[A] を事象 A が生起する確率とする．連続型確率変数 X に対し，X が x 以下である確率，

(1. 3)

で定義される関数 F(x) を累積分布関数（cumulative distribution function）という．この関数を用いると，確率変数 X が区間 (a, b) に落ちる確率は，

(1. 4)

で表せる．なお，離散型確率変数でも積分を和に変えることにより，同様に階段型の累積分布関数が定義できる．

２．分布の代表値

平均

　平均（mean）μ は，分布の中心的な位置（location）座標を表す．確率変数 X に対しては， X の期待値（expectation） E[X ] とも表記し，

(2. 1)

と定義される．離散確率変数では，積分を総和に変えることにより括弧内のように定義できる．

分散

　分散（variance） σ² は，分布の拡がり（dispersion）の程度を表す．確率変数 X に対しては，Var[X ] と表記する．確率変数 X の関数 (X - μ)² の期待値でもあり，

(2. 2)

と定義される．
　また，σ を標準偏差（SD : Standard Deviation）といい，平均と同じ次元で分布の拡がりの大きさを表す量である．

積率母関数

　確率変数 X の関数 exp(tX ) の期待値

(2. 3)

を積率母関数（moment generating function）という．積率母関数を t で r 回微分すると，

となるので，t → 0，とすれば，

(2. 4)

となり，確率変数 X の r 次の積率（moment）E[X ^r] が比較的簡単に求まる．

３．統計的独立

独立

　２つの確率変数 X ～ f(x)，Y ～ g(y) に対し，その同時分布（joint distribution）の密度関数を h(x, y) とする． X と Y が互いに独立（independent）であるのは，同時分布が h(x, y) = f(x)g(y) と変数分離される場合である．

独立な確率変数の平均と分散

　２つの独立な確率変数 X，Y の平均と分散がそれぞれ，

E[X ] = μ_x，Var[X ] = σ_x²， E[Y ] = μ_y，Var[Y ] = σ_y²，

であるとき，

E[X + Y ] = E[X ] + E[Y ] = μ_x + μ_y， Var[X + Y ] = Var[X ] + Var[Y ] = σ_x² + σ_y²
E[X - Y ] = E[X ] - E[Y ] = μ_x - μ_y， Var[X - Y ] = Var[X ] + Var[Y ] = σ_x² + σ_y²

である．一般に，スカラー a，b に対して

E[aX + bY ] = aE[X ] + bE[Y ] = aμ_x + bμ_y， Var[aX + bY ] = a²Var[X ] + b²Var[Y ] = a²σ_x² + b²σ_y² (3. 1)

である．

無作為標本

　分布 f(x) から大きさ n の標本 X₁，…，X_n を抽出したとき，それらが互いに独立であればその同時分布は

(3. 2)

と因数分解される．このような標本を無作為標本（random sample）という．また，標本全体を母集団（population）といい，f(x) を母集団分布ともいう．母集団分布の平均や分散などを母集団母数（parameter）という．

大数の法則

　平均が μ である分布をもつ母集団から大きさ n の無作為標本 X₁，…，X_n を抽出したときに，

(3. 3)

を標本平均（sample mean）という．このとき，標本の大きさ（サンプルサイズ(sample size)）n を大きくしていくと，標本平均 X^- は母集団平均 μ に近づく（確率収束する）．すなわち，どんな小さな正の数 ε に対して，

(3. 4)

が成り立つ．これを大数の法則（law of large numbers）という．
　大数の法則は，母集団からたくさんの標本を取れば取るほど，より正確に母集団分布についての推論が行えることを保証している．また，母集団分布に従う乱数が正しく生成できるならば，コンピュータシミュレーションにより母集団母数ははぼ正確に求めることができる．さらに，母集団分布に従う乱数で生成される経験分布関数は，母集団分布関数に近づいていく．

４．１変量分布

4-1．離散分布の例

離散一様分布（discrete uniform distribution)

　N 個のセルがあり，各セルの生起確率が互いに等しく，Pr[X = i] = 1/N である分布．
　平均： E[X ] = (N + 1)/2，分散： Var[X ] = (N² - 1)/12．
たとえば，サイコロの出る目の分布．

表 2 ：離散一様分布
　変数 X 　1　　2　　3　　4　　5　　6　

　確率 P 　1/6 　1/6 　1/6 　1/6 　1/6 　1/6

表 2 ：離散一様分布
変数 X	1	2	3	4	5	6
確率 P	1/6	1/6	1/6	1/6	1/6	1/6

ベルヌイ分布（Bernoulli distribution)

　成功（X = 1）確率が p，失敗（X = 0）確率が q = 1 - p である分布．
　　平均： E[X ] = 0・q + 1・p = p，分散： Var[X ] = E[X²] - (E[X ])² = 0²・q + 1²・p - p² = p(1 - p) = pq

二項分布（binomial distribution)

　独立な n 回のベルヌイ試行を行ったときの成功回数 X の分布．その確率密度は，

(4. 1)

である．X ～ B(n，p) と書くこともある．
m(t) = E[e^tx] = Σe^tx_nC_x p^xq^n-x = Σ_nC_x (pe^t)^xq^n-x = (pe^t + q)ⁿ
より，m'(t) = npe^t(pe^t + q)^n-1， m''(t) = n(n - 1)(pe^t)²(pe^t + q)^n-2 + npe^t(pe^t + q)^n-1
平均：　E[X ] = m'(0) = np，分散： Var[X ] = m''(0) - (np)² = n(n - 1)p² + np - (np)² = np(1 - p) = npq

３割バッターが１試合（５打席）で打つヒット数の分布
３割バッターが５打席で打つヒット数は，成功確率 p = 0.3 のベルヌイ試行を n = 5 回行ったときに成功する回数の分布で，n = 5，p = 0.3 の二項分布 B(5, 0.3) に従う．そのグラフは，以下のようになる．

**二項分布のグラフの R スクリプト**
n <- 5	#試行回数
x <- 0:n	#回数
p <- 0.3	#成功確率
hit <- dbinom(x, size=n, prob=0.3)	#二項確率
plot(x, hit, type="h", ylim=c(0,0.4), xlim=c(0,n), cex.lab=0.8, xlab="ヒット数", ylab="確率密度")
title(main="ヒット数の分布（n=5，p=0.3）")	#タイトル
sum(x*hit)	#平均（np = 1.5）
sum(hit(x - np)^2)	#分散（np(1 - p) = 1.05）

Weldon のサイコロ実験
イギリスの統計学者 Weldon は，12個のサイコロを同時に投げ，5か6の目が出た個数をカウントする実験を26306回行った．その結果以下のデータを得た．

   5,6の個数       0    1    2    3    4    5    6

出た回数    185       1149       3265       5475       6114       5194       3067

   5,6の個数       7    8    9    10    11    12    合計

出た回数    1331       403       105       18       0       0    26306

5,6の個数	0	1	2	3	4	5	6
出た回数	185	1149	3265	5475	6114	5194	3067
5,6の個数	7	8	9	10	11	12	合計
出た回数	1331	403	105	18	0	0	26306

まず，データから統計量を求める．
平均：x^- ＝ (0*185＋1*1149＋2*3265＋…＋12*0)/26306＝4.052，
分散：s² ＝ {(0－4.052)2*185＋(1－4.052)2*1149＋ …＋(12－4.052)2*0}/26306＝2.696．
　1つのサイコロが5か6の目を出す確率を p とすると，12個のサイコロを同時に振って，5か6の目が出る個数は 2 項分布 B(12, p) に従うはずである．
2 項分布 B(12, p) の平均は 12p，分散は 12p(1－p) であるので，p のデータからの推定値p^{^}は，データの平均値を用いて，

12p^{^}＝4.052， p^{^}＝4.052/12＝0.338

と推定される．このときの分散は，

分散＝12p^{^}(1－p^{^})＝ 12*0.338*(1－0.338)＝2.685

この分散はデータが2項分布 B(12, 0.338) に従っていれば取る値である．これがデータの分散と大きく違わないので，データはほぼ2項分布に従っていると判断される．
　なお，サイコロが完全に正しければ，12個のサイコロを同時に振って，5か6の目が出る個数は 2 項分布 B(12, 1/3) に従うはずである．この仮定のもとでは，平均と分散はそれぞれ，

平均＝12*1/3＝4，分散＝12*1/3*(1－1/3)＝4*2/3＝2.667

となるはずである．データの統計量とは異なる度合いが大きいので，サイコロは完全に正しくなく，大きな目（5か6の目）の出る確率の方がほんの少し高いと言える．
　以上を表にまとめると，

   モデル       平均       分散

   データ       4.052       2.696

   2項分布（p＝0.338）       4.052       2.685

   2項分布（p＝0.333）       4       2.667

モデル	平均	分散
データ	4.052	2.696
2項分布（p＝0.338）	4.052	2.685
2項分布（p＝0.333）	4	2.667

となる．なお，サイコロが本当に正しくないのかどうかの検定は，後で詳しく行う．
各モデルのもとでの確率分布のグラフは，

であり，データと 2項分布（p＝0.338）の確率モデルがよく適合しているようにみえる．

Weldon のサイコロ実験の R スクリプト
x <- 0:12 #個数　
dice <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18,0,0) #回数データ　
sum(dice) #試行回数　
pdice <- dice/sum(dice) #回数の確率　
m <- sum(x*pdice) #平均　
p <- m/12 #5，6の出る確率　
s2 <- sum(pdice*(x-m)^2) #分散　
v <- 12*p*(1-p) #二項分布のもとでの分散　
h1 <- dbinom(x, 12, 1/3) #正しいサイコロのもとでの二項確率分布　
h2 <- dbinom(x, 12, p) #推定確率からの二項確率分布　
dicedis <- rbind(pdice,h2,h1) #行ベクトル－＞行列　
colnames(dicedis) <- as.character(0:12) #列の名前　
barplot(dicedis, beside=TRUE, cex.axis=0.8, cex.lab=1.0, xlab="5,6の個数", ylab="確率", legend=c("データ","p=0.338", "p=0.333"))
title(main="Weldon のサイコロ実験の分布") #グラフタイトル　

**Weldon のサイコロ実験の R スクリプト**
x <- 0:12	#個数
dice <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18,0,0)	#回数データ
sum(dice)	#試行回数
pdice <- dice/sum(dice)	#回数の確率
m <- sum(x*pdice)	#平均
p <- m/12	#5，6の出る確率
s2 <- sum(pdice*(x-m)^2)	#分散
v <- 12p(1-p)	#二項分布のもとでの分散
h1 <- dbinom(x, 12, 1/3)	#正しいサイコロのもとでの二項確率分布
h2 <- dbinom(x, 12, p)	#推定確率からの二項確率分布
dicedis <- rbind(pdice,h2,h1)	#行ベクトル－＞行列
colnames(dicedis) <- as.character(0:12)	#列の名前
barplot(dicedis, beside=TRUE, cex.axis=0.8, cex.lab=1.0, xlab="5,6の個数", ylab="確率", legend=c("データ","p=0.338", "p=0.333"))
title(main="Weldon のサイコロ実験の分布")	#グラフタイトル

12人の兄弟中の女児数のデータ
19世紀末のドイツの病院のデータによると，同じ両親で12人きょうだいがいる6155家族の女児数の数は以下のようであった．

   女児数       0    1    2    3    4    5    6

度数    7       45       181       478       829       1112       1343

   女児数       7    8    9    10    11    12    合計

度数    1033       670       286       104       24       3    6155

女児数	0	1	2	3	4	5	6
度数	7	45	181	478	829	1112	1343
女児数	7	8	9	10	11	12	合計
度数	1033	670	286	104	24	3	6155

課題：「Weldon のサイコロ実験」と同様な解析を「12人のきょうだい中の女児数のデータ」に対して行え．

幾何分布（geometric distribution)

　成功確率を p とし，0 以上の整数 X に対し，確率密度が，

(4. 2)

となる分布．

　下の負の二項分布で n = 1 とおけば，
m(t) = p/(1 - qe^t)，q = 1 - p
E[X ] = q/p = (1 - p)/p，Var[X ] = q/p² = (1 - p)/p²

　成功確率 p のベルヌイ試行において，最初の成功が起こるまでの失敗の回数の分布．すなわち，離散的時間を考えた場合，初めて成功するまでの待ち時間の分布．

１の目が出るまで振るサイコロの回数の分布
　サイコロを振って，１の目が出るまで振り続けるとする．１の目が出るまで無駄に振り続ける回数は，成功確率 1/6 の幾何分布に従い，その確率密度のグラフは以下のようになる．

**幾何分布のグラフの R スクリプト**
x <- 0:100	#本当は無限大まで必要（永遠に１が出ない）
p <- 1/6	#成功確率
y <- dgeom(x, p)	#幾何分布の確率密度）
# x は 0 から 20 まで表示
plot(x, y, type="h", cex.lab=0.8, xlim=c(0,20), xlab="サイコロを無駄に振った回数", ylab="確率密度")
title(main="1 の目が出るまで振る回数の分布")	#タイトル
sum(x*y)	#平均（(1 - p)/p = (5/6)/(1/6) = 5）
sum(y*(x - (1-p)/p)^2)	#分散（(1 - p)/p² = (5/6)/(1/36) = 30）

超幾何分布（hypergeometric distribution)

　m 個の白石と n 個の黒石が入った袋から k 個の石を無作為に取り出したとき，白石の個数 X の従う確率密度は，

(4. 3)

で与えられる．
平均：E[X ] = kp，分散：Var[X ] = kp(1 - p)(m + n - k)，p = m/(m + n)

課題：超幾何分布の確率密度は dhyper(x, m, n, k)　で与えられる．m = 50，n = 450，k = 50，のときの上図のような超幾何分布の確率密度のグラフを描け．

復元抽出と非復元抽出
　超幾何分布は，標本抽出においての非復元抽出（sampling without replacement）に対する確率分布である．すなわち，m 個の白石と n 個の黒石が入った袋から石を１つ取り出し，それを戻さないで k 回繰り返したときの白石の個数 X の分布が超幾何分布になる．
　これに対し，石を取り出したらそれを袋に戻してから石をまた取り出す場合を復元抽出（sampling with replacement）という．このときは，白石が選ばれる確率が p_w = m/(m + n) と一定になるので，これを k 回繰り返したときの白石の個数 X の分布は２項分布 B(k, p_w) に従う．
　これより，石の個数 m + n が十分大きければ，超幾何分布は２項分布に近づくことがわかる．

**超幾何分布と二項分布の R スクリプト**
x <- 0:20	# 抽出された白石の個数
m <- 50	# 白石の個数
n <- 450	# 黒石の個数
k <- 50	# 抽出回数
y <- dhyper(x, m, n, k)	# 超幾何分布確率分布
plot(x,y,type="h")	#
p <- m/(m+n)	# 白石が選ばれる確率
yd <- dbinom(x, k, p)	# 二項確率
x1 <- x+0.1	# 0.1 ずらして表示
points(x1, yd, type="h", col="red")	# 二項確率のグラフ
title(main="超幾何分布と二項分布")	#
legend(locator(1),c("超幾何分布","二項分布"), lty=1, col=c("black", "red"))

生息数の推定（捕獲再捕獲法）
　超幾何分布は，生態学で動物の生息数を推定するのに用いられる．m 匹の動物を捕獲して目印をつけたのち放す．十分時間が経過して，捕獲した動物とそうでない動物が混ざり合ったと考えられてから，k 匹を再捕獲する．再捕獲で捕まった目印つきの個体数を X とすると，対象となっている動物の生息数を m + n と考えると，X は超幾何分布に従う．これを捕獲再捕獲法（capture - recapture method），もしくは，標識・再捕獲法（mark - recapture method）という．
　いま，再捕獲した k 匹のうち，x 匹が目印をつけたものだったとする．つまり，超幾何分布に従う確率変数 X の実現値が x であった（X = x と表記）とする．母集団全体での目印つき動物の頻度は m/(m + n) なので，再捕獲したサンプル中でも同じ割合で目印つきの動物を捕まえる可能性が高いと考えるのが，最も単純な考え方である．これより，動物の生息数 m + n は，
と推定される．これを Petersen 法という．

捕獲再捕獲法の精度（シミュレーション）
　課題の状況では，生息数500頭の生物種のうち 10 ％の m = 50 頭に標識をつけ，k = 50 頭を再捕獲して標識した頭数 X の分布を求めた．これをみると，期待される X = 5 が最頻値（モード）になっているが，X は， 0 から 10 くらいまでばらついていた．このため，X の実現値 x により生息数の推定値 km/x が変動する．この変動の程度は，シミュレーションにより求めることができる．すなわち，パラメータ m，n，k，の超幾何分布に従う乱数を発生させ，各乱数ごとに得られる生息数の推定値の分布をみればだいたいの様子がわかる． 10000回のシミュレーションを行ったところ，44回は標識のついた生物が再捕獲されなかったので，生息数の推定値は無限大に発散してしまった．このため，推定値分布の平均や分散を求めることはできなかったが，メディアンや分位点についての考察することはできる．
　それによると，推定値分布のメディアンは500であり，推定値分布の95％は312.5から1250の間に入っているので，このような調査はそれなりに現実を推定できることがわかる．

**生息数推定値分布の R スクリプト**
m <- 50	# 標識をつけた数
n <- 450	# 標識をつけられていない数
k <- 50	# 再捕獲数
estimate <- NULL	# 個体数推定値の定義
for (i in 1:10000) {	# 10000回のシミュレーション
x <- rhyper(1, m, n, k)	# 超幾何分布に従う乱数１つ抽出
estimate <- c(estimate, k * m / x)	# 個体数推定値列ベクトル
}	#
hist(estimate, breaks=seq(0, 2600, by=100), main="")	#推定値のヒストグラム
title(main="生息数の推定値の分布")	#タイトル
table(estimate)	#推定値の階級分け
ord <- order(estimate)	#推定値の順位
estimate[ord[500]]	#小さい方から500番目（5％点）
estimate[ord[1000]]	#小さい方から1000番目（10％点）
median(estimate)	#メディアン（中央値）
estimate[ord[9000]]	#小さい方から9000番目（90％点）
estimate[ord[9500]]	#小さい方から9500番目（95％点）

課題：上の例では捕獲（m = 50），再捕獲（k = 50）合わせて100頭を捕獲している．捕獲数を半分にした場合の生息数の推定精度はどうなるか．

ポアソン分布（Poisson distribution）

　正のパラメータ λ と，0 以上の整数 X に対し，確率密度が

(4. 4)

となる分布．

poisson

二項分布の極限分布としてのポアソン分布
　いま，単位時間を n 等分して，この時間間隔である事象が１回生起する確率を p とする．ここで，

np = λ

とおいて λ の値を固定して n を大きくして時間間隔を小さくしていくと，生起確率 p も小さくなるので，１つの時間間隔で，事象が２回起こる確率は十分小さくなり，無視できるとする．すると，単位時間内（n　回の試行）で事象が x 回生起する確率は成功確率 p の二項分布に従い，単位時間内での平均生起回数は，np = λ となる．これより，

となり，n の極限でこの二項分布は平均 λ のポアソン分布に従うことがわかる．

二項分布がポアソン分布に近づく，の R スクリプト
x <- 0:8 #グラフのx軸の範囲　
lam <- 2 #λの定義　
yp <- dpois(x,lam) #ポアソン分布の確率密度　
y1 <- dbinom(x, 5, 0.4) #二項分布（n = 5，p = 0.4）の確率密度　
y2 <- dbinom(x, 10, 0.2) #二項分布（n = 10，p = 0.2）の確率密度　
y3 <- dbinom(x, 20, 1/10) #二項分布（n = 20，p = 0.1）の確率密度　
y4 <- dbinom(x, 40, 1/20) #二項分布（n = 40，p = 0.05）の確率密度　
plot(x, y1, type="b", ylab="確率") #二項分布（n = 5，p = 0.4）のプロット（黒）　
points(x, y2, type="b", col="green") #二項分布（n = 10，p = 0.2）のプロット（赤）
points(x, y3, type="b", col="blue") #二項分布（n = 20，p = 0.1）のプロット（青）
points(x, y4, type="b", col="purple") #二項分布（n = 40，p = 0.05）のプロット（紫）
points(x, yp, type="b", col="red") #ポアソン分布のプロット（赤）
title(main="二項分布（np = 2）がポアソン分布に近づく様子")
# 凡例の記述（locator(1)は，凡例の記述場所をクリックで指定）
legend(locator(1), c("p=0.4", "p=0.2", "p=0.1", "p=0.05", "ポアソン"),
lty=1, col=c("black", "green", "blue", "purple", "red"))

**二項分布がポアソン分布に近づく，の R スクリプト**
x <- 0:8	#グラフのx軸の範囲
lam <- 2	#λの定義
yp <- dpois(x,lam)	#ポアソン分布の確率密度
y1 <- dbinom(x, 5, 0.4)	#二項分布（n = 5，p = 0.4）の確率密度
y2 <- dbinom(x, 10, 0.2)	#二項分布（n = 10，p = 0.2）の確率密度
y3 <- dbinom(x, 20, 1/10)	#二項分布（n = 20，p = 0.1）の確率密度
y4 <- dbinom(x, 40, 1/20)	#二項分布（n = 40，p = 0.05）の確率密度
plot(x, y1, type="b", ylab="確率")	#二項分布（n = 5，p = 0.4）のプロット（黒）
points(x, y2, type="b", col="green")	#二項分布（n = 10，p = 0.2）のプロット（赤）
points(x, y3, type="b", col="blue")	#二項分布（n = 20，p = 0.1）のプロット（青）
points(x, y4, type="b", col="purple")	#二項分布（n = 40，p = 0.05）のプロット（紫）
points(x, yp, type="b", col="red")	#ポアソン分布のプロット（赤）
title(main="二項分布（np = 2）がポアソン分布に近づく様子")
# 凡例の記述（locator(1)は，凡例の記述場所をクリックで指定）
legend(locator(1), c("p=0.4", "p=0.2", "p=0.1", "p=0.05", "ポアソン"),
lty=1, col=c("black", "green", "blue", "purple", "red"))

死亡記事件数
　ポアソン分布は，稀な事象の生起モデル（自動車事故，機械の故障，DNAの塩基置換など）に用いられる．下の表は，ロンドンの新聞記事（1096日間）に載った１日あたりの85才以上の死亡記事の件数である．（L. Whitaker "On Poisson's law of small numbers" Biometrika, Vol. 10, p36- , 1914）

死亡記事件数 0 1 2 3 4 5 6 以上　

日数 484 391 164 45 11 1 0

死亡記事件数	0	1	2	3	4	5	6 以上
日数	484	391	164	45	11	1	0

　死亡記事件数データの平均は0.8239，分散は0.8294，であった．このデータがポアソン分布に従っていると考える．ポアソン分布の平均は λ なので，λ = 0.8239 のポアソン分布にあてはめてみたところ，非常によく一致していた．
　また，ポアソン分布は，平均と分散が等しいという特徴がある．データの平均と分散の値が近いことから，データはポアソン分布によく適合していることを示している．

Poisson
データとモデルとの統計量の比較
   モデル    平均    分散
   死亡記事件数    0.8239    0.8294
   ポアソン分布    0.8239    0.8239

死亡記事件数の R スクリプト
x <- 0:6 #グラフのx軸の範囲　
y <- c(484, 391, 164, 45, 11, 1, 0) #死亡記事件数データ　
s <- sum(y) #データ総数　
m <- sum(x*y/s) #データ分布の平均　
v <- sum((x-m)^2*y/s) #データ分布の分散　
yp <- dpois(x, m) #平均 m のポアソン分布確率密度　
plot(x, y/s, type="h", ylab="確率") #データの棒グラフ表示　
points(x, yp, type="b", col="red") #ポアソン分布の重ねがき（赤）　
title(main="死亡記事件数へのポアソン分布のあてはめ")
legend(3.5, 0.4, c("データ", "ポアソン分布"), lty=1, col=c("black","red"))

**死亡記事件数の R スクリプト**
x <- 0:6	#グラフのx軸の範囲
y <- c(484, 391, 164, 45, 11, 1, 0)	#死亡記事件数データ
s <- sum(y)	#データ総数
m <- sum(x*y/s)	#データ分布の平均
v <- sum((x-m)^2*y/s)	#データ分布の分散
yp <- dpois(x, m)	#平均 m のポアソン分布確率密度
plot(x, y/s, type="h", ylab="確率")	#データの棒グラフ表示
points(x, yp, type="b", col="red")	#ポアソン分布の重ねがき（赤）
title(main="死亡記事件数へのポアソン分布のあてはめ")
legend(3.5, 0.4, c("データ", "ポアソン分布"), lty=1, col=c("black","red"))

商品在庫
　ある商店において，ある商品の１日で売れる数は，平均 2 のポアソン分布にほぼ従っていた．１週間で商品の在庫を切らさない確率を95％以上にするためには，少なくとも在庫をいくつ抱えておく必要があるか．

　7 日間で売れる商品の個数は，平均 2×7 = 14 のポアソン分布に従うと考えられる．ポアソン分布の95％（分位）点を求めればよい．これより，必要な在庫は20個とわかる．

**商品在庫の R スクリプト**
m <- 14	# ポアソン分布のパラメータ（平均）
x <- 0:30	# グラフのx軸の範囲
yp <- dpois(x, m)	# ポアソン分布確率密度
cyp <- ppois(x,m)	# ポアソン分布累積確率
stok <- qpois(0.95, m)	# 95％（分位）点
stok	# 答えの表示
op <- par(mfrow = c(1, 2))	# 横に２つのグラフを並べる
plot(x, yp, type="h", ylab="確率密度")	# 確率密度グラフ
points(x[2:21],yp[2:21], type="h", col="red")	# x = 20 まで赤色表示
plot(x,cyp,type="s", ylab="累積確率")	# 累積確率グラフ
arrows(stok,0.95,stok,0, length=0.1, col="red")	# 赤矢印
segments(0,0.95, stok,0.95, col="red")	#
par(op)	# グラフ表示もとに戻す
title(main="ポアソン分布（λ = 14）の95％点")	# グラフタイトル

生物個体の分布シミュレーション
　ある領域内である生物種が一様に分布しているとする．この領域をメッシュで区切ると，メッシュ内で観測される生物の個体数はポアソン分布に従う．このような点の配置は，ポアソン配置（Poison configuration）と呼ばれている．ポアソン分布は平均と分散が等しいことから，カウント数分布の平均と分散の値を計算すれば，ポアソン分布に従いそうかがわかる．
　点の配置パターンは，なわばりを持つような鳥の巣の配置や，コンビニなどの店舗配置パターンなどに応用される．

**生物個体分布シミュレーションの R スクリプト**
n <- 200	# 個体数
m <- 10	# メッシュ（m² 個）
x <- runif(n)	# 一様乱数n個
y <- runif(n)	# 一様乱数n個
#	#
count <- NULL	# count の定義
for(i in 1:m){	# m 回の繰り返し
n1 <- (1:n)[x < (i-1)/m]	# (i-1)/m 以下の乱数である番号
n2 <- (1:n)[x < i/m]	# i/m 以下の乱数である番号
nin <- n2[!n2 %in% n1]	# (i-1)/m から i/m の番号
yy <- y[nin]	# 上記 x 座標に対する y 座標
a <- hist(yy, breaks=0:m/m)	# yy を 0 から 1 まで 1/m きざみで区切る
count <- c(count, a$counts)	# 区切った領域に入った個体の個数のベクトル
}	# 繰り返しここまで
mc <- max(count)	# メッシュ内の個数の最大値
xp <- 0:mc	# 個数の定義域
d <- factor(count, levels=xp)	# 個数が 0 の階級も含める
table(d)	# メッシュ内個数の区分
s <- sum(table(d))	#　総個数
m1 <- sum(xp*table(d)/s)	# カウント分布の平均
m2 <- mean(count)	# カウントデータの標本平均
v1 <- sum(table(d)/s*(xp-m1)^2)	# カウント分布の分散
v2 <- var(count)	# カウントデータの標本分散
#	#
op <- par(mfrow = c(1, 2))	# 横に２つのグラフを並べる
plot(x,y, col="red")	# 個体分布の表示
abline(h=0, v=0)	# 外枠
abline(h=1, v=1)	# 外枠
for(i in 1:m) abline(h=i/m, v=i/m, lty=2)	# メッシュ区分線
plot(xp, table(d)/s, type="h")	# 区分された分布のグラフ
lam <- n/(m*m)	# 平均
yp <- dpois(xp, lam)	# ポアソン分布確率密度
points(xp, yp, type="b", col="red")	# ポアソン分布グラフ表示
par(op)	# グラフ表示もとに戻す
title(main="区画内個体数へのポアソン分布のあてはめ（平均：2）")	# タイトル
m1	# カウント分布平均
v1	# カウント分布分散

負の二項分布（negative binomial distribution）

　成功確率が p で，正のパラメータ n をもち，0 以上の整数 X に対し，確率密度が

(4. 5)

となる分布．Γ関数で確率密度を定義すれば，n は必ずしも整数である必要はない．

negb

　成功率 p のベルヌイ試行において，n 回の成功が起こるまでの失敗の回数の分布で，幾何分布は負の二項分布で n = 1 とおいたものである．

虫歯数データ
　ある小学校の一人あたりの虫歯の数のデータ（鳥居敏雄ら　「医学・生物学のための推計学」　東京大学出版会）

虫歯の数 0 1 2 3 4 5 6 7 8

児童の数 4 9 16 13 9 7 5 4 3

虫歯の数	0	1	2	3	4	5	6	7	8
児童の数	4	9	16	13	9	7	5	4	3

　虫歯数データが負の二項分布に従っていると仮定する．虫歯データの平均は3.3286，分散は4.3063，であった．負の二項分布は２つのパラメータ n，p を持つのでそれをデータから推定する必要がある．分布パラメータの推定法は後に詳しく議論するが，ここでは直感的に理解しやすいモーメント法を説明する．これは，分布の平均や分散などのモーメントがデータの平均や分散と等しい，とおくことでパラメータの推定を行う方法である．

課題：負の二項分布の確率密度は dnbinom(x, n, p) で与えられる．モーメント法を用いてパラメータ n，p の値を推定し，この推定値を用いて，児童の虫歯データに負の二項分布をあてはめたグラフ（下図のようなもの）を描け．

　モーメント法による推定値を用いて，虫歯数データに負の二項分布をあてはめた．データは分布にまあまあ適合しているようにみえる．

負の二項分布パラメータ推定（モーメント法）の R スクリプトの一部
x <- 0:10 #グラフのx軸の範囲　
y <- c(4,9,16,13,9,7,5,4,3,0,0) #虫歯数　
s <- sum(y) #虫歯数の総和　
m <- sum(x*y/s) #平均　
v <- sum((x-m)^2*y/s) #分散　

**負の二項分布パラメータ推定（モーメント法）の R スクリプトの一部**
x <- 0:10	#グラフのx軸の範囲
y <- c(4,9,16,13,9,7,5,4,3,0,0)	#虫歯数
s <- sum(y)	#虫歯数の総和
m <- sum(x*y/s)	#平均
v <- sum((x-m)^2*y/s)	#分散

生物個体の分布２シミュレーション
　前節のポアソン配置は，個体間が互いにまったく影響を受けずにランダムに配置している場合であった．個体の配置が他の個体の配置に影響を与える配置として，個体がグループをなしている場合を考える．
　親個体をランダムに配置させる（下図の緑十字）．子ども個体の数はポアソン分布に従い，親を中心とした正規分布（後述）に従って配置させるとする．こうすると，親個体のまわりに子ども個体が集まって分布するようになるので，パッチをつくったような配置ができる．先ほどと同様にメッシュで区切り，その中の子ども個体数をカウントした分布は，分散が平均より大きくなるので，ポアソン分布にはあてはまらず，負の二項分布の方があてはまりがよい．

**生物個体分布２シミュレーションの R スクリプト**
n <- 200	# 個体数（予定）
m <- 10	# メッシュ（m² 個）
p <- 4	# 平均子ども数
sig <- 0.05	# 正規分布標準偏差
xx <- NULL	# xx（子ども座標）の定義
xx0 <- NULL	# xx0（親座標）の定義
np <- round(n/p)	# 親の数（round() は 5 捨 6 入）
for(i in 1:np){	# np 回の繰り返し
x0 <- runif(2)	# 親個体の座標を一様乱数で生成
n0 <- rpois(1, p)	# 子どもの数を平均 p のポアソン乱数で生成
for(j in 1:n0){	# n0 回の繰り返し
xd <- rnorm(2, m=x0, sd=sig)	# 子どもの座標を，正規乱数 N（x0, sig²）で生成
if(xd[1] > 1) xd[1] <- xd[1] - 1	# x 座標が 1 を超えたとき区画の左端に
if(xd[1] < 0) xd[1] <- xd[1] + 1	# x 座標が 0 未満のとき区画の右端に
if(xd[2] > 1) xd[2] <- xd[2] - 1	# y 座標が 1 を超えたとき区画の下辺に
if(xd[2] < 0) xd[2] <- xd[2] + 1	# y 座標が 0 未満のとき区画の上辺に
xx <- rbind(xx, xd)	# 個体座標行列の行の追加
xx0 <- rbind(xx0, x0)	# 個体座標行列の行の追加
}	#
}	#
# 区画内個体数	#
x <- xx[,1]; y <- xx[,2]	# x 座標ベクトル，y 座標ベクトル
count <- NULL	# count の定義
for(i in 1:m){	# m 回の繰り返し
n1 <- (1:n)[x < (i-1)/m]	# (i-1)/m 以下の乱数である番号
n2 <- (1:n)[x < i/m]	# i/m 以下の乱数である番号
nin <- n2[!n2 %in% n1]	# (i-1)/m から i/m の番号
yy <- y[nin]	# 上記 x 座標に対する y 座標
a <- hist(yy, breaks=0:m/m)	# yy を 0 から 1 まで 1/m きざみで区切る
count <- c(count, a$counts)	# 区切った領域に入った個体の個数のベクトル
}	#
mc <- max(count)	# メッシュ内の個数の最大値
xp <- 0:mc	# 個数の定義域
d <- factor(count, levels=xp)	# 個数が 0 の階級も含める
table(d)	# メッシュ内個数の階級区分
s <- sum(table(d))	# 総個体数
m1 <- sum(xp*table(d)/s)	# カウント分布の平均
v1 <- sum(table(d)/s*(xp-m1)^2)	# カウント分布の分散
#	#
op <- par(mfrow = c(1, 2))	# 横に２つのグラフを並べる
plot(x,y, col="red")	# 個体分布の表示
points(xx0, pch="+", col="green")	# 個体分布の表示
abline(h=0, v=0)	# 外枠
abline(h=1, v=1)	# 外枠
for(i in 1:m) abline(h=i/m, v=i/m, lty=3)	# メッシュ区分線
plot(xp, table(d)/s, type="h", ylim=c(0,0.35))	# 区分された分布のグラフ
yp <- dpois(xp, m1)	# ポアソン分布確率
points(xp, yp, type="b", col="red")	# ポアソン分布表示
nbp <- m1/v1	#
nbn <- m1*nbp/(1-nbp)	#
ynb <- dnbinom(xp, nbn, nbp)	# 負の二項分布確率
points(xp, ynb, type="b", col="green")	# 負の二項分布表示
legend(2.5, 0.33, c("データ", "ポアソン分布", "負の二項分布"), lty=1, col=c("black","red","green"))
par(op)	# 個体数
title(main="区画内個数へのポアソン分布と負の二項分布のあてはめ")
m1	# 平均
v1	# 分散

課題：平均子ども数である p の値を増やすとどうなるか．

離散分布あてはめのまとめ

　記事数や虫歯数などのカウントデータに離散分布をあてはめるとき留意する点は，データの平均と分散の大きさを比較することである．

データ適応する分布   平均     分散

  平均＞分散     二項分布 np np(1 - p)

  平均≒分散     ポアソン分布 λ λ

  平均＜分散     負の二項分布      n(1 - p)/p       n(1 - p)/p²

データ	適応する分布	平均	分散
平均＞分散	二項分布	np	np(1 - p)
平均≒分散	ポアソン分布	λ	λ
平均＜分散	負の二項分布	n(1 - p)/p	n(1 - p)/p²

4-2．連続型分布の例

一様分布（uniform distribution)

　2 つのパラメータ a，b（a＜b）をもつ確率密度関数が

(4. 6)

で表される分布．

poisson

（0, 1）一様乱数（uniform random number）
　（0, 1）区間の一様乱数はパソコンに標準装備されているが，周期性があるなどの擬似乱数になっている． R の擬似乱数はかなり性質がよいことが知られている．runif(n) で n 個の（0，1）一様乱数が生成される．

一様乱数のRスクリプト
x <- runif(10000) #（0，1）一様乱数1000個列
hist(x, main="(0, 1) 一様乱数 10000個") #ヒストグラム表示

一様乱数のRスクリプト
x <- runif(10000)	#（0，1）一様乱数1000個列
hist(x, main="(0, 1) 一様乱数 10000個")	#ヒストグラム表示

π の値の推定
　一様乱数を用いて π の値を推定することができる．すなわち，区間（-1，1）の一様乱数 2 個で， -1＜x＜1，-1＜y＜1 の正方形内の１点が定義できる．この点を多数生成させ，その内 x² + y² ＜ 1，を満たす点の個数の割合を計算すれば π の値が求まる．

**一様乱数による π の近似の R スクリプト**
n <- 10000	#一様乱数の個数
x <- runif(n, -1, 1)	#(-1, 1) の範囲の一様乱数 n 個生成
y <- runif(n, -1, 1)	#
r <- x^2 + y^2	#原点からの距離の２乗
plot(x,y, type="n", xlim=c(-1,1), ylim=c(-1,1))	#グラフの表示範囲の指定
abline(h=0)	# x 軸の表示
abline(v=0)	# y 軸の表示
segments(-1, 1, 1, 1)	#(-1, 1)から(1, 1)までの直線
segments(1, 1, 1, -1)	#
segments(-1, -1, 1, -1)	#
segments(-1, -1, -1, 1)	#
pin <- (1:n)[r<1]	#乱数のうち単位円内に入る乱数の番号
points(x[-pin], y[-pin], pch=".", col="green")	#単位円の外の乱数を緑点で表示
points(x[pin], y[pin], pch=".", col="red")	#単位円内の乱数を赤い点で表示
s <- 0:360	# 0 度から 360 度
theta <- s*pi/180	#度をラジアンに変換
xp = sin(theta)	#単位円の x 座標
yp = cos(theta)	#単位円の y 座標
points(xp,yp, type="l")	#単位円を表示
title(main="（-1，1）一様乱数による点列と単位円")
length(pin)	#単位円内に入った乱数の個数

課題：単位円内に落ちた乱数の個数から π の近似値を求めよ．また，π の近似の精度を上げて π の近似値を再計算せよ．

β（ベータ）分布（beta distribution)

　2 つの正のパラメータ a，b をもつ確率密度関数が

(4. 7)

で表される分布．

であるので，

パラメータによる形状の違い


	a	b
－	0.5	0.5
－	1	1
－	2	2
－	1	3
－	2	4

β分布密度関数の R スクリプト
x <- seq(0,1, by=0.01) # x の定義 0 から 1 まで 0.01 きざみ　
y <- dbeta(x, 0.5, 0.5) # a, b = 0.5, 0.5 のβ分布
plot(x, y, type="l", ylim=c(0,3), col="red") # y を 0 から 3 に指定（赤）　
abline(v=0, h=0) # y 軸と x 軸の表示　
curve(dbeta(x, 1, 1), 0, 1, add=T) # a, b = 1, 1 のβ分布　
curve(dbeta(x, 2, 2), 0, 1, add=T, col="blue") # a, b = 2, 2 のβ分布（青）　
curve(dbeta(x, 1, 3), 0, 1, add=T, col="green") # a, b = 1, 3 のβ分布（緑）　
curve(dbeta(x, 2, 4), 0, 1, add=T, col="purple") # a, b = 2, 4 のβ分布（紫）　
title(main="β分布") #　

指数分布（exponential distribution)

　正のパラメータ λ をもつ確率密度関数 f(x) が

(4. 8)

で表される分布．
　下のガンマ分布で，a = 1，s = 1/λ，とおいて，

指数分布密度関数の R スクリプト
curve(dexp(x, 1), 0, 10) #指数分布の密度関数　
abline(v=0, h=0) # y 軸，x 軸表示　
title(main="指数分布") #　

Γ（ガンマ）分布（gamma distribution)

　正の形状（シェープ）パラメータ a，正のスケールパラメータ s をもつ確率密度関数が

(4. 9)

で表される分布．形状パラメータが a = 1 で，λ = 1/s とおくと指数分布に一致する．

パラメータによる形状の違い

a s
－ 1 1 指数分布
－ 2 1
－ 2 2
－ 4 1

課題： Γ 分布の確率密度は dgamma(x, a, 1/s)，もしくは，dgamma(x, shape=a, scale=s) で与えられる．Γ 分布のパラメータを変えて，上の図のようなグラフを描け．

伝染分布（contagious distribution）
　ある分布のパラメータが，ある分布に従っているとき伝染分布という．ポアソン分布のパラメータ θ がガンマ分布（シェープパラメータ r，レイトパラメータ λ）に従っているとき，その分布は負の二項分布に従う．すなわち，

であるとすると，伝染分布は，

となる．シェープパラメータ r = 2，レイトパラメータ λ = 1，としたガンマ分布から乱数を発生させ，生成した乱数 θ をパラメータにもつポアソン分布から乱数を生成させる．このような乱数を多数（n = 1000）生成させ，得られたカウント数の分布に理論上考えられる負の二項分布をあてはめてみた

**伝染分布シミュレーションの R スクリプト**
n <- 1000	# 伝染分布乱数の個数
count <- NULL	# count の定義
for(i in 1:n){	# n 回の繰り返し
r <- 2	# ガンマ分布のシェープパラメータ
lam <- 1	# ガンマ分布のレイトパラメータ
theta <- rgamma(1, r, lam)	# ガンマ分布乱数１つ生成
count <- c(count, rpois(1, theta))	# 生成乱数をパラメータにするポアソン乱数列
}	#
mc <- max(count)	# ポアソン・ガンマ乱数列の最大値
xp <- 0:mc	# 個数の定義域
d <- factor(count, levels=xp)	# 個数が 0 の階級も含める
table(d)	# メッシュ内個数の階級区分
s <- sum(table(d))	# 総個体数
m1 <- sum(xp*table(d)/s)	# カウント分布の平均
v1 <- sum(table(d)/s*(xp-m1)^2)	# カウント分布の分散
plot(xp, table(d)/s, type="h", ylim=c(0,0.3), ylab="確率")	# ポアソン・ガンマ乱数のグラフ
p <- lam/(lam+1)	# 負の二項分布確率パラメータ
ynb <- dnbinom(xp, r, p)	# 負の二項分布確率密度
points(xp, ynb, type="b", col="red")	# 負の二項分布のグラフ
title(main="ポアソン・ガンマ伝染分布乱数への\n負の二項分布のあてはめ", cex.main=0.8)

コーシー分布（Cauchy distribution)

　位置パラメータ a と正のスケールパラメータ s をもつ確率密度関数 f(x) が

(4.10)

で表される分布．積率母関数が計算できないので，平均も分散も存在しない．

コーシー分布密度関数の R スクリプト
curve(dcauchy(x), -5, 5) #コーシー分布の密度関数　
abline(v=0, h=0) # y 軸，x 軸表示　
title(main="コーシー分布") #　

正規分布（normal distribution)

　平均 μ，分散 σ² の２つのパラメータをもつ確率密度関数が

(4.11)

で表される分布で，N（μ，σ²）と表記する．μ は位置パラメータ（location parameter）で，スケールパラメータ σ を標準偏差（standard deviation）という．

パラメータによる形状の違い

**正規分布密度関数の R スクリプト**
# 平均の異なる正規分布
curve(dnorm(x, 40, 4), 30, 70, ylim=c(0,0.2), xlab="",ylab="確率密度")	# 平均：40，標準偏差：4
curve(dnorm(x, 50, 4), add=TRUE, col="red")	# 平均：50，標準偏差：4
curve(dnorm(x, 60, 4), add=TRUE, col="blue")	# 平均：60，標準偏差：4
title(main="正規分布（異なる平均(μ)）\n平均＝40, 50, 60，σ = 4")	# タイトル
legend(52,0.19,c("μ = 40","μ = 50","μ = 60"), lty=1, col=c("black","red","blue"))	# 凡例
# 標準偏差（分散）の異なる正規分布
curve(dnorm(x, 50, sd=2), 30, 70, ylim=c(0,0.2), xlab="",ylab="確率密度")	# 平均：50，標準偏差：2
curve(dnorm(x, mean=50, sd=4), add=TRUE, col="red")	# 平均：50，標準偏差：4
curve(dnorm(x, mean=50, sd=6), add=TRUE, col="green")	# 平均：50，標準偏差：6
title(main="正規分布（異なる標準偏差(σ)）\n平均＝50，σ = 2, 4, 6")	# タイトル
legend(55,0.19,c("σ = 2","σ = 4","σ = 6"), lty=1, col=c("black","red","green"))	# 凡例

標準正規分布 N（0，1）
平均 0，分散 1 の正規分布を標準正規分布（Standard normal distribution）という．
確率変数 X が平均 μ，分散 σ² の正規分布に従っている，すなわち， X ～ N（μ，σ²），であるとき，

(4.12)

と標準化すると，確率変数 Z は，標準正規分布 N（0，1）に従い，その確率密度関数は，φ(z)

(4.13)

で表され，累積分布関数は，Φ(z)

(4.14)

で表現される．標準正規分布では，-1 ≦ z ≦ 1 の範囲に全体の68.3％が含まれ， -2 ≦ z ≦ 2 の範囲に全体の95.4％が含まれる（下左図）．

**# 標準正規分位点の R スクリプト**
pnorm(-1)	# = 0.16（赤矢印）
pnorm(1) - pnorm(-1)	# = 0.683
pnorm(2) - pnorm(-2)	# = 0.954
qnorm(0.975)	# = 1.96（青矢印），両側 5 ％点

英国成人男子身長データの正規分布へのあてはめ
　身長や体重などの身体データや得点データなどは正規分布に従うことが多い．たとえば，英国成人男子身長（インチ）のデータに正規分布をあてはめてみる．

英国男子身長データ（インチ）
   身長       57    58    59    60    61    62    63    64    65    66

人数    2       4       14       41       83       169       394       669       990       1223

   67    68    69    70    71    72    73    74    75    76    77
   1329       1230       1063       646       392       202       79       32       16       5       2

まず，データから統計量を求める．データ総数は 8585 名である．
平均：x^- ＝ (57*2＋58*4＋59*14＋…＋77*2)/8585 ＝ 67.02，
分散：s² ＝ {(57－67.02)²*2＋(58－67.02)²*4 ＋ …＋(77－67.02)²*2}/8585 ＝ 6.62．
これより，平均 μ ＝ 67.02，分散 σ² ＝ 6.62，の正規分布にあてはめ，赤線でグラフ表示したところ，データ分布によく一致していた．

# 英国成人男子身長データの R スクリプト
x <- 57:77 # 身長（x）の範囲　
y <- c(2, 4, 14, 41, 83, 169, 394, 669, 990, 1223, 1329, # 身長ごとのデータ　
1230, 1063, 646, 392, 202, 79, 32, 16, 5, 2) # 　
s <- sum(y) # データ総数　
m <- sum(x*y/s); m # データの平均　
v <- sum(y/s*(x-m)^2); v # データの分散　
plot(x, y/s, type="h", xlab="身長（インチ）", ylab="頻度") # データの棒グラフ表示　
curve(dnorm(x, m, sqrt(v)), 57, 77, add=T, col="red") # 正規密度のグラフ表示　
title(main="身長データに対する正規分布のあてはめ") # タイトル　

シェパードの補正（Sheppard's correction）
　連続値データ（X）を階級に分けた場合（X_c），階級の中央値にデータを丸めた誤差（E）が発生する．すなわち，
X_c = X + E
これより，階級に丸めた値の分散は，
Var[X_c] = Var[X] + Var[E]
となり，元のデータより，階級データの方が丸め誤差の分だけ分散が大きい．
　いま，階級幅を h とすると，階級の中央値はその値から h/2 小さいデータから h/2 大きいデータまでが丸められる．この丸められた値は区間 h に一様に分布すると考えると，丸め誤差は幅 h の一様分布の分散 h²/12 になると考えられる．これより，元データの分散は，手元にある丸めた階級データの分散から，
Var[X] = Var[X_c] - h²/12
となり，実際のデータは少し分散が小さかったと推定される．
　先ほど取り上げた英国人男子身長データは，1 インチにデータを丸めたものなので，身長の本当の分散は，インチ単位で丸めた値から得られる分散より 1/12 小さいと推定される．すなわち，分散推定値は，
Var[X] = Var[X_c] - Var[E] = 6.617 - 0.083 = 6.534
となり，標準偏差の推定値は，√6.534 = 2.556，と推定される．

正規分布から得られる確率的推論
　英国人成人男子身長データは，平均 μ = 67.02，標準偏差 σ = 2.556 の正規分布，N(67.02, 6.534)，に従っていることがわかった．このことから，

英国人成人男子で 70 インチ以上の人の比率は．
英国人成人男子で背の高い方から 10％以内に入るには何インチ以上であればよいか．
英国人成人男子で，65 インチから 70 インチまでの間の比率は．

といった集団に対する確率的問に簡単に答えることができる．

# 正規分布から得られる確率的推論の R スクリプト
m <- 67.02; s <- 2.556 # 平均と標準偏差の指定
1 - pnorm(70, mean=m, sd=s) # 1. 70 までの累積確率を 1 から引く
qnorm(0.9, mean=m, sd=s) # 2. 累積確率が 0.9 となる身長
pnorm(70, mean=m, sd=s) - pnorm(65, mean=m, sd=s) # 3. （70 までの累積確率）－（60 までの累積確率）

二項分布が正規分布に近づく様子
　成功確率 p の二項分布は，試行回数 n を増やしていくと，平均 np，分散 np(1 - p) の正規分布に近づく．

# 二項分布が正規分布に近づく様子の R スクリプト
n <- 5 # 打数　
x <- 0:n # xの範囲　
p <- 0.3 # 打率　
hit <- dbinom(x, size=n, prob=0.3) # 二項確率　
y <- pbinom(x, size=n, prob=0.3) # 二項累積確率　
m <- n*p # 平均　
sd <- sqrt(n*p*(1-p)) # 標準偏差　
op <- par(mfrow = c(1, 2)) # 　
plot(x, hit, type="h", ylim=c(0,0.4), xlim=c(0,7), xlab="ヒット数", ylab="確率密度")
curve(dnorm(x, mean=m, sd=sd), add=TRUE, col="red") # 確率密度　
plot(x, y, type="s", ylim=c(0,1), xlim=c(0,7), xlab="ヒット数", ylab="累積確率")
curve(pnorm(x, mean=m, sd=sd), add=TRUE, col="red") # 累積確率　
par(op) # 　
title(main ="二項分布：n = 5 打数，打率 p = 0.3；正規分布：N(1.5, 1.05) ")

課題：打数 n を大きくして，二項分布が正規分布に近づく様子を確かめよ．

正規 Q - Q （Quantile - Quantile）プロット
　正規分布の分位点と標本（サンプル）の分位点との関係を２次元上にプロットしたもの，標本分布が正規分布に従っていれば直線上に分布する．直線からの「ずれ」で正規分布からの隔たりが視覚的に表現される．なお，直線を表示する関数 qqline() は，対応する正規分布とデータの 1/4 分位点（25％点）と 3/4 分位点（75％点）とを結んだ直線である．すなわち，データと正規分布との中央部分（四分位範囲）をそろえた場合の直線である．
　標準正規乱数 n = 1000 個を発生させ，そのヒストグラムに標準正規分布をあてはめたグラフと，正規分布との適合性をみるため正規 Q - Q プロットを表示した．これをみると，正規乱数は，正規分布にピッタリと適合していることがわかる．

**# 標準正規乱数の正規 Q - Q プロットの R スクリプト**
n <- 10000	# 乱数列の長さ
x <- rnorm(n)	# 標準正規乱数
op <- par(mfrow = c(1, 2))	#
hist(x, breaks=seq(-10,10, by=0.2), xlim=c(-5,5),freq=F, main="")	# 乱数のヒストグラム
curve(dnorm(x), add=T, col=2)	# 標準正規分布の重ね合わせ
title(main="正規乱数のヒストグラム")	# タイトル
qqnorm(x, xlab="正規分布分位点", ylab="正規乱数分位点", main="")	# 正規 Q - Q プロット
qqline(x, col=2)	# 正規分布の四分位範囲直線表示
title(main="正規 Q - Q プロット")	# タイトル
par(op)	#

　次に，英国人成人身長データと二項分布を正規分布にあてはめた場合について，正規分布との適合性を正規 Q - Q プロットでみてみよう．左下図は，身長データの正規 Q - Q プロットで，正規分布から少しずれている様子がわかる．右下図は，打率 p = 0.3 の選手の n = 100 打席でのヒット数の分布で，正規分布によくフィットしているのがよくわかる．

課題：英国人身長データと二項分布を正規分布にあてはめたときの正規 Q - Q プロットを描け．

中心極限定理（central limit theorem）
　X₁，X₂，…，X_n，を平均 μ，分散 σ² である分布からの無作為標本であるとすると，標本平均 X^-_n の分布は，サンプルサイズ n を大きくしていくと平均 E[X^-_n ] = μ，分散 Var[X^-_n ] = σ²/n の正規分布に近づく．
　中心極限定理により，平均と分散をもっていれば母集団の分布が何であっても，標本平均の分布はサンプルサイズを大きく（サンプル数を多く）すれば正規分布に従うので，正規分布に基づいた確率的推論を行ってもよいことを保証している．
　中心極限定理が成り立つ様子を一様分布でみてみる．一様乱数 2 個の標本平均 X^-₂ 10000個の分布は三角形型をしていて，正規分布とは似ていない．しかし，10個の標本平均 X^-₁₀ 10000個の分布は正規分布と近づいたが，尾（テイル）の部分のあてはまりは良くない．30個の標本平均 X^-₃₀ 10000個の分布をみると，尾の部分のあてはまりも改善されてくる．

**# 一様乱数による中心極限定理の R スクリプト**
N <- 10000	# 乱数列の長さ
n <- 2	# 標本平均のサイズ
u <- matrix(data=runif(n*N), ncol=n)	# N×n の一様乱数行列
um <- apply(u, 1, mean)	# 行ごとの平均
op <- par(mfrow = c(1, 2))	# 標本平均のヒストグラム
hist(um, breaks=seq(0,1,by=0.02), freq=FALSE, ylim=c(0, 2.5), xlab="", ylab="頻度", main="")
m <- mean(um)	# 標本平均列の平均
s <- sd(um)	# 標本平均列の標準偏差
curve(dnorm(x, m, s), 0, 1, add=TRUE, col="red")	# 正規分布の重ねがき
qqnorm(um, xlab="正規分布分位点", ylab="データ分位点", main="")	# 正規 Q - Q プロット
qqline(um, col="red")	# 正規分布の四分位範囲直線表示
par(op)	#
title(main="一様乱数 2 個の標本平均分布に正規分布 N(0.5, 0.042) を重ね書き")

課題：サンプルサイズ n を大きくして，標本平均の分布が正規分布に近づく様子を確かめよ．

中心極限定理が成り立たないコーシー分布
　コーシー分布は，平均，分散の存在しない分布なので，中心極限定理が成立しない．すなわち，コーシー分布乱数列の標本平均の分布は正規分布に収束しない．
　いま，n = 100 のコーシー分布乱数を発生させ，その標本平均を出す．これを N = 10000回繰り返して平均値の分布を生成する．平均値の平均は -0.17 であり，0 に近いが，平均値の標準偏差は 41 になり，平均値の分布（ヒストグラム）と赤線で表示した正規分布 N(-0.17, 41²) はまったく異なっている．
　これは，コーシー分布は極端な値が出やすいので，平均値も大きくふれるからである．

# コーシー分布標本平均分布の R スクリプト

N <- 10000  			# 乱数列の長さ 　 
n <- 100  				# 標本平均のサイズ 　 
u <- matrix(data=rcauchy(n*N), ncol=n)  # N×n のコーシー分布乱数行列 　 
um <- apply(u, 1, mean)  	# 行ごとの平均 
min(um)					# 平均値の最小値
max(um)					# 平均値の最大値　
uma <- floor(min(um))		
umb <- ceiling(max(um))		 
op <- par(mfrow = c(1, 2))  # 標本平均のヒストグラム 　 
hist(um, breaks=seq(uma,umb,by=0.2), xlim=c(-5,5), freq=FALSE, xlab="", ylab="頻度", main="")  
m <- mean(um); m  			# 標本平均列の平均 　 
s <- sd(um); s  			# 標本平均列の標準偏差 　 
curve(dnorm(x, m, s), -5, 5, add=TRUE, col="red")  # 正規分布の重ねがき 　 
qqnorm(um, xlab="正規分布分位点", ylim=c(-6,6), ylab="データ分位点", main="")  # 正規 Q - Q プロット 　 
qqline(um, col="red")  		# 正規分布の四分位範囲直線表示 　 
par(op)   　 
title(main="コーシー分布乱数の標本平均分布に正規分布を重ね書き")

課題：コーシー分布乱数のメデイアンの分布を正規分布にあてはめたときはどうなるか．

π の近似の誤差分布
　前節で，n = 10000 の一様乱数から生成される点を用いて π の近似値の計算を行った．n を増やしていくと近似の精度がよくなることは明らかであるが，どの程度良くなるかはわからない．これをシミュレーションで考えてみよう．
　まず，n = 1000 の乱数点から π の近似値を出す．乱数点のセットにより π の近似値がふれるので，このセットを多数回（たとえば，N = 10000）繰り返すと「ふれ」の程度が推定できる．n = 1000 からの近似値では「ふれ」の標準偏差は 0.05 程度であり，π の推定値は 3.1 ～ 3.3 に分布し，その近似値の分布は正規分布に非常に近いことがわかる．

**# π の近似の誤差分布の R スクリプト**
N <- 10000	# 点セット列の長さ
n <- 1000	# 乱数点の個数
pi <- NULL	# pi の定義
for(i in 1:N){	# N 回の繰り返し
x <- runif(n, -1, 1)	# (-1, 1) の範囲の一様乱数 n 個生成
y <- runif(n, -1, 1)	#
r <- x^2 + y^2	# 原点からの距離の２乗
pi <- c(pi, 4*length(r[r<1])/n)	# πの近似値ベクトル
}	#
m = mean(pi)	# 近似値の平均
s = sd(pi)	# 近似値の標準偏差
m	# 平均の表示
s	# 標準偏差の表示
op <- par(mfrow = c(1, 2))	#
hist(pi, freq=F, main="")	# 近似値のヒストグラム
curve(dnorm(x, m, s), 2.9, 3.4, add=T, col=2)	# 正規分布の重ね合わせ
qqnorm(pi, xlab="正規分布分位点", ylab="データ分位点", main="")	# 正規 Q - Q プロット
qqline(pi, col="red")	# 正規分布の直線表示
par(op)	#
title(main="一様乱数による π の近似値の誤差分布")	# タイトル

レポート１：乱数点の個数 n を大きくしたとき，誤差の大きさの減少の程度を n のオーダーで表せ．

χ²（カイ２乗）分布（chi-squared distribution)

　正の自由度パラメータ n をもつ確率密度関数が

(4.15)

で表される分布．ガンマ（Γ）分布で，シェープパラメータを a = n/2，スケールパラメータを s = 2，とおいた分布．
平均：E[X ] = as = n/2×2 = n，
分散：Var[X ] = as² = n/2×4 = 2n．

自由度パラメータ n による形状の違い

# 自由度パラメータ n による形状の違いの R スクリプト
curve(dchisq(x, 1), 0, 20) # 自由度 1 の χ² 分布のグラフの表示　
abline(v=0, h=0) # x 軸と y 軸の表示　
curve(dchisq(x, 4), add=T, col=2) # 自由度 4 の χ² 分布のグラフを色 2（赤）で追加
curve(dchisq(x, 10), add=T, col=3) # 自由度 10 の χ² 分布のグラフを色 3（緑）で追加
legend(10, 0.7, c("n = 1", "n = 4", "n = 10"), lty=1, col=c(1, 2, 3))
title(main="χ2 分布の自由度 n による形状の違い") # タイトル　

正規分布する確率変数の２乗和（誤差２乗和）の分布
　Z₁，…，Z_n を互いに独立で同一の標準正規分布 N(0，1) に従う確率変数列とする．このように独立で同一の分布に従う確率変数を iid rv (independent identically distributed random variables)，もしくは無作為標本（random sample）という．このとき，Z_i の n 個の２乗和 U_n は自由度 n の χ² 分布，χ²(n)，に従う．すなわち，

(4.16)

である．標準正規乱数 n = 2 個の２乗和を N = 10000個発生させ，そのヒストグラムをつくり，自由度 n = 2 の χ² 分布を重ね合わせたところよく一致していた（左下図）．さらに，標本累積分布関数に χ² 累積分布関数を重ね合わせた図（右下図）をみると，両者はほとんど一致していることがよくわかる．

# 自由度パラメータ n による χ² 分布の形状の違いの R スクリプト
N <- 10000 # 乱数列の長さ　
n <- 2 # 自由度　
u <- matrix(rnorm(n*N), ncol=n) # N×n の標準正規乱数乱数行列　
u2 <- u^2 # 行列の要素の２乗　
un <- apply(u2, 1, sum) # 行ごとの和　
umx <- ceiling(max(un)) # 最大値を超える整数　
op <- par(mfrow = c(1, 2)) # 　
hist(un, breaks=seq(0,umx,by=0.5), freq=FALSE, xlim=c(0,15), xlab="標準正規乱数（n = 2）の２乗和", ylab="頻度", main="")
curve(dchisq(x, n), 0, 15, add=T, col=2) # 自由度 2 の χ² 分布の重ね合わせ　
plot(ecdf(un), do.points=F, verticals=T, xlim=c(0,12), ylab="累積確率", main="")
curve(pchisq(x,2), 0, 15, add=T, col=2) # 自由度 2 の χ² 累積分布関数の重ね合わせ　
par(op) # 　
title(main="標準正規乱数の２乗和に χ2 分布の重ね合わせ") # タイトル　

課題：２乗和する数 n を大きくした場合も，標準正規乱数の２乗和の分布が自由度 n の χ² 分布に従うことを確かめよ．

F 分布（F distribution)

　正の２つの自由度パラメータ m，n をもつ確率密度関数が

(4.17)

で表される分布．

分子，分母の自由度パラメータ m，n による形状の違い

# 分子，分母の自由度パラメータによる形状の違いの R スクリプト
op <- par(mfrow = c(1, 2)) # 　
curve(df(x, 1, 10), 0, 5, ylim=c(0,1.5), ylab="確率密度", xlab="n = 10") # m = 1，n = 10 の F 分布　
abline(v=0, h=0) # x 軸，y 軸　
curve(df(x, 2, 10), 0, 5, add=T, col=2) # m = 2，n = 10 の F 分布　
curve(df(x, 4, 10), 0, 5, add=T, col=3) # m = 4，n = 10 の F 分布　
curve(df(x, 8, 10), 0, 5, add=T, col=4) # m = 8，n = 10 の F 分布　
legend(2.5, 1.4, c("m = 1", "m = 2", "m = 4", "m = 8"), lty=1, col=1:4) # 凡例　
curve(df(x, 4, 50), 0, 5, col=2, ylab="確率密度", xlab="m = 4") # m = 4，n = 50 の F 分布　
abline(v=0, h=0) # x 軸，y 軸　
curve(df(x, 4, 10), 0, 5, add=T) # m = 4，n = 10 の F 分布　
legend(2.5, 0.7, c("n =10", "n = 50"), lty=1, col=c("black","red")) # 凡例　
par(op) # 　
title(main="F 分布の分子と分母の自由度の違いによる形状") # タイトル　

分散比（variance ratio）の分布としての F 分布
　U と V をそれぞれ独立に自由度 m と n の χ² 分布に従う確率変数とする．このとき，U と V をそれぞれの自由度で割った量の比は，分散比もしくは F 値（F value）と呼ばれ，自由度 m，n の F 分布，F(m，n)，に従う．すなわち，

(4.18)

である．標準正規乱数 m = 4 の２乗和と n = 10 の２乗和をそれぞれ N = 10000個発生させ，自由度 4 と 10 の χ² 分布に従う乱数列を生成する．この χ² 分布乱数を自由度で割った量の比を取った乱数列を生成させたところ，そのヒストグラムは自由度 4，10 の F 分布とよく一致していた．なお，χ² 分布乱数は rchisq(num, df) で発生させることができるが，データが従うと想定されることが多い正規分布に従う確率変数から F 分布が生成されることを実感するため，正規乱数から F 分布ヒストグラムを構成した．

# 独立な χ² 分布乱数の比に F 分布をあてはめる R スクリプト
N <- 10000 # 乱数列の長さ　
m <- 4 # 分子自由度　　
n <- 10 # 分母自由度
um0 <- matrix(rnorm(m*N), ncol=m) # N×m の標準正規乱数乱数行列
um2 <- um0^2 # 行列の要素の２乗　
um <- apply(um2, 1, sum) # 行ごとの和　
um <- um/m # 自由度で割る　
un0 <- matrix(rnorm(n*N), ncol=n) # 　
un2 <- un0^2 # 　
un <- apply(un2, 1, sum) # 　
un <- un/n # 　
fv <- um/un # χ² 分布乱数の比
fmx <- ceiling(max(fv)) # fv の最大値を超える整数
hist(fv, breaks=seq(0,fmx,by=0.2), freq=FALSE, xlim=c(0,6), main="")
curve(df(x, m, n), 0, 6, add=T, col=2) # 自由度 4，10 の F 分布の重ね合わせ
title(main="独立な χ2 分布乱数の比に自由度 4，10 の F 分布のあてはめ", cex.main=0.9)

課題：分子，分母の自由度を変えて，χ² 分布乱数の比が F 分布に従うことを確かめよ．

t 分布（t distribution)

　正の自由度パラメータ n をもつ確率密度関数が

(4.19)

で表される分布．標準正規分布より裾が重く（x が 0 より離れてもなかなか確率密度が 0 に近づかない），自由度が小さいほど裾が重くなる．自由度 n = 1 のときはコーシー分布になり，n = ∞ のときは標準正規分布となる．

自由度パラメータ n による形状の違い

課題： t 分布の確率密度関数は，自由度を n として dt(x, n) で与えられる．n の値を変えることにより，上の図のようなグラフを描け．

自由度パラメータ n による正規分布とのずれ
　t 分布の正規分布との「ずれ」は分布の端の方が顕著となる．このため，自由度による 97.5％点の違いが重要となる．分布の端のずれは，正規 Q - Q プロットでよくわかる．これをみると，自由度 2 の t 分布では，分布の端は正規分布と大きく異なるが，自由度 10 の t 分布ではずれが小さくなっている．

**# 正規分布とのずれの R スクリプト**
N <- 10000; n <- 10	# 乱数列の長さと自由度の設定
y <- rt(N, df=n)	# 自由度 10 の t 分布乱数1000個生成
qqnorm(y, xlab="正規分布分位点", ylab="t 分布分位点", main="")	# 正規 Q - Q プロット
qqline(y, col=2)	# 正規分布の四分位範囲直線表示
title(main="t(10) の正規 Q - Q プロット")	#

課題： 自由度を n の値を変えて，標準正規分布とのずれの様子を正規 Q - Q プロットで確かめよ．

標準正規分布確率変数と χ² 分布確率変数との比
　Z を標準正規分布に従う確率変数とし，U を自由度 n の χ² 分布に従う確率変数で，Z と U は互いに独立であるとする．このとき，Z と U をその自由度 n で割った量の平方根との比は，t 値（t value）と呼ばれ，自由度 n の t 分布，t(n)，に従う．すなわち，

(4.20)

である．なおこの関係より，t² は，自由度 1，n の F 分布，F(1，n)，に従うことがわかる．
　標準正規乱数 Z と，n = 10 の２乗和を発生させてできる自由度 n の χ² 分布乱数をそれぞれ N = 10000個発生させる．この列から標準正規乱数と，自由度 n の χ² 分布乱数をその自由度で割った量の平方根との比である t 値を生成し，そのヒストグラムをみると，自由度 n の t 分布，t(n)，によく一致していた．

**# 標準正規乱数から生成した t 値と t 分布の R スクリプト**
N <- 10000	# シミュレーション回数
n <- 10	# １回のサンプルサイズ（自由度）
un0 <- matrix(rnorm(n*N), ncol=n)	# N×n の標準正規乱数行列
un2 <- un0^2	# 標準正規乱数行列の要素の２乗
un <- apply(un2, 1, sum)	# 要素の２乗の各行の和（自由度 n の χ² 分布乱数 N 個）
unr <- sqrt(un/n)	# 自由度 n の χ² 分布乱数を n で割った平方根
z <- rnorm(N)	# 標準正規乱数 N 個
tv <- z/unr	# t 値 N 個
tmx <- ceiling(max(abs(tv)))	# t 値の絶対値の最大
hist(tv, breaks=seq(-tmx,tmx,by=0.2), freq=FALSE, xlim=c(-5,5), main="")
curve(dt(x, n), add=T, col=2)	# 自由度 10 の t 分布の重ねがき
title(main="標準正規乱数から生成した t 値に t 分布の重ね合わせ")

参考文献

工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004，数理工学社
実践生物統計学－分子から生態まで－（第 1 章，第 2 章），東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社

変数 X	x₁	x₂	x₃	x₄	x₅
確率 P	p₁	p₂	p₃	p₄	p₅

x <- 0:12	#個数
dice <- c(185,1149,3265,5475,6114,5194,3067,1331,403,105,18,0,0)	#回数データ
sum(dice)	#試行回数
pdice <- dice/sum(dice)	#回数の確率
m <- sum(x*pdice)	#平均
p <- m/12	#5，6の出る確率
s2 <- sum(pdice*(x-m)^2)	#分散
v <- 12p(1-p)	#二項分布のもとでの分散
h1 <- dbinom(x, 12, 1/3)	#正しいサイコロのもとでの二項確率分布
h2 <- dbinom(x, 12, p)	#推定確率からの二項確率分布
dicedis <- rbind(pdice,h2,h1)	#行ベクトル－＞行列
colnames(dicedis) <- as.character(0:12)	#列の名前
barplot(dicedis, beside=TRUE, cex.axis=0.8, cex.lab=1.0, xlab="5,6の個数", ylab="確率", legend=c("データ","p=0.338", "p=0.333"))
title(main="Weldon のサイコロ実験の分布")	#グラフタイトル

x <- 0:8	#グラフのx軸の範囲
lam <- 2	#λの定義
yp <- dpois(x,lam)	#ポアソン分布の確率密度
y1 <- dbinom(x, 5, 0.4)	#二項分布（n = 5，p = 0.4）の確率密度
y2 <- dbinom(x, 10, 0.2)	#二項分布（n = 10，p = 0.2）の確率密度
y3 <- dbinom(x, 20, 1/10)	#二項分布（n = 20，p = 0.1）の確率密度
y4 <- dbinom(x, 40, 1/20)	#二項分布（n = 40，p = 0.05）の確率密度
plot(x, y1, type="b", ylab="確率")	#二項分布（n = 5，p = 0.4）のプロット（黒）
points(x, y2, type="b", col="green")	#二項分布（n = 10，p = 0.2）のプロット（赤）
points(x, y3, type="b", col="blue")	#二項分布（n = 20，p = 0.1）のプロット（青）
points(x, y4, type="b", col="purple")	#二項分布（n = 40，p = 0.05）のプロット（紫）
points(x, yp, type="b", col="red")	#ポアソン分布のプロット（赤）
title(main="二項分布（np = 2）がポアソン分布に近づく様子")
# 凡例の記述（locator(1)は，凡例の記述場所をクリックで指定）
legend(locator(1), c("p=0.4", "p=0.2", "p=0.1", "p=0.05", "ポアソン"),
lty=1, col=c("black", "green", "blue", "purple", "red"))

x <- seq(0,1, by=0.01)	# x の定義 0 から 1 まで 0.01 きざみ
y <- dbeta(x, 0.5, 0.5)	# a, b = 0.5, 0.5 のβ分布
plot(x, y, type="l", ylim=c(0,3), col="red")	# y を 0 から 3 に指定（赤）
abline(v=0, h=0)	# y 軸と x 軸の表示
curve(dbeta(x, 1, 1), 0, 1, add=T)	# a, b = 1, 1 のβ分布
curve(dbeta(x, 2, 2), 0, 1, add=T, col="blue")	# a, b = 2, 2 のβ分布（青）
curve(dbeta(x, 1, 3), 0, 1, add=T, col="green")	# a, b = 1, 3 のβ分布（緑）
curve(dbeta(x, 2, 4), 0, 1, add=T, col="purple")	# a, b = 2, 4 のβ分布（紫）
title(main="β分布")	#

curve(dexp(x, 1), 0, 10)	#指数分布の密度関数
abline(v=0, h=0)	# y 軸，x 軸表示
title(main="指数分布")	#

curve(dcauchy(x), -5, 5)	#コーシー分布の密度関数
abline(v=0, h=0)	# y 軸，x 軸表示
title(main="コーシー分布")	#

身長	57	58	59	60	61	62	63	64	65	66
人数	2	4	14	41	83	169	394	669	990	1223
67	68	69	70	71	72	73	74	75	76	77
1329	1230	1063	646	392	202	79	32	16	5	2

x <- 57:77	# 身長（x）の範囲
y <- c(2, 4, 14, 41, 83, 169, 394, 669, 990, 1223, 1329,	# 身長ごとのデータ
1230, 1063, 646, 392, 202, 79, 32, 16, 5, 2)	#
s <- sum(y)	# データ総数
m <- sum(x*y/s); m	# データの平均
v <- sum(y/s*(x-m)^2); v	# データの分散
plot(x, y/s, type="h", xlab="身長（インチ）", ylab="頻度")	# データの棒グラフ表示
curve(dnorm(x, m, sqrt(v)), 57, 77, add=T, col="red")	# 正規密度のグラフ表示
title(main="身長データに対する正規分布のあてはめ")	# タイトル

m <- 67.02; s <- 2.556	# 平均と標準偏差の指定
1 - pnorm(70, mean=m, sd=s)	# 1. 70 までの累積確率を 1 から引く
qnorm(0.9, mean=m, sd=s)	# 2. 累積確率が 0.9 となる身長
pnorm(70, mean=m, sd=s) - pnorm(65, mean=m, sd=s)	# 3. （70 までの累積確率）－（60 までの累積確率）

n <- 5	# 打数
x <- 0:n	# xの範囲
p <- 0.3	# 打率
hit <- dbinom(x, size=n, prob=0.3)	# 二項確率
y <- pbinom(x, size=n, prob=0.3)	# 二項累積確率
m <- n*p	# 平均
sd <- sqrt(np(1-p))	# 標準偏差
op <- par(mfrow = c(1, 2))	#
plot(x, hit, type="h", ylim=c(0,0.4), xlim=c(0,7), xlab="ヒット数", ylab="確率密度")
curve(dnorm(x, mean=m, sd=sd), add=TRUE, col="red")	# 確率密度
plot(x, y, type="s", ylim=c(0,1), xlim=c(0,7), xlab="ヒット数", ylab="累積確率")
curve(pnorm(x, mean=m, sd=sd), add=TRUE, col="red")	# 累積確率
par(op)	#
title(main ="二項分布：n = 5 打数，打率 p = 0.3；正規分布：N(1.5, 1.05) ")

curve(dchisq(x, 1), 0, 20)	# 自由度 1 の χ² 分布のグラフの表示
abline(v=0, h=0)	# x 軸と y 軸の表示
curve(dchisq(x, 4), add=T, col=2)	# 自由度 4 の χ² 分布のグラフを色 2（赤）で追加
curve(dchisq(x, 10), add=T, col=3)	# 自由度 10 の χ² 分布のグラフを色 3（緑）で追加
legend(10, 0.7, c("n = 1", "n = 4", "n = 10"), lty=1, col=c(1, 2, 3))
title(main="χ2 分布の自由度 n による形状の違い")	# タイトル

op <- par(mfrow = c(1, 2))	#
curve(df(x, 1, 10), 0, 5, ylim=c(0,1.5), ylab="確率密度", xlab="n = 10")	# m = 1，n = 10 の F 分布
abline(v=0, h=0)	# x 軸，y 軸
curve(df(x, 2, 10), 0, 5, add=T, col=2)	# m = 2，n = 10 の F 分布
curve(df(x, 4, 10), 0, 5, add=T, col=3)	# m = 4，n = 10 の F 分布
curve(df(x, 8, 10), 0, 5, add=T, col=4)	# m = 8，n = 10 の F 分布
legend(2.5, 1.4, c("m = 1", "m = 2", "m = 4", "m = 8"), lty=1, col=1:4)	# 凡例
curve(df(x, 4, 50), 0, 5, col=2, ylab="確率密度", xlab="m = 4")	# m = 4，n = 50 の F 分布
abline(v=0, h=0)	# x 軸，y 軸
curve(df(x, 4, 10), 0, 5, add=T)	# m = 4，n = 10 の F 分布
legend(2.5, 0.7, c("n =10", "n = 50"), lty=1, col=c("black","red"))	# 凡例
par(op)	#
title(main="F 分布の分子と分母の自由度の違いによる形状")	# タイトル

N <- 10000	# 乱数列の長さ
m <- 4	# 分子自由度
n <- 10	# 分母自由度
um0 <- matrix(rnorm(m*N), ncol=m)	# N×m の標準正規乱数乱数行列
um2 <- um0^2	# 行列の要素の２乗
um <- apply(um2, 1, sum)	# 行ごとの和
um <- um/m	# 自由度で割る
un0 <- matrix(rnorm(n*N), ncol=n)	#
un2 <- un0^2	#
un <- apply(un2, 1, sum)	#
un <- un/n	#
fv <- um/un	# χ² 分布乱数の比
fmx <- ceiling(max(fv))	# fv の最大値を超える整数
hist(fv, breaks=seq(0,fmx,by=0.2), freq=FALSE, xlim=c(0,6), main="")
curve(df(x, m, n), 0, 6, add=T, col=2)	# 自由度 4，10 の F 分布の重ね合わせ
title(main="独立な χ2 分布乱数の比に自由度 4，10 の F 分布のあてはめ", cex.main=0.9)

明治大学新領域創造

統計特論３

東京大学大学院農学生命科学研究科 大森宏

この講義の目的

Rの基本的使い方を学ぶ．

０．方程式の根

ニュートン法

多項式方程式の根

１．確率分布

離散確率分布

離散確率変数

連続型確率分布

連続型確率変数

累積分布関数

２．分布の代表値

平均

分散

積率母関数

３．統計的独立

独立

独立な確率変数の平均と分散

無作為標本

大数の法則

４．１変量分布

4-1．離散分布の例

離散一様分布（discrete uniform distribution)

ベルヌイ分布（Bernoulli distribution)

二項分布（binomial distribution)

幾何分布（geometric distribution)

超幾何分布（hypergeometric distribution)

ポアソン分布（Poisson distribution）

負の二項分布（negative binomial distribution）

離散分布あてはめのまとめ

4-2．連続型分布の例

一様分布（uniform distribution)

β（ベータ）分布（beta distribution)

指数分布（exponential distribution)

Γ（ガンマ）分布（gamma distribution)

コーシー分布（Cauchy distribution)

正規分布（normal distribution)

χ2（カイ２乗）分布（chi-squared distribution)

F 分布（F distribution)

t 分布（t distribution)

参考文献

東京大学大学院農学生命科学研究科　大森宏

χ²（カイ２乗）分布（chi-squared distribution)