数学いらずの医科統計学PART2 CHAPTER6で、ポアソン分布に従う事象の話がある。
そこで、それぞれの事象は1回だけ数えられないといけないのだが、飛行機のニアミスがどの程度生じたか調べた研究では、お互いの飛行機の操縦士と副操縦士合わせて4人が各々報告してしまい、1回のニアミスが4回と報告されていたらしい。
NASA blows millions on flawed airline safety survey. 2007
The problem is that NASA appears to have counted some incidents more than once. Pilots were given anonymity, so NASA can't tell when several reports of an incident refer to the same event.
というわけで、単位時間内における本当のニアミス回数をポアソン分布で発生させて、ニアミスがあれば絶対に1人は報告して残りはランダムに報告する状況と、ニアミスがあってもみんなランダムにしか報告しない状況をシミュレーションする。
ポアソン分布への適合度検定をパクって使う。
# 単位時間i番目の最中に本当にニアミスする回数 true_miss <- 1.5 No_unit <- 2000 nearmiss <- rpois(No_unit, true_miss) barplot(table(nearmiss)) # 一回のニアミスでそれぞれのパイロットが報告する確率 pilots_prob <- c(0.5, 0.5, 0.5, 0.5) # 真のニアミス回数 # 一人は報告して、残りはランダム # 残りのランダムな報告 # みんな適当 reports <- matrix(0, No_unit, 4) reports[, 1:2] <- nearmiss for(i in seq(nearmiss)){ if(nearmiss[i] != 0){ # ニアミスがあれば for(j in seq(nearmiss[i])){ # ニアミスした回数だけ愚直にloop回す # 絶対に報告するパイロット who_report <- sample(seq(pilots_prob), size=1) # 他のパイロットが報告する 1 かしない 0 か duplicate_reports <- mapply(function(p) sample(0:1, size=1, prob=c(1-p, p)), pilots_prob) reports[i, 2] <- reports[i, 2] + sum(duplicate_reports[-who_report]) # 一人は報告して、残りはランダム reports[i, 3] <- reports[i, 3] + sum(duplicate_reports[-who_report]) # 残りのランダムな報告 reports[i, 4] <- reports[i, 4] + sum(duplicate_reports) # みんな適当 } } } label0 <- c("真のニアミス回数", "一人は報告して、残りはランダム", "残りパイロットのランダムな報告回数", "みんな適当") xl <- max(reports) yl <- max(sapply(apply(reports, 2, table), max)) par(mfrow=c(1, 4)) for(i in seq(ncol(reports))){ barplot(table(reports[, i]), xlim=c(0, xl), ylim=c(0, yl), main=label0[i]) }
lapply(apply(reports, 2, table), poissondist)
[[1]] ポアソン分布への適合度の検定 data: X[[1L]] X-squared = 4.7681, df = 6, p-value = 0.5739 sample estimates: n lambda 2000.0000 1.4585 [[2]] ポアソン分布への適合度の検定 data: X[[2L]] X-squared = 5032.405, df = 10, p-value < 2.2e-16 sample estimates: n lambda 2000.0000 3.6565 [[3]] ポアソン分布への適合度の検定 data: X[[3L]] X-squared = 1119.898, df = 7, p-value < 2.2e-16 sample estimates: n lambda 2000.0000 2.1975 [[4]] ポアソン分布への適合度の検定 data: X[[4L]] X-squared = 5741.804, df = 10, p-value < 2.2e-16 sample estimates: n lambda 2000.000 3.651
20万単位時間あったとしてシミュレーションしたら、1人は確実に報告する場合と、みんな適当なときと期待値は同じになるようだが、分布の形はけっこう違うようである。
colMeans(reports) [1] 1.508110 3.767765 2.259655 3.768885
