2015年度生物測定基礎実験

統計解析１

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

この実験の目的

　統計解析ソフトＲを用いて，統計解析の理論と実践を学ぶ

R の基本的使い方

# #以下はコメント文なので，R には読み込まれず，無視される．
#
# 英語得点データを用いて，データの基本統計量の計算演習を行う．
#
#英語の得点	
eigo <- c(	#
36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67,	
50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65,	
65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64,	
63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67,	
47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55)	
length(eigo)	     	 #データ数　
mean(eigo)	     	 #標本平均　
var(eigo)	     	 #標本分散　
sd(eigo)	     	 #標本標準偏差　
boxplot(eigo, main="英語得点の箱ヒゲ図")	     	 #箱ひげ図
boxplot.stats(eigo)	     	 #箱ひげ図用統計量　
summary(eigo)				#英語得点データの要約
hist(eigo, breaks=seq(0, 100, by=5), xlab="English score", ylab="Frequency", main="")
title(main = "英語得点のヒストグラム")	     	 #グラフタイトル
stem(eigo, scale=2)	     	 #幹葉表示
#
# 散布図と回帰直線
#
x <- 1:10  	#連続した自然数　 
y <- c(3,5,6,6,9,10,6,9,10,7)   
plot(x, y)  	# データの散布図　 
reg <- lm(y ~ x)	# 回帰の計算
summary(reg)		# 回帰の結果表示
reg$coefficients	# 回帰係数
abline(reg, col="red")		# 回帰直線の表示
cor(x, y)		# 相関係数

回帰直線：y = 4.2 + 0.527x

主要な確率分布

１．確率分布

離散確率分布

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

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

連続型確率分布

　関数 f(x) が，

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

累積分布関数

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

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

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

２．分布の代表値

平均

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

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

分散

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

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

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

　２つの独立な確率変数 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²

である．

３．二項分布（Binomial distribution）

　独立な n 回のベルヌイ試行（成功確率 p，失敗確率 q = 1 - p ）を行ったときの成功回数 X の分布．
　その確率密度は，

である．X ～ Binom(n，p) と書くこともある． R では二項分布の確率密度は，dbinom(x, n, p) で計算される．
　二項分布の平均と分散はそれぞれ，

${\rm E}[X] = np, \ {\rm Var}[X] = np(1-p)=npq$

　これより，成功確率 p の推定値 X/n の平均と分散は，

${\rm E}[X/n] = {\rm E}[X]/n = np/n = p, \ {\rm Var}[X/n] = {\rm Var}[X]/n^2 = np(1-p)/n^2 = p(1-p)/n$

となる．

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

まず，データから統計量を求める．

平均： $\bar{x}=\frac{1}{26306}(0\cdot185+1\cdot1149+2\cdot3265+\ \cdots \ +12\cdot0)=4.052$ ，

分散： $s^2 = \frac{1}{26306}\{(0-4.052)^2\cdot185+(1-4.052)^2\cdot1149+ \ \cdots \ +(12-4.025)^2\cdot0 \}=2.696$

　1つのサイコロが5か6の目を出す確率を p とすると，12個のサイコロを同時に振って，5 か 6 の目が出る個数は 2 項分布 Binom(12, p) に従うはずである．
2 項分布 Binom(12, p) の平均は 12p，分散は 12p(1－p) であるので， p のデータからの推定値 $\hat{p}$ は，データの平均値を用いて，

$12\hat{p}=4.052, \ \hat{p}=\frac{4.052}{12}=0.338$

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

$12\hat{p}(1-\hat{p})=12\cdot0.338\cdot(1-0.338)=2.685$

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

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

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

Weldon サイコロ実験の解析結果
   モデル       平均       分散

   データ       4.052       2.696

   2項分布（p＝0.338）       4.052       2.685

   2項分布（p＝0.333）       4       2.667

Weldon サイコロ実験の解析結果
モデル	平均	分散
データ	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); m #平均　
p <- m/12 #5，6の出る確率　
s2 <- sum(pdice*(x-m)^2); s2 #分散　
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="Number of 5,6", ylab="Probability", legend=c("data","p=0.338", "p=0.333"))
title(main="Distribution of the experiment of dices by 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); m	#平均
p <- m/12	#5，6の出る確率
s2 <- sum(pdice*(x-m)^2); s2	#分散
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="Number of 5,6", ylab="Probability", legend=c("data","p=0.338", "p=0.333"))
title(main="Distribution of the experiment of dices by Weldon")	#グラフタイトル

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

   女児数       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    6115

