sec2

変数 X	1	2	3	4	5	6
確率 P	1/6	1/6	1/6	1/6	1/6	1/6

$f(x; \ n, \ p) = {}_n {\rm C}_x p^x (1-p)^{n-x} = {}_n {\rm C}_x p^x q^{n-x} , \ x=0,1,2, \cdots, n$

　 $m(t) ={\rm E}[e^{tX}]=\sum_x e^{tx} {}_n {\rm C}_x p^xq^{n-x} =\sum_x {}_n {\rm C}_x (pe^t )^xq^{n-x}=(pe^t + q)^n$

　 $m'(t)= npe^t(pe^t+q)^{n-1}, \ m''(t)=n(n-1)(pe^t)^2(pe^t+q)^{n-2}+npe^t(pe^t+q)^{n-1}$

幾何分布（geometric distribution)

**二項分布のグラフの 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）

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

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

**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 のサイコロ実験の分布")	#グラフタイトル

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

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

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

超幾何分布（hypergeometric distribution)

**幾何分布のグラフの 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）

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

復元抽出と非復元抽出
　超幾何分布は，標本抽出においての非復元抽出（sampling without replacement）に対する確率分布である．すなわち，m 個の白石と n 個の黒石が入った袋から石を 1 つ取り出し，それを戻さないで 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 は，
$\frac{m}{m+n}=\frac{x}{k}, \ m+n=\frac{km}{x}$
と推定される．これを 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) #推定値の階級分け　 quantile(estimate, c(0.05, 0.1, 0.5, 0.9, 0.95)) # 分位点（パーセンタイル）

生息数推定値分布の 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)  				#推定値の階級分け　 
quantile(estimate, c(0.05, 0.1, 0.5, 0.9, 0.95))	# 分位点（パーセンタイル）

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

ポアソン分布（Poisson distribution）

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

$f(x; \ \lambda)=e^{-\lambda} \ \frac{\lambda^x}{x!}, \ x=0,1,2, \ \cdots$

　
となる分布．積率母関数とその微分は，
　

　 $m(t)={\rm E}[e^{tX}]= \sum^\infty_{x=0} \frac{e^{tx}e^{-\lambda}\lambda^x}{x!} = e^{-\lambda} \sum^\infty_{x=0} \frac{(\lambda e^t)^x}{x!} = e^{-\lambda} e^{\lambda e^t}$

　 $m'(t)=\lambda e^{-\lambda} e^t e^{\lambda e^t}, \ m''(t)=\lambda e^{-\lambda} e^t e^{\lambda e^t} ( \lambda e^t +1)$

　
なので，平均と分散は以下のようになる．ポアソン分布の平均と分散は等しいことに注意．
　

　平均： ${\rm E}[X] = m'(0)=\lambda$ ，

　分散： ${\rm Var}[X]={\rm E}[X^2] - ({\rm E}[X] )^2 = m''(0)-\lambda^2 = \lambda(\lambda+1)-\lambda^2=\lambda$

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

np = λ

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

$f(x)= {}_n {\rm C}_x p^x (1-p)^{n-x} = \frac{n!}{x!(n-x)!} \left( \frac{\lambda}{n} \right)^x \left( 1-\frac{\lambda}{n} \right)^{n-x}$

$= \frac{\lambda^x}{x!} \Bigl[ \bigl(1- \frac{1}{n} \big) \cdot \ \cdots \ \cdot \bigl( 1- \frac{n-1}{n} \bigr) \Bigr] \bigl(1-\frac{\lambda}{n} \bigr)^{-x} \bigl(1-\frac{\lambda}{n} \bigr)^n$

$\to e^{-\lambda} \frac{\lambda^x}{x!}, \ n \to \infty, \qquad \Bigl( \bigl(1-\frac{\lambda}{n} \bigr)^n \to e^{-\lambda} \Bigr)$

となり，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"))

死亡記事件数
　ポアソン分布は，稀な事象の生起モデル（自動車事故，機械の故障，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"))

商品在庫
　ある商店において，ある商品の１日で売れる数は，平均 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 に対し，確率密度が
　

$f(x; \ n, \ p) = {}_{n+x-1} {\rm C}_x p^n (1-p)^x =\frac{\Gamma(x+n)}{\Gamma(n)x!} p^n (1-p)^x, \ x =0,1,2,\cdots,$

$\Gamma(s) =\int_0^{\infty} x^{s-1} e^{-x} dx, \ \Gamma(s+1)=s\Gamma(s), \ \Gamma(n) =(n-1)!$

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

　 $m(t) = {\rm E}[e^{tX}] = \sum_{x=0}^{\infty} e^{tx} {-r \choose x}p^n(-q)^x = \sum_{x=0}^{\infty} {-r \choose x}p^n(-qe^t)^x = \left(\frac{p}{1-qe^t} \right)^n$

　 $m'(t)=p^n(-n)(1-qe^t)^{-n-1}(-qe^t)$

$m''(t)=nqp^n [q(n+1)e^{2t}(1-qe^t)^{-n-2} + e^t (1-qe^t)^{-n-1} ]$

　
となるので，平均と分散は以下のようになる．
　

　平均： ${\rm E}[X] = m'(0) = \frac{nq}{p}$

　分散： ${\rm Var}[X] = {\rm E}[X^2] -({\rm E}[X])^2 = m''(0)-\Bigl( \frac{nq}{p} \Bigr)^2$

$=nqp^n[qp^{-n-2}(n+1)+p^{-n-1}] - \Bigl(\frac{nq}{p} \Bigr)^2 = \frac{nq^2}{p^2} + \frac{nq}{p} = \frac{nq}{p^2}$

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

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

虫歯の数	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)	#分散
nbp <- m/v # 　 nbn <- m*nbp/(1-nbp) plot(x, y/s, type="h", xlab="虫歯数", ylab="確率") yp <- dpois(x, m) points(x, yp, type="b", col="red") ynb <- dnbinom(x, nbn, nbp) points(x, ynb, type="b", col="green") # 負の二項分布表示　 legend(locator(1), c("データ", "ポアソン分布", "負の二項分布"), lty=1, col=c("black","red","green")) title(main="児童虫歯数への負の二項分布のあてはめ")

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

**生物個体分布２シミュレーションの 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²

４．１変量分布

4-1．離散分布の例

離散一様分布（discrete uniform distribution)

ベルヌイ分布（Bernoulli distribution)

二項分布（binomial distribution)

幾何分布（geometric distribution)

超幾何分布（hypergeometric distribution)

ポアソン分布（Poisson distribution）

負の二項分布（negative binomial distribution）

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

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: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）