2008年度生物測定学基礎実験
2008年6月18日

生物個体の分布シミュレーション（続き）

　先週のスクリプトを再掲する．

生物個体分布シミュレーションの R スクリプト（再掲）
n <- 200 # 個体数　
m <- 10 # メッシュ（m2 個）　
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) # カウント分布の平均　
v1 <- sum(table(d)/s*(xp-m1)^2) # カウント分布の分散　
# # 　
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 # カウント分布分散　

カウント数分布（離散確率分布）の平均と分散

離散確率分布（n = 5）
　変数 X （カウント数） 　x1　 　x2　 　x3　 　x4　 　x5　
　確率 P    　p1 　p2 　p3 　p4 　p5

平均：μ = Σi xipi， 分散：σ2 = Σi(xi - μ)2pi

ポアソン分布（Poisson distribution）

　正のパラメータ λ と，0 以上の整数 x に対し，確率密度が
fig16
となる分布．
　平均と分散はともに λ である．

課題４の結果紹介

先週の課題４の，個体が互いに影響を与える例をつくれ，の解答を紹介する．

なわばりモデル
各個体は少なくとも r = 0.05 は離れていなければならばいモデル．これは，半径 r の円盤を敷き詰める ことに相当する．

n <- 200
m <- 10
r <- 0.05
count <- 1
x <- runif(1)
y <- runif(1)
while(count < n)
{
tx <- runif(1)
ty <- runif(1)
for(j in 1:count) {
if( (x[j]-tx)^2 + (y[j]-ty)^2 < r^2 ) break
}
if(j==count) #どの近くにもない
{
x <- c(x,tx)
y <- c(y,ty)
count <- count+1
}
}
　このなわばりモデルではカウント数の分布はポアソン分布に従っていない．ポアソン分布は， 平均と分散が等しいが，このカウント分布は分散が小さいので，ポアソン分布ではないことがわかる． 対応する分布はないように思われる．
ハーレムモデル
全個体数は n = 200 頭，n1 = 20 頭いるオスのまわり（幅 d の一様分布）に複数のメスが住むという想定．
提出されたスクリプトではメスの位置がオスのまわりに上下左右の幅が d = 0.01 と狭かったので， d = 0.05 に拡げた． この例では，各ハーレムはオス１頭とメス９頭の10頭からなっていて，ハーレムの大きさは変動していない．
また，もとのスクリプトでは，(0, 1) 区間からはみ出る個体の処理をしていないので，はみ出た部分は反対側に 表示させるスクリプトを付け加えた．

n1 <- 20
n <- 200
n2 <- n -n1
m <- 10
d <- 0.05 # ハーレムの幅
x1 <- runif(n1)
y1 <- runif(n1)
x2 <- runif(n2 ,x1-d, x1+d)
y2 <- runif(n2 ,y1-d, y1+d)
x <- c(x1, x2)
y <- c(y1, y2)
# (0, 1)区間からはみ出た点の処理
nx <- (1:n)[x>1]
x[nx] <- x[nx] - 1
nx <- (1:n)[x<0]
x[nx] <- x[nx] + 1
ny <- (1:n)[y>1]
y[ny] <- y[ny] - 1
ny <- (1:n)[y<0]
y[ny] <- y[ny] + 1
　このハーレムモデルではカウント数分布は，平均より分散が大きいので，ポアソン分布に従っていない． すなわち，オスがいるメッシュでは個体数が多くなるが， オスがいないメッシュでは個体がほとんどいないからである．

家族モデル
　上のハーレムモデルを少し改変したモデルとして家族モデルを考えることができる．n1 = 40 頭いる親のまわりに 平均 p = 4 のポアソン分布に従う子どもがいるモデル．子どもは親のまわりに標準偏差 d = 0.05 の正規分布 で配置させる．このように設定すると，各家族ごとに子どもの数が異なり，家族分布も親の周りにある程度の 領域に散らばるようになる．