女児数	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	6115

問題１：「Weldon のサイコロ実験」と同様な解析を「12人のきょうだい中の女児数のデータ」に対して行え．
具体的には，Weldon サイコロ実験の解析結果の表のようなものをつくり，結果の解釈を行う．

４．超幾何分布（Hypergeometric distribution)

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

で与えられる．R ではこの確率密度は，dhyper(x, m, n, k) で与えられる．

　平均： ${\rm E}[X]=kp$ ，分散： ${\rm Var}[X]=kp(1-p)(m+n-k), \ p=\frac{m}{m+n}$

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

**超幾何分布と二項分布の 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="Hypergeometric and binomial distribution")	#
legend(locator(1),c("Hypergeometric","Binomial"), lty=1, col=c("black", "red"))
# グラフの適当なところをクリックすると，凡例がそこに表示される（locator(1)）．

生息数の推定（捕獲再捕獲法）
　超幾何分布は，生態学で動物の生息数を推定するのに用いられる．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であり，推定値分布の90％は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="Distribution of habitat estimates ")	#タイトル
table(estimate)	#推定値の階級分け
quantile(estimate, c(0.05, 0.1, 0.5, 0.9, 0.95))	# 分位点（パーセンタイル）

問題２：上の例では捕獲（m = 50），再捕獲（k = 50）合わせて100頭を捕獲している．捕獲総数100頭を固定して，捕獲頭数と再捕獲頭数を変えたときの推定精度はどうなるか．また，予算の関係で捕獲総数を半分の50頭にした場合の推定精度はどうなるか．

注）推定精度は，推定値の平均 mv や標準偏差（標準誤差）sd を求め，そこから正規近似により，未知母数 θ の 95 ％信頼区間として，

mv - 1.96×sd < θ < mv + 1.96×sd

を出すのが一般的である．
　しかし，この例では，推定値が無限大に発散してしまうことがあるので推定値の平均や標準誤差を求めることができない．このため，シミュレーションにより行った．

ヒント：推定精度に必要な項目として，推定値の発散割合（table(estimate)），メディアン（中央値），90％もしくは95％信頼区間（quantile(estimate, c(0.05, 0.1, 0.5, 0.9, 0.95)）を比較する．

５．ポアソン分布（Poisson distribution）

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

となる分布．R ではこの確率密度は平均を m として，dpois(x, m) で与えられる．

平均：E[X ] = λ，分散：Var[X ] = λ

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

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

死亡記事件数
　ポアソン分布は，稀な事象の生起モデル（自動車事故，機械の故障，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.8239，分散は0.8294，であった．このデータがポアソン分布に従っていると考える．ポアソン分布の平均は λ なので，λ = 0.8239 のポアソン分布にあてはめてみたところ，非常によく一致していた．
　また，ポアソン分布は，平均と分散が等しいという特徴がある．データの平均と分散の値が近いことから，データはポアソン分布によく適合していることを示している．

**データとモデルとの統計量の比較**
モデル	平均	分散
死亡記事件数	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"))

生物個体の分布シミュレーション
　ある領域内である生物種が一様に分布しているとする．この領域をメッシュで区切ると，メッシュ内で観測される生物の個体数はポアソン分布に従う．このような点の配置は，ポアソン配置（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="Distribution of number of individuals in mesh partitions")	# タイトル
m1	# カウント分布平均
v1	# カウント分布分散

６．負の二項分布（negative binomial distribution）

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

となる分布．Γ関数で確率密度を定義すれば，n は必ずしも整数である必要はない．
　成功確率 p のベルヌイ試行において，n 回の成功が起こるまでの失敗の回数の分布． R ではこの確率密度は，dnbinom(x, n, p) で与えられる．

　平均，分散はそれぞれ，

${\rm E}[X] = \frac{n(1-p)}{p}, \ {\rm Var}[X] = \frac{n(1-p)}{p^2}$

である．

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

**生物個体分布２シミュレーションの 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("data", "poisson", "negativebino"), lty=1, col=c("black","red","green"))
par(op)	# 個体数
title(main="区画内個数へのポアソン分布と負の二項分布のあてはめ")
m1	# 平均
v1	# 分散

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

