information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その3です。
TokyoR発表資料
その1 「"The R Foundation for Statistical Computing"を実感する」
その2 「Rの統計関連(dpqr)関数チートシート」
その4 「誕生日のパラドックスをRの関数で求める方法」
問題（上記リンク「その1」の再掲）
ネット通販サイトへの注文が1分間に平均3件入るポアソン過程を仮定する。1時間(60分)で150件以上の注文がある確率を求めよ。
正規分布近似による解答
60分間では平均で$3 \times 60 = 180$件の注文が入る。よってポアソン分布$X \sim Po(180)$に従う。ポアソン分布$X \sim Po(180)$は正規分布$Y \sim N(180, 180)$に近似するので連続補正を使い
P(X \geq 150) \approx P(Y > 149.5) = P\left(Z > \dfrac{149.5 - 180}{\sqrt{180}}\right)
= 1 - P\left(Z \leq \dfrac{149.5 - 180}{\sqrt{180}}\right)
これをRコードで解くとき、チートシートからpnorm()関数を使うことはわかったが()内に何を入力すればよいかわからない。そんなときは、?を使うと便利です。
> ?pnorm
と入力しEnterキーを押下すると、ヘルプ画面が展開されます。Usageを見ると何となく何を入力すべきかイメージができるはずです。
dnorm(x, mean = 0, sd = 1, log = FALSE)
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
rnorm(n, mean = 0, sd = 1)
X,q,p,nは引数です。nは生成する乱数の数なのでわかりやすいと思いますが、あとは実際にコードを実行して試しながら引数を探っていくことでRに慣れていくことができます。logやlog.pは表現のとおり自然対数を取るかどうかのオプションです。lower.tailについては後述します。
mean = 0, sd = 1の記載のとおり、normファミリーは平均0、標準偏差1の標準正規分布を基本としています。ですから、()内は引数だけでも下記のとおり標準正規分布における値を返してくれます。
> pnorm(-1.96)
[1] 0.0249979
本問の場合は$mean = \lambda, sd = \sqrt{\lambda}$を()内に指定すればよいです。
# 正規分布近似での計算
p_norm <- 1 - pnorm(149.5, mean = 180, sd = sqrt(180))
> p_norm
[1] 0.988497
ポアソン分布として考えたときの解答
ポアソン分布は離散型確率分布であり、$x$は非負の整数なので、余事象の計算として以下のように表現できます。
P(X \geq 150) = 1 - P(X \leq 149) = 1 - \sum_{k = 0}^{149}\dfrac{\exp(-\lambda)\lambda^k}{k!}
ここでppois()のUsageを見てみると、
dpois(x, lambda, log = FALSE)
ppois(q, lambda, lower.tail = TRUE, log.p = FALSE)
qpois(p, lambda, lower.tail = TRUE, log.p = FALSE)
rpois(n, lambda)
ポアソン分布は平均も分散も$\lambda$なので、引数と$\lambda$だけ指定すれば良いことがわかります。
p_pois <- 1 - ppois(149, lambda = 180)
> p_pois
[1] 0.9900899
lower.tail =オプションの注意点
stats公式のArgumentsに以下の記載があります。
lower.taillogical; if TRUE (default), probabilities are $P[X \leq x]$, otherwise, $P[X > x]$.
lower.tail = FALSEを使うときは$P[X > x]$が適用されます。
連続型確率分布のときにはあまり問題にならないと思いますが、離散型確率分布でlower.tail = FALSEを適用するときは$P[X \geq x]$でない（$P(X = x)$を含まない）点に注意が必要です。$X \sim Po (180)$のポアソン分布について$P(X \geq 150)$の確率を求めるときは
ppois(149, 180, lower.tail = FALSE)
と、引数から1を引く必要があります。人間的な思考を勘案つつ間違えないようにするために私はlower.tail =オプションは使用せず、1から余事象を引く方法で計算しています。
1 - ppois(149, 180)
慣れてくればとても便利なオプションだと思いますので、うまく使いこなしましょう。
単位時間1秒の二項分布として考えたときの解答
ポアソン分布は「非常に短い時間の中で注文が入る確率」を考えますが、「非常に短い時間」を1秒と定義した二項分布として考えると、
1時間=3,600秒なので、$n = 3600$
1時間に平均180件の注文が入るので1秒間に注文が入る確率は$p = \frac{180}{3600}=\frac{1}{20}$より、
P(X \geq 150) = 1 - P(X \leq 149) = 1 - \sum_{k = 0}^{149}
\binom{3600}{k}\left(\dfrac{1}{20}\right)^{k} \left(\dfrac{19}{20}\right)^{3600-k}
例えば$\displaystyle\binom{3600}{149}$を計算するだけでも私には無理ですが、これも関数一発で計算してくれるのはうれしいことです。
ここでpbinom()のUsageを見ると、
dbinom(x, size, prob, log = FALSE)
pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)
qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE)
rbinom(n, size, prob)
引数以外はsizeとprobが必要ですね。
p_binom_second <- 1 - pbinom(149, size = 3600, prob = 1/20)
> p_binom_second
[1] 0.9915267
（おまけ）二項分布解答の拡張版
二項分布と考えたときは「非常に短い時間」の設定は任意です。上記では1秒単位で取りましたが、もっと短い時間設定にしたほうが精度は上がります。1秒単位だとsize = 3600ですが、これを2倍、3倍にして単位時間を0.5秒ごと、0.33秒ごと…に設定するとどうなるでしょうか？1時間の注文数は単位時間に関係なく180件ですので、以下のようにコーディングできます。
library(tidyverse)
df_binom <- dplyr::tibble(n = 1:20) |>
mutate(n_size = 3600 * n) |>
mutate(p_binom = 1 - pbinom(149, size = n_size, prob = 180 / n_size))
> df_binom
> df_binom$p_binom
# A tibble: 20 × 3
n n_size p_binom
<int> <dbl> <dbl>
1 1 3600 0.992
2 2 7200 0.991
3 3 10800 0.991
4 4 14400 0.990
5 5 18000 0.990
6 6 21600 0.990
7 7 25200 0.990
8 8 28800 0.990
9 9 32400 0.990
10 10 36000 0.990
11 11 39600 0.990
12 12 43200 0.990
13 13 46800 0.990
14 14 50400 0.990
15 15 54000 0.990
16 16 57600 0.990
17 17 61200 0.990
18 18 64800 0.990
19 19 68400 0.990
20 20 72000 0.990
[1] 0.9915267 0.9908200 0.9905791 0.9904578 0.9903847 0.9903358 0.9903008 0.9902745 0.9902541 0.9902377 0.9902243
[12] 0.9902131 0.9902037 0.9901956 0.9901885 0.9901824 0.9901770 0.9901721 0.9901678 0.9901639
単位時間を短くすればするほどポアソン分布で求めた値0.9900899に近づいていくことがわかります。
information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その3です。時間があれば資料その4に続きます。
TokyoR発表資料
その1 「"The R Foundation for Statistical Computing"を実感する」
その2 「Rの統計関連(dpqr)関数チートシート」
その4 「誕生日のパラドックスをRの関数で求める方法」