В случае биномиального распределения при p далеких от 0.5 “качество” доверительного интервала падает… В общем, нужно большее значение n, при симуляции надо бросать монетку не 20 раз для усреднения результата, а хотя бы 100. Иногда неясно, сколько раз “бросать монетку”. Чтобы это выяснить, нужно построить “раститровку” по количеству бросков. Значение среднего будет куда-то стремиться и мы увидим, при каком n оно уже пришло к концу. Вот, например, симуляция расчета среднего при распределении Poisson(λ=0.005,t=?). Доля событий небольшая, всего 0.5%, поэтому нужно наблюдать 200 раз, чтобы увидеть в среднем одно событие. А мы хотим понять, как увеличением t меняется точность предсказания времени, т.е. при каком t в 95% симуляций вычисленный из данных доверительный интервал содержит истинное среднее. Возьмем небольшое количество симуляций (100 раз будем считать, сколько выпало событий за данное время) сначала, чтобы не грузить сильно комп и прикинем необходимое t, чтобы устремиться к некому значению. Видно, что сигнал слегка шумный, но тренд ясен: надо около t=600, чтобы среднее попадало в 95% доверительных интервалов. На большем количестве симуляций качественно картина не меняется, просто сигнал менее шумный.

library(scales)
lambda = 0.005

# first set
nosim = 100
t_r=20:2000

par(mfrow=c(1,2))

means = c()
for (t in t_r){
lhats = rpois(nosim,lambda=lambda*t)/t;
ll = lhats - qnorm(0.975)*sqrt(lhats/t);
ul = lhats + qnorm(0.975)*sqrt(lhats/t);
means = c(means,mean(lhats<ul & lhats > ll))
}

plot(t_r,means,type="l",log="x",col=alpha("black",0.2), main="nosim=100, t=20:2000")
abline(h=1,col="darkgray", lty="dashed")
lines(smooth.spline(t_r,means),col="red")
abline(v=mean(t_r[means==0.9]),col="blue", lty="dotted")
abline(v=mean(t_r[means==0.95]),col="blue", lty="dashed")


# getting on nosim = 100 a value of 500 to get 
nosim = 1000
means = c()
for (t in t_r){
lhats = rpois(nosim,lambda=lambda*t)/t;
ll = lhats - qnorm(0.975)*sqrt(lhats/t);
ul = lhats + qnorm(0.975)*sqrt(lhats/t);
means = c(means,mean(lhats<ul & lhats > ll))
}

plot(t_r,means,type="l",log="x",col=alpha("black",0.2), main="nosim=1000, t=20:2000")
abline(h=1,col="darkgray", lty="dashed")
lines(smooth.spline(t_r,means),col="red")
abline(v=mean(t_r[means==0.9]),col="blue", lty="dotted")
abline(v=mean(t_r[means==0.95]),col="blue", lty="dashed")