工場労働者の事故数の分布
　工場労働者一人当たりの遭遇事故数の分布データ（Greenwood et al. 1920）が以下のようであった．

問題３：上のデータをポアソン分布と負の二項分布にあてはめ，どちらの分布のあてはまりがよいか考察せよ．ただし，'5 以上'は 5 とみなすとする．（正確ではない．）

ヒント：ポアソン分布へのあてはめは，上の死亡記事件数の例と同じである．
　負の二項分布へのあてはめは，データの平均と分散が負の二項分布の平均と分散と等しいと考えることにより，負の二項分布のパラメータの推定ができる．これをモーメント法による分布パラメータの推定という．

発展問題：'5 以上'は 5 とみなさない正確な推定法を考えよ．（できなくてよい）

７．一様分布（uniform distribution)

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

で表される分布．

poisson

π の値の推定
　一様乱数を用いて π の値を推定することができる．すなわち，区間（-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)	#単位円内に入った乱数の個数
4*length(pin)/n	# πの推定値．

課題：n を増やして近似の精度を上げてみよう．また，近似の精度をどのように評価すればよいか考えよ．

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

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

で表される分布で，区間 [0, 1] 上の分布である．
　平均，分散はそれぞれ，

${\rm E}[X] = \frac{a}{a b}, \ {\rm Var}[X] = \frac{ab}{(a b 1)(a b)^2}$

である．

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


	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="β分布") #　

**β分布密度関数の 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="β分布")	#

ベータ二項分布

　二項分布に従うと思われていたデータが想定よりも大きな分散をもっていたため，二項分布モデルで説明できないことがある．これを過分散（over disparsion）という．このようなとき，二項分布の成功確率 p がベータ分布に従うというモデルが考えられる．この分布をベータ二項分布という．このように，パラメータに分布を想定したときに生成される分布を伝染分布（contagious distribution）という．これは，パラメータ分布を事前分布とするベイズ的な考え方と近いように見えるが，パラメータの事後分布という概念が無いときはベイズではない。
　ベータ二項分布が出てくるモデルとしては，正しいコインを n 回トスしたときに表が出る回数 x の分布は二項分布 Binom(n, p = 0.5) に従うが，いま，正しくないコインがたくさん袋に入っており，そこからコインを取りだしてトスをしてコインを袋に戻す．これを n 回繰り返したときの表の出る回数 x の分布である．すなわち，袋の中のコインの表が出る確率がコインごとに異なり，それがベータ分布に従っていると想定している．

密度関数

　いま，二項分布 Binom(n, p) の密度関数を f(x |p )，成功確率 p の密度関数を g(p |a, b) とすると，これらはそれぞれ，

$f(x|p) = {}_n{\rm C}_x \ p^x (1-p)^{n-x}$

$g(p|a,b) = \frac{1}{B(a,b)}p^{a-1} (1-p)^{b-1}$

となる．このとき，

$\mu = \frac{a}{a + b}, \ \theta = \frac{1}{a + b}$

とおくと，ベータ分布の平均と分散はそれぞれ，

${\rm E}[p] = \mu, {\rm Var}[p] = \frac{\mu(1-\mu)\theta}{1 + \theta}$

となる．
　ベータ二項分布は x と p の同時分布を p で積分して周辺化した x の周辺分布（marginal distribution）である．すなわち，

$f(x|a, b) = \int_0^1 f(x, p|a, b)dp = \int_0^1 f(x|p)g(p|a,b) dp$
$= \int_0^1 {}_n{\rm C}_x \ p^x (1-p)^{n-x} \frac{1}{B(a,b)}p^{a-1} (1-p)^{b-1} dp$
$= \frac{ {}_n{\rm C}_x}{B(a,b)}\ \int_0^1 p^{x + a -1} (1-p)^{n- x + b - 1}dp$
$= {}_n{\rm C}_x \frac{B(x + a,n-x + b)}{B(a,b)}$

がベータ二項分布の密度関数である．
　ベータ二項分布の平均と分散はそれぞれ，

${\rm E}[x] = n\mu, {\rm Var}[x] = n\mu(1-\mu)\Bigl[1 + \frac{(n-1)\theta}{1 + \theta} \Bigr]$

