適当に疑問に思った事などを計算して確かめてはすぐ捨てるみたいな感じを続けてたんだけどたまには記録に残していこうかなと。
Qiita にしようか迷ったけど何となくこっちに。
浜田 宏「その問題、数理モデルが解決します の第1章で喫煙率の推定が題材になっていて面白い内容だったので考えた事や試した事などを書いてみる。
書籍では喫煙やドラッグなどのように社会的に望ましくない質問への解答はバイアスが掛かるという事で下記のルールによるアンケートを提案している(回答のランダム化)。
- コインを投げる(回答者しか分からない)
- 表が出たら質問へ(正直に)回答
- 裏が出たら質問の内容に依らず「はい」と回答
上記を行う事で回答者単位での質問への回答は保護されるが、全体としての「はい」を回答する比率は算出可能というのはとても上手いと思う。以前に何かの本で読んでた内容だが今回ちゃんと考えてみて改めて上手いなと感じられる。
その上でもう少し踏み込んで考えてみるという事で、得られた回答から喫煙率の範囲を推定する事を試してみたいと思う。
まず下記の設定の下で考えてみる。
\begin{align} N \times p \times q + N \times (1 - p) = K \end{align}
により
\begin{align} q = \frac{1}{p} \left\{ \frac{K}{N} - \left( 1 - p \right) \right\} \end{align} となるので N, p, K にそれぞれ 1000, 0.5, 600 を代入する事で喫煙率 q は 20% と算出される。
前提として質問への回答に嘘がないというのは良いとして、コイン投げによる試行にはブレがあり p = 1/2 の下でも表の出る回数は 500 とは限らずその結果によって推定される q は範囲を持つ事になると考えられる。そこで以下では二項分布に従う変数 H を導入し、それを用いて喫煙率 q の範囲を考えてみる。
\begin{align} H \sim Binom(N, p) \end{align} を満たしつつ
\begin{align} q = \frac{K - (N - H)}{H} \end{align} であり、H を 500 に設定すると N, K がそれぞれ 1000, 600 の下でやはり q は 20% という結果を得る。
以下に実際に H の分布を記載してみると 500 を中心に主に 400〜600 付近に分布している事が分かる。
その上でもう少し踏み込んで考えてみるという事で、得られた回答から喫煙率の範囲を推定する事を試してみたいと思う。
まず下記の設定の下で考えてみる。
- アンケート回答数: N
- コイン投げで表の出る確率: p
- アンケートで「はい(=喫煙者)」と回答された数: K
- 推定される「はい」の比率(=喫煙率): q
\begin{align} N \times p \times q + N \times (1 - p) = K \end{align}
により
\begin{align} q = \frac{1}{p} \left\{ \frac{K}{N} - \left( 1 - p \right) \right\} \end{align} となるので N, p, K にそれぞれ 1000, 0.5, 600 を代入する事で喫煙率 q は 20% と算出される。
前提として質問への回答に嘘がないというのは良いとして、コイン投げによる試行にはブレがあり p = 1/2 の下でも表の出る回数は 500 とは限らずその結果によって推定される q は範囲を持つ事になると考えられる。そこで以下では二項分布に従う変数 H を導入し、それを用いて喫煙率 q の範囲を考えてみる。
- アンケート回答数: N
- コイン投げで表の出る確率: p
- コイン投げで表の出た回数(推定値): H
- アンケートで「はい(=喫煙者)」と回答された数: K
- 推定される「はい」の比率(=喫煙率): q
\begin{align} H \sim Binom(N, p) \end{align} を満たしつつ
\begin{align} q = \frac{K - (N - H)}{H} \end{align} であり、H を 500 に設定すると N, K がそれぞれ 1000, 600 の下でやはり q は 20% という結果を得る。
以下に実際に H の分布を記載してみると 500 を中心に主に 400〜600 付近に分布している事が分かる。
H の 2.5%〜97.5% タイルを算出してみるとその範囲は 469〜531 であり、q が H の単調増加な関数である事からq の範囲(95%)は式(4) の H に 469 と 531 を代入した 14.7%〜24.7% となる。
# H の 2.5%〜97.5% タイル H.lower <- qbinom(p = 0.025, size = 1000, prob = 1/2) H.upper <- qbinom(p = 0.975, size = 1000, prob = 1/2) H <- c(H.lower, H.upper) # q (600 - (1000 - H)) / H # 0.1471215 0.2467043上記が理論上の q の範囲(95%)であり、以下ではシミュレーションによる q の推定を行ってみる。
N <- 1000
K <- 600
# 二項分布に従う乱数の生成
H <- rbinom(n = 1000, size = N, prob = 1/2)
# q の 2.5%〜97.5% タイル
q <- (K - (N - H)) / H
quantile(q, probs = c(0.025, 0.975)) # 0.1471215 0.2467043
# 可視化
tibble::enframe(q) %>%
ggplot(aes(value)) +
geom_histogram(aes(y=..density..), binwidth = 0.005) +
geom_density(fill = "red", alpha = 1/3) +
geom_vline(xintercept = quantile(q, probs = c(0.025, 0.975)), linetype = 2, colour = "blue", alpha = 1/2)
理論値とそんなに離れておらず精度はそんなに悪くない。 さらにアンケート数を増やす(Nを大きくする)方向で q を算出してみると対数の法則に従って q の 95% 区間が狭まりながら 20% に収束していく様子が見て取れる。
上記で見てきたようにコイン投げの回数が二項分布に従ってばらつく事を加味すると喫煙率 q もばらつく事となり、その分布はアンケート数 N を増加させる事で 20% に近づいていく(大数の法則)。
0 件のコメント:
コメントを投稿