information
本記事は2025年4月26日（土）開催のTokyoRにて発表予定の初心者セッション資料その1です。
TokyoR発表資料
その2 「Rの統計関連(dpqr)関数チートシート」
その3 「Rの統計関連(dpqr)関数の書き方」
その4 「誕生日のパラドックスをRの関数で求める方法」
Rの起動画面の謳い文句を改めて実感する出来事
Rを起動するとこんな表示が出てきます。
R version 4.4.3 (2025-02-28 ucrt) -- "Trophy Case"
Copyright (C) 2025 The R Foundation for Statistical Computing
統計に特化したプラットフォームであることを謳っており、それがために統計を学び統計で生計を立てている人はRを利用するのですが、それが紛う事無き真実であると改めて実感した出来事がありました。
統計検定のポアソン分布の問題を解いているときに…
統計検定準１級受けるための勉強しなきゃいかんなあとウダウダ考えているときに下記の問題に出会いました。
問題
ネット通販サイトへの注文が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)
= P(Z > -2.27) = 0.988
でもRならポアソン分布を直接計算できる
おそらく上記が求められている模範解答ですが、せっかくRを使っているのだから、ポアソン分布の計算をしてみようと思いました。確率変数$X$がポアソン分布に従う場合、ある期間に平均$\lambda$回おこる現象が$k$回起こる確率は、
P(X=k) = \frac{\exp(-\lambda)\lambda^k}{k!}
ポアソン分布は離散型確率分布であり、$k$は0以上の整数なので、計算式をコーディングして149件以下の注文確率の総和を求めて、1から引けばいいじゃん、と思いました。計算式にすると、
P(X \geq 150) = 1 - P(X \leq 149) = 1 - \displaystyle\sum_{k = 0}^{149} P(X = k) = 1 - \sum_{k = 0}^{149} \frac{\exp(-\lambda)\lambda^k}{k!}
ということでこんなRコードを作ってみました。
library(conflicted)
library(tidyverse)
lambda <- 180
k <- 149
start_num <- 0
end_num <- k
df <- dplyr::tibble(num = seq(start_num, end_num)) |>
mutate(x_n = str_c("P(X=", num, ")")) |>
mutate(po = (exp(-lambda) * lambda^num) / factorial(num))
> head(df)
# A tibble: 6 × 3
num x_n po
<int> <chr> <dbl>
1 0 P(X=0) 6.71e-79
2 1 P(X=1) 1.21e-76
3 2 P(X=2) 1.09e-74
4 3 P(X=3) 6.53e-73
5 4 P(X=4) 2.94e-71
6 5 P(X=5) 1.06e-69
vec_total <- sum(df$po)
1 - vec_total
しかし結果は
> vec_total
[1] Inf
> 1 - vec_total
[1] -Inf
なぜ???
> tail(df, n = 15)
# A tibble: 15 × 3
num x_n po
<int> <chr> <dbl>
1 135 P(X=135) 0.0000723
2 136 P(X=136) 0.0000957
3 137 P(X=137) Inf
4 138 P(X=138) Inf
5 139 P(X=139) Inf
6 140 P(X=140) Inf
7 141 P(X=141) Inf
8 142 P(X=142) Inf
9 143 P(X=143) Inf
10 144 P(X=144) Inf
11 145 P(X=145) Inf
12 146 P(X=146) Inf
13 147 P(X=147) Inf
14 148 P(X=148) Inf
15 149 P(X=149) Inf
po列の値から見て$X \geq 137$以降の値がinfinityになるとは考えられない。
しかしですねー、冷静に考えればこれは計算中のオーバーフローに違いない。そりゃ$180^{137}$や$137!$は相当な値だからして！実際Rに計算させてみると、
> 180^136
[1] 5.212676e+306
> 180^137
[1] Inf
むしろ、よくぞ$180^{136}$まで計算していただきました。
2025年4月22日追記
私の環境（windows10、R4.5.0)では、下記がdouble型で取り扱える最大値と最小値のようです。
> .Machine$double.xmax
[1] 1.797693e+308
> .Machine$double.xmin
[1] 2.225074e-308
「自然対数をとる」の意味を知る
そこで、どうすれば良いか生成AIに尋ねてみると、自然対数とれば良いとの回答。理系の皆様には当然のことかもしれませんが、先人たちが確立した数学の世界は素晴らしい！ということを文系のわたしくめも50歳となりようやく感じ入りました。
P(X=k) = \frac{\exp(-\lambda)\lambda^k}{k!} \rightarrow \ln\{P(X=k)\} = -\lambda + k \ln \lambda - \ln k!
上記の式でもまだ$k!$があるので、そのときはガンマ関数を使えば良いとの生成AIのご託宣がありました。ガンマ関数とは正の実数$x>0$について、
\Gamma(x) = \int^{\infty}_0 t^{x-1} e^{-t} dt
と定義され、自然数$n$の場合
\Gamma(n) = (n-1)!
が成り立つそうです（ここは理解していないのでそのまま写しています）。
つまり$\ln(k!) = \ln \Gamma(n + 1)$であるので、ガンマ関数の自然対数をとる関数 lgamma(k+1)を使用することで、階乗のオーバーフローも防げるということですね。で、最後に求めた$-\lambda + k \ln(\lambda) - \ln(k!)$の値をexp()関数で元に戻せば良いです。以上をまとめると
\ln P(X = k) = -\lambda + k \ln \lambda - \ln \Gamma(k + 1)
P(X = k) = \exp\{\ln P(X = k)\}
df_improve <- dplyr::tibble(num = seq(start_num, end_num)) |>
mutate(x_n = str_c("P(X=", num, ")")) |>
mutate(ln_po = - lambda + num * log(lambda) - lgamma(num + 1)) |>
mutate(po = exp(ln_po))
vec_total_improve <- sum(df_improve$po)
> vec_total_improve
[1] 0.009910119
> 1 - vec_total_improve
[1] 0.9900899
バッチリ解答が出ました。
♪こんな計算しなくても
しかし、生成AIより更にこんなご託宣も。
> ppois(149, 180)
[1] 0.009910119
> ppois(149, 180, lower.tail = FALSE)
[1] 0.9900899
なんと関数一発じゃないですか！！！
そりゃそーですよね。Rは統計のために作られたものですからこんなことは朝飯前のはずでした…。
今の生成AIは本当に気が利いていて、ついでに149件の注文が入る確率$P(X = 149)$を求める関数までお示ししてくださいました。
> dpois(149, 180)
[1] 0.00191335
ここまできてやっと「あー、prorm()とdnorm()の関係かあ」と気付くのでした。
Rの統計に関する関数が便利なことを改めて思い知らされたので、せっかくだから一覧にしてまとめようというのが趣旨です。が、前置きが長くなりましたので、dpqrチートシートの記事は以下のリンクへ。
information
TokyoRの発表はその2に続きます。
その2 「Rの統計関連(dpqr)関数チートシート」
その3 「Rの統計関連(dpqr)関数の書き方」
その4 「誕生日のパラドックスをRの関数で求める方法」