となることが知られている．ベータ分布で θ → 0 とすると，Var[p] → 0 となり，ベータ分布は μ に集中した分布に退化する．このとき，Var[x] → nμ(1-μ) になり，ベータ二項分布は二項分布 Binom(n, μ) に収束する．

12 人きょうだい中の女児数のデータ（再掲）
1 9世紀末のドイツの病院のデータによると，同じ両親で 12 人きょうだいがいる 6115家族の女児数の数は以下のようであった．問題１では，二項分布モデルにあてはめてみたが，あてはまりが良くなかった．

女児数	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	6115

このデータの平均，分散はそれぞれ，

平均：m = 5.77，分散：v = 3.49

であった．単純な二項分布モデルを考えると女児の生まれる確率の推定値 $\hat{p}$ は，

平均： $m=n \hat{p} \to \hat{p}=m/n$ = 5.77/12 = 0.481

となる．データが二項分布に従っているなら，その平均，分散はそれぞれ，

平均： $m=n \hat{p} = 5.77$ ，分散： $v' = n\hat{p}(1-\hat{p}) = 2.966$

となるはずである．しかし，データの分散 3.49 はこの想定される分散より大きく，過分散になっていると考えられる．

問題４：12 人のきょうだい中の女児数のデータにベータ二項分布をあてはめ，単純な二項分布（p = 0.481）のあてはめと比較せよ．なお，ベータ二項分布のあてはめには，データの平均・分散と分布の平均・分散とが等しいとおくモーメント法を用いよ．
　また，あてはめたベータ分布から女児出生確率の95％信用区間を求めよ．

ヒント１：ベータ二項分布の密度関数のスクリプトは，

y <- lchoose(n, x) + lbeta(x+a, n-x+b) - lbeta(a,b)
h2 <- exp(y)

を用いよ．

ヒント２：パラメータ a，b のベータ分布の95％区間のスクリプトは，

qbeta(c(0.025, 0.975), a, b)

である．

９．正規分布（normal distribution)

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

で表される分布で，N（μ，σ²）と表記する．μ は位置パラメータ（location parameter）で，スケールパラメータ σ を標準偏差（standard deviation）という． R ではこの密度関数は，平均を m，標準偏差を s として，dnorm(x, m, s) で与えられる．なお，dnorm(x) とすると標準正規分布になる．

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

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

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

で表現される．標準正規分布では，-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 ％点

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

**# 標準正規乱数の正規 Q - Q プロットの R スクリプト**
n <- 1000	# 乱数列の長さ
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="Histgram of Normal random numbers")	# タイトル
qqnorm(x, main="")	# 正規 Q - Q プロット
qqline(x, col=2)	# 正規分布の四分位範囲直線表示
title(main="Normal Q - Q plot")	# タイトル
par(op)	#

英語得点データの正規分布へのあてはめ
　英語得点データの標本平均と標本分散を，正規分布の平均と分散とみなして正規分布をあてはめた． Q - Q プロットでみると，下の方の得点分布とあてはまりが悪いようにみえる．

