No.11073 Re: もう一台の件数分布 【青木繁伸】 2009/10/15(Thu) 23:20
なるべく愚直な方法でシミュレーションしてみました。
func <- function(n)n を色々変えて（N はもう少し少なくてもよい）試してみて下さい。
{
na <- nb <- 0
repeat {
if (sample(LETTERS[1:2], 1) == "A") {
na <- na+1
if (na == n) {
return(nb)
}
}
else {
nb <- nb+1
}
}
}
set.seed(777)
n <- 10 # A のアクセス数
N <- 100000 # 試行回数
result <- replicate(N, func(n)) # B へのアクセス数
result <- table(result) # 集計
sufix <- as.integer(names(result))+1
freq <- integer(max(sufix)) # カウント0のセルをスキップしないように
freq[sufix] <- result
x <- 0:(max(sufix)-1)
y <- N*dnbinom(x, n, p=0.5)
barplot(freq, ylim=c(0, max(y*1.05))) # B へのアクセス数
axis(1, x*1.2+0.7, x)
lines(x*1.2+0.7, y, col="red", lwd=2) # 理論度数多角形
図は，クリックすると原寸表示されます