m <- 10
n1 <- 40
p <- 4
n <- n1*(p+1)
d <- 0.05
x0 <- runif(n1); x <- x0
y0 <- runif(n1); y <- y0
for(i in 1:n1){
nf <- rpois(1, p) # 子どもの数を平均 p のポアソン乱数で生成　
for(j in 1:nf){
x <- c(x, rnorm(1, mean=x0[i], sd=d))
y <- c(y, rnorm(1, mean=y0[i], sd=d))
}
}
# (0, 1)区間からはみ出た点の処理
nx <- (1:n)[x>1]
x[nx] <- x[nx] - 1
nx <- (1:n)[x<0]
x[nx] <- x[nx] + 1
ny <- (1:n)[y>1]
y[ny] <- y[ny] - 1
ny <- (1:n)[y<0]
y[ny] <- y[ny] + 1

課題５：なわばりモデル，ハーレムモデル，家族モデルでパラメータの値を変えることにより， 配置が見た目でどのように

負の二項分布（negative binomial distribution）

　成功確率が p で，正のパラメータ n をもち，0 以上の整数 x に対し，確率密度が
fig17
となる分布．Γ関数で確率密度を定義すれば，n は必ずしも整数である必要はない．
　平均は，n(1 - p)/p，分散は，n(1 - p)/p2 である．
　成功確率 p のベルヌイ試行において，n 回の成功が起こるまでの失敗の回数の分布．
　平均に比べて分散が大きいので，分散が大きなカウント分布でのあてはまりがよい場合がある．

課題６：負の二項分布の確率密度は dnbinom(x, n, p) で与えられる．モーメント法を用いてパラメータ n，p の値を推定し，この推定値を用いて，ハーレムモデルと家族モデルのカウント数分布に負の二項分布を あてはめてみよ．
　注）モーメント法：カウントデータの平均や分散が確率分布の平均や分散と等しいとおくことに より，分布パラメータの値を推定すること．ポアソン分布のあてはめでは，カウントデータの平均が ポアソン分布の平均 λ と等しいと考えて，ポアソン分布をあてはめた．

競争モデル
n <- 150
x <- runif(n) ; y <- runif(n)
z <- runif(n) ; w <- runif(n)
t <- NULL
for(i in 1:n){
s <- (z[i]-x[1:n])^2+(w[i]-y[1:n])^2
t <- rbind(t,min(s))
}
r <- 0.0036
num1 <- (1:n)[t > r];num2 <- (1:n)[t < r]
cau <- data.frame(z[num1],w[num1]) #cau=P.caudatum
aur <- data.frame(x,y) #aru=P.aurelia
m <- 10
count <- NULL
for(i in 1:m){
n1 <- (1:n)[cau < (i-1)/m]
n2 <- (1:n)[cau < i/m]
nin <- n2[!n2 %in% n1]
ww <- w[nin]
a <- hist(ww, breaks=0:m/m)
count <- c(count, a$counts)
}
mc <- max(count)
xp <- 0:mc
d <- factor(count, levels=xp)
table(d)
su <- sum(table(d))
m1 <- sum(xp*table(d)/su)
m2 <- mean(count)
v1 <- sum(table(d)/su*(xp-m1)^2)
v2 <- var(count)
count2 <- NULL
for(i in 1:m){
n3 <- (1:n)[aur < (i-1)/m]
n4 <- (1:n)[aur < i/m]
nin2 <- n4[!n4 %in% n3]
yy <- y[nin2]
b <- hist(yy, breaks=0:m/m)
count2 <- c(count2, b$counts)
}
mc2 <- max(count2)
xp2 <- 0:mc2
d2 <- factor(count2, levels=xp2)
table(d2)
su2 <- sum(table(d2))
m3 <- sum(xp2*table(d2)/su2)
m4 <- mean(count2)
v3 <- sum(table(d2)/su2*(xp2-m3)^2)
v4 <- var(count2)
plot(cau,xlim=c(0,1),ylim=c(0,1),xlab="",ylab="",col="blue",cex=2);par(new=T)
plot(aur,xlim=c(0,1),ylim=c(0,1),xlab="",ylab="",col="red",cex=3)
abline(h=0, v=0)
abline(h=1, v=1)
for(i in 1:m) abline(h=i/m, v=i/m, lty=2)
title(main="P.caudatumとP.aureliaの種間競争")
title(sub="BLUE=P.caudatum,RED=P.aurelia")
m1;v1
m3;v3
Copyright (C) 2005, Hiroshi Omori. 最終更新：2085年 6月16日