**# 英語得点の正規分布適合度検定の R スクリプト**
eigo <- c( 36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67, # 英語得点データ 50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65, 65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64, 63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67, 47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55) d <- 5 # ヒストグラムの階級幅　 n <- length(eigo) # データ数 m <- mean(eigo) # 平均 s <- sd(eigo) # 標準偏差 op <- par(mfrow = c(1, 2)) # グラフを横に２つ並べて表示 #hist(eigo, breaks=seq(0, 100, by=d), xlab="英語得点", ylab="頻度", main="") # 頻度ヒストグラム #curve(nddnorm(x, m, s), 0, 100, add=TRUE, col="red") # 推定正規分布重ねて表示 hist(eigo, breaks=seq(0, 100, by=d), freq=F, xlab="英語得点", ylab="密度", main="") # 密度ヒストグラム curve(dnorm(x, m, s), 0, 100, add=TRUE, col="red") # 推定正規分布重ねて表示 title(main="英語得点のヒストグラム") # qqnorm(eigo, xlab="正規分布分位点", ylab="英語得点分位点", main="") qqline(eigo, col=2) # Q-Q プロット title(main="正規 Q - Q プロット") # par(op)

# 英語得点の正規分布適合度検定の R スクリプト

eigo <- c( 36,70,56,68,76,60,50,63,62,42,64,60,50,68,71,67,	# 英語得点データ
50,65,67,57,72,64,61,66,46,80,46,51,59,32,55,65, 65,52,57,64,23,57,53,54,38,71,57,69,77,61,51,64,
63,43,65,61,51,69,72,68,53,66,68,58,73,65,62,67, 47,81,47,52,59,33,56,66,67,52,58,65,24,58,54,55)
d <- 5	# ヒストグラムの階級幅 　
n <- length(eigo)	# データ数
m <- mean(eigo)	# 平均
s <- sd(eigo)	# 標準偏差
op <- par(mfrow = c(1, 2))	# グラフを横に２つ並べて表示 
#hist(eigo, breaks=seq(0, 100, by=d), xlab="英語得点", ylab="頻度", main="")	# 頻度ヒストグラム
#curve(n*d*dnorm(x, m, s), 0, 100, add=TRUE, col="red")	# 推定正規分布重ねて表示
hist(eigo, breaks=seq(0, 100, by=d), freq=F, xlab="英語得点", ylab="密度", main="")	  # 密度ヒストグラム
curve(dnorm(x, m, s), 0, 100, add=TRUE, col="red")	# 推定正規分布重ねて表示
title(main="英語得点のヒストグラム")	# 
qqnorm(eigo, xlab="正規分布分位点", ylab="英語得点分位点", main="")
qqline(eigo, col=2)	# Q-Q プロット 
title(main="正規 Q - Q プロット")	# 
par(op)

中心極限定理（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="Normal quantile", ylab="Data quantile", main="")	# 正規 Q - Q プロット
qqline(um, col="red")	# 正規分布の四分位範囲直線表示
par(op)	#
title(main="一様乱数 2 個の標本平均分布に正規分布 N(0.5, 0.042) を重ね書き")

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

ダイズ地上部乾物重データの正規分布へのあてはめ
　長野県中信試験場（現在は，長野県野菜・花卉試験場）で，畦間 75 cm，株間 15 cm で栽培されたダイズ品種エンレイの個体ごとの地上部乾物重データを正規分布にあてはめてみよう．
　データの平均と分散を正規分布のパラメータとおけば，正規分布をあてはめることができる．

データダウンロード

ダイズデータの正規分布あてはめの R スクリプト

dai <- read.csv("daizu.csv"); dai #ダイズ地上部乾物重データ読み込み dim(dai) #データの大きさ表示 head(dai) # データの部分表示 height <- as.vector(as.matrix(dai)) #行列データのベクトル化 boxplot(height) #箱ヒゲ図 boxplot.stats(height) #箱ひげ図用統計量 # m <- mean(height, na.rm=TRUE); m #標本平均 s <- sd(height, na.rm=TRUE); s #標本標準偏差 op <- par(mfrow = c(1, 2)) hist(height, breaks=seq(0,200, by=10), freq=FALSE) curve(dnorm(x, mean=m, sd=s), 0, 200, add=TRUE, col="red") #正規分布 qqnorm(height) #正規 Q - Q プロット qqline(height, col="red") par(op)

ダイズデータの正規分布あてはめの R スクリプト
dai <- read.csv("daizu.csv"); dai #ダイズ地上部乾物重データ読み込み dim(dai) #データの大きさ表示 head(dai) # データの部分表示 height <- as.vector(as.matrix(dai)) #行列データのベクトル化 boxplot(height) #箱ヒゲ図 boxplot.stats(height) #箱ひげ図用統計量 # m <- mean(height, na.rm=TRUE); m #標本平均 s <- sd(height, na.rm=TRUE); s #標本標準偏差 op <- par(mfrow = c(1, 2)) hist(height, breaks=seq(0,200, by=10), freq=FALSE) curve(dnorm(x, mean=m, sd=s), 0, 200, add=TRUE, col="red") #正規分布 qqnorm(height) #正規 Q - Q プロット qqline(height, col="red") par(op)

