Rの統計関連(dpqr)関数の書き方

:::note info
information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その3です。
TokyoR発表資料
[その1 「"The R Foundation for Statistical Computing"を実感する」](https://qiita.com/hige_megane/private/e49ef15222bba6be14d7)
[その2 「Rの統計関連(dpqr)関数チートシート」](https://qiita.com/hige_megane/private/8383b8926cb97f5040cf)
[その4 「誕生日のパラドックスをRの関数で求める方法」](https://qiita.com/hige_megane/private/617064410e1f8008c45d)
:::

#### 問題（上記リンク「その1」の再掲）

ネット通販サイトへの注文が1分間に平均3件入るポアソン過程を仮定する。1時間(60分)で150件以上の注文がある確率を求めよ。


## 正規分布近似による解答
60分間では平均で$3 \times 60 = 180$件の注文が入る。よってポアソン分布$X \sim Po(180)$に従う。ポアソン分布$X \sim Po(180)$は正規分布$Y \sim N(180, 180)$に近似するので連続補正を使い

```math
P(X \geq 150) \approx P(Y > 149.5) = P\left(Z > \dfrac{149.5 - 180}{\sqrt{180}}\right)
```

```math
= 1 - P\left(Z \leq \dfrac{149.5 - 180}{\sqrt{180}}\right)
```

これをRコードで解くとき、チートシートから`pnorm()`関数を使うことはわかったが()内に何を入力すればよいかわからない。そんなときは、`?`を使うと便利です。

```console
> ?pnorm
```

と入力しEnterキーを押下すると、ヘルプ画面が展開されます。Usageを見ると何となく何を入力すべきかイメージができるはずです。

```norm_usage.r
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.r
> pnorm(-1.96)
[1] 0.0249979
```
本問の場合は$mean = \lambda, sd = \sqrt{\lambda}$を()内に指定すればよいです。

```pnorm_answer.r
# 正規分布近似での計算
p_norm <- 1 - pnorm(149.5, mean = 180, sd = sqrt(180))
> p_norm
[1] 0.988497
```

## ポアソン分布として考えたときの解答

ポアソン分布は離散型確率分布であり、$x$は非負の整数なので、余事象の計算として以下のように表現できます。
```math
P(X \geq 150) = 1 - P(X \leq 149) = 1 - \sum_{k = 0}^{149}\dfrac{\exp(-\lambda)\lambda^k}{k!}
```

ここで`ppois()`のUsageを見てみると、

```pois_usage.r
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$だけ指定すれば良いことがわかります。

```ppois.r
p_pois <- 1 - ppois(149, lambda = 180)
> p_pois
[1] 0.9900899
```
:::note warn
`lower.tail =`オプションの注意点

stats公式のArgumentsに以下の記載があります。
> `lower.tail` logical; 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.r
ppois(149, 180, lower.tail = FALSE)
```
と、引数から1を引く必要があります。人間的な思考を勘案つつ間違えないようにするために私は`lower.tail =`オプションは使用せず、1から余事象を引く方法で計算しています。

```ppois.r
1 - ppois(149, 180)
```
慣れてくればとても便利なオプションだと思いますので、うまく使いこなしましょう。
:::

## 単位時間1秒の二項分布として考えたときの解答

ポアソン分布は「非常に短い時間の中で注文が入る確率」を考えますが、「非常に短い時間」を1秒と定義した二項分布として考えると、
1時間=3,600秒なので、$n = 3600$
1時間に平均180件の注文が入るので1秒間に注文が入る確率は$p = \frac{180}{3600}=\frac{1}{20}$より、

```math
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を見ると、

```binom_family.r
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`が必要ですね。
```pbinom.r
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件ですので、以下のようにコーディングできます。

```df_binom.r
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に近づいていくことがわかります。

:::note info
information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その3です。時間があれば資料その4に続きます。
TokyoR発表資料
[その1 「"The R Foundation for Statistical Computing"を実感する」](https://qiita.com/hige_megane/private/e49ef15222bba6be14d7)
[その2 「Rの統計関連(dpqr)関数チートシート」](https://qiita.com/hige_megane/private/8383b8926cb97f5040cf)
[その4 「誕生日のパラドックスをRの関数で求める方法」](https://qiita.com/hige_megane/private/617064410e1f8008c45d)
:::