問題５：夏作物では周辺部分は，隣接個体が少ないことから競合が少なく，また受光や風通しがよくなったりして圃場内部の個体より生育条件が有利になることがある．これを周辺効果という．
　周辺効果の影響を避けるため，圃場実験では周囲の１～２畦や４～５株のデータは解析から外すことが多い．ここで取り上げたデータは，畦１が西側，１行目が北側境界であった．
　周辺部分のデータを除き，地上部乾物重データを正規分布にあてはめ，あてはまりの変化を考察せよ．
　また，異常に生育の悪い個体が一部ある．たとえば，乾物重が20以下の個体を除いてみると正規分布へのあてはまりはどうなるか．

注）データ行列 dai の一部を取り出すときは，例えば

dai2 <- dai[3:25,3:12]
height2 <- as.vector(as.matrix(dai2))

のようにすればよい．また，異常に生育の悪い個体を除くには，

bad <- which(height2 < 20)
height3 <- height2[-bad]

のようにすればよい．

10．χ² 分布（chi-squared distribution)

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

で表される分布．R ではこの密度関数は，自由度を n として，dchisq(x, n) で与えられる．
　χ² 分布は，ガンマ（Γ）分布で，シェープパラメータを a = n/2，スケールパラメータを s = 2，とおいた分布．

平均：E[X ] = as = n/2×2 = n，分散：Var[X ] = as² = n/2×4 = 2n．

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

である．
　標準正規乱数 n = 2 個の２乗和を N = 10000個発生させ，そのヒストグラムをつくり，自由度 n = 2 の χ² 分布を重ね合わせたところよく一致していた

**# 標準正規乱数の２乗和のヒストグラムに χ² 分布をあてはめた 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))	# 最大値を超える整数
hist(un, breaks=seq(0,umx,by=0.5), freq=FALSE, xlim=c(0,15), xlab="Sum of squares of n normal random numbers", main="")
curve(dchisq(x, n), 0, 15, add=T, col=2)	# 自由度 2 の χ² 分布の重ね合わせ
title(main="Histogram of sum of squares of \n normal random numbers")	# タイトル

11．F 分布（F distribution)

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

で表される分布．R ではこの密度関数は，分子自由度を m，分母自由度を n として，df(x, m, n) で与えられる．

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

である．
　標準正規乱数 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="Histgram of ratio of \n independent chi-squared random numbers")

12．t 分布（t distribution)

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

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

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

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

である．なおこの関係より，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="Histgram of t - values ")

参考文献（古い順）

Introduction to the Theory of Statistics, Mood, A. M., Graubill, F. A. & Boes, D. C., 1974, McGRAW-HILL
「実験」生産環境生物学，東京大学大学院農学生命科学研究科生産・環境生物学専攻編，1999，朝倉書店
工学のためのデータサイエンス入門－フリーな統計環境Rを用いたデータ解析－，間瀬茂ら，2004，数理工学社
実践生物統計学－分子から生態まで－（第 1 章，第 2 章），東京大学生物測定学研究室編（大森宏ら）， 2004，朝倉書店
The R Tips データ解析環境 R の基本技・グラフィックス活用集，船尾暢男，2005，九天社
R で学ぶデータマインニング I －データ解析の視点から－，熊谷悦生・船尾暢男，2007，九天社
R で学ぶデータマインニング II －シミュレーションの視点から－，熊谷悦生・船尾暢男，2007，九天社

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); m	#平均
p <- m/12	#5，6の出る確率
s2 <- sum(pdice*(x-m)^2); s2	#分散
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="Number of 5,6", ylab="Probability", legend=c("data","p=0.338", "p=0.333"))
title(main="Distribution of the experiment of dices by Weldon")	#グラフタイトル

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="β分布")	#

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="Histgram of ratio of \n independent chi-squared random numbers")

2015年度生物測定基礎実験

統計解析１

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

この実験の目的

R の基本的使い方

主要な確率分布

１．確率分布

離散確率分布

連続型確率分布

累積分布関数

２．分布の代表値

平均

分散

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

３．二項分布（Binomial distribution）

４．超幾何分布（Hypergeometric distribution)

５．ポアソン分布（Poisson distribution）

６．負の二項分布（negative binomial distribution）

７．一様分布（uniform distribution)

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

ベータ二項分布

密度関数

９．正規分布（normal distribution)

10．χ2 分布（chi-squared distribution)

11．F 分布（F distribution)

12．t 分布（t distribution)

参考文献（古い順）

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

10．χ² 分布（chi-squared distribution)