カテゴリー: 統計学

前回の記事で名前だけ登場した指数分布についてついでに整理しておきます。
参考: 幾何分布の無記憶性について

指数分布は幾何分布の連続分布版のような確率分布です。
(古さに関係なく一定確率で壊れる機械について)機械が故障するまでの時間や、
(単位時間あたり一定確率で発生する災害について)災害が発生するまでの時間など、
一定確率で発生する何かしらの事象が、次に発生するまでの時間が従う分布です。

数学的には次のように定義されます。
パラメーター$\lambda > 0$に対して、次の確率密度関数を持つ分布を指数分布と呼び、$Exp(\lambda)$と書きます。
$$
f(x;\lambda) = \left\{
\begin{align}
&\lambda e^{-\lambda x} \quad & (x \geq 0)\\
&0 \quad & (x < 0)
\end{align}
\right.
$$
期待値は$\frac{1}{\lambda}$、分散は$\frac{1}{\lambda^2}$です。

モーメント母関数を使うと簡単に導出できますので見ておきましょう。
まず、モーメント母関数は$t<\lambda$の範囲で次のように定義されます。
$$
\begin{align}
M_X(t) &= E(e^{tX})\\
&= \int_{0}^{\infty} e^{tx}\lambda e^{-\lambda x} dx\\
&= \lambda \int_{0}^{\infty} e^{(t-\lambda)x} dx\\
&= \frac{\lambda}{\lambda -t}.
\end{align}
$$

これの微分は簡単ですね。1回微分と2回微分はそれぞれ次のようになります。
$$
\begin{align}
\frac{d}{dt}M_X(t) &= \frac{\lambda}{(\lambda-t)^2}\\
\frac{d^2}{dt^2}M_X(t) &= \frac{2\lambda}{(\lambda-t)^3}.
\end{align}
$$

これを使うと期待値と分散は次のように計算できます。
$$
\begin{align}
E(X) &= \left.\frac{d}{dt}M_X(t)\right|_{t=0}\\
&= \frac{\lambda}{(\lambda-0)^2}\\
&= \frac{1}{\lambda}.
\end{align}
$$
$$
\begin{align}
E(X^2) &= \left.\frac{d^2}{dt^2}M_X(t)\right|_{t=0}\\
&= \frac{2\lambda}{(\lambda-0)^3}\\
&= \frac{2}{\lambda^2}
\end{align}
$$
より、
$$
\begin{align}
V(X) &= E(X^2) – E(X)^2\\
&= \frac{2}{(\lambda-0)^2} – \left(\frac{1}{\lambda}\right)^2\\
&= \frac{1}{\lambda^2}.
\end{align}
$$

前回の記事でも触れました通り、指数分布は無記憶性を持つ連続分布です。
$x_1, x_2 \geq 0$に対して、$P(X\geq x_1+x_2|X\geq x_1) = P(X\geq x_2)$が成り立ちます。
実際、$x>0$とすると、 $P(X\geq x)=e^{-\lambda x}$ ですから、
$$
\begin{align}
P(X\geq x_1+x_2|X\geq x_1) &= \frac{P(X\geq x_1+x_2)}{P(X\geq x_1)}\\
&= \frac{e^{-\lambda (x_1+x_2)}}{e^{-\lambda x_1}}\\
&= e^{-\lambda x_2}\\
&= P(X\geq x_2)
\end{align}
$$
となります。

ここ数回の記事で幾何分布に関連する話を取り上げているので、ついでに幾何分布が持つ無記憶性という性質について紹介します。
これは条件付き確率を用いて、次の数式で表される性質です。

$$
P(X > m+n|X > m) = P(X > ​n) \quad \text{ただし}m, n\geq 0.
$$

まず、幾何分布について上の数式が成り立つとを確認しておきましょう。
$P(X=k) = p(1-p)^{k-1}$ ですから、
$$
\begin{align}
P(X>n) &= \sum_{k=n+1}^{\infty}p(1-p)^{k-1}\\
&= p\cdot\frac{(1-p)^{n}}{1-(1-p)}\\
&= (1-p)^{n}
\end{align}
$$
となります。

よって、
$$
\begin{align}
P(X> m+n|X > m) &= \frac{P(X> m+n \land X > m)}{P(X > m)}\\
&= \frac{P(X > m+n)}{P(X > m)}\\
&= \frac{(1-p)^{m+n}}{(1-p)^{m}}\\
&= (1-p)^{n}\\
&= P(X>n)
\end{align}
$$
となり、幾何分布が冒頭の数式を満たすことが示されました。

これはどういうことか説明します。
幾何分布は確率$p$で成功する独立な試行を、初めて成功するまで繰り返すときに要した回数の分布ですから、
$P(X>n)$というのは、初めて成功するまでに$n+1$回以上かかる確率、言い換えると初めて成功するまでに$n$回以上失敗する確率になります。
これに対して、$P(X > m+n|X > m)$はどういうことかというと、成功するまでに$m+1$回以上かかる、つまりすでに$m$回失敗したという条件のもとで、
成功するのに$m+n+1$回以上かかる、つまり追加で$n$回以上失敗し成功するまでに$n+1$回以上かかる確率を意味します。

この二つが等しいということはどういうことかというと、
成功するまでに$n$回以上失敗する確率は、今の時点で何回失敗しているかという事実に全く影響を受けないということです。

例えば、$1/20$の確率で当たりが出るクジで、連続して20回ハズレを引くと、
そろそろ当たるんじゃないかなという気がしてくる人も多いと思うのですが、
そんなことは全くなく、この先あたりを引くまでにかかる回数の期待値は全く変わってないということを示しています。

この無記憶性は、離散分布の中では幾何分布だけが持つ性質です。
(逆にいうと、離散分布で、無記憶性を持っていたらそれは幾何分布だと言えます。)
このほか、連続分布まで範囲を広げると、指数分布が幾何分布同様に無記憶性を持ちます。

前回の記事で、コンプガチャの期待値と分散を求めましたが、いまいち自信がなかったのでシミュレーションしてみました。
参考: 全種類の景品を集めるのに必要な回数の期待値

おさらいしておくと、$n$種類の景品があるクジを景品が全種類揃うまで引く回数は、
期待値が$n\sum_{k=1}^{n}\frac{1}{k}$, 分散が$n\sum_{k=1}^{n-1}\frac{k}{(n-k)^2}$です。

実際にそのような結果になるのか、仮に$n=20$として、プログラムで繰り返し実行してみましょう。
ちなみに、$n=20$の場合の期待値と分散は次ようになります。

import numpy as np


n = 20
# 期待値
print(sum(n/np.arange(1, n+1)))
# 71.95479314287363

# 分散
print(sum([n*k/(n-k)**2 for k in range(1, n)]))
# 566.5105044223357

シミュレーションに使うために、景品が全種類揃うまでクジを引く関数を実装します。

def complete_gacha(n):
    # 揃ったアイテムの配列
    item_list = []

    # n種類揃うまでクジを引く
    while len(set(item_list)) < n:
        item_list.append(np.random.randint(n))

    return item_list

ためしに$n=5$で実行すると次のような結果が得られます。

print(complete_gacha(5))
# [4, 2, 0, 2, 4, 1, 3]

それでは、この関数を100000回実行し、回数(=返された配列の長さ)のリストを作って、平均値と不変分散を出してみましよう。

result_list = np.array([len(complete_gacha(20)) for _ in range(100000)])

# 期待値
print(result_list.mean())
# 71.86858

# 不偏分散
print(result_list.var(ddof=1))
# 566.9716985005849

試行回数がかなり大きいのもあって、理論値にかなり近い結果が得られましたね。
どうやら前回の記事の結果は一応正しそうです。

いわゆるコンプガチャの問題です。

$n$種類の景品があるクジにおいて、全ての景品を揃えるためには何回程度引けば良いのか(=期待値)を考えていたところ、
うまく解けたので記事にすることにしました。

まず改めて問題の前提条件を整理しておきます。
– クジには$n$種類の景品がある。
– どの景品も当たる確率は等しく$1/n$である。
– クジは無限にあり、過去の景品は将来のクジの確率に影響しない。
– 全集類の景品を最低1回引くまでクジを引き続け、その回数の期待値を求める。

3番目の条件は重要です。要するに景品Aがなかなか出なかったからといって、そのあと景品Aを引ける確率が上がったりしないということです。

当初、場合分けを色々考えてアプローチしていたのですが、次のように考えるとすんなり解けました

この問題は、全ての景品が揃うまでのクジの回数を次のように分解して考えます。

全ての景品が揃うまでの回数 =
1種類目の景品を引くまでの回数
+ 1種類持っている状態から2種類目の景品を引くまでの回数
+ 2種類持っている状態から3種類目の景品を引くまでの回数
・・・
+ $n-1$種類持っている状態から$n$種類目の景品を引くまでの回数

すると、期待値の線形性から次のようになります。
$$
E(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}E(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数})
$$

あとは、$k-1$種類持っている状態から$k$種類目の景品を引くまでの回数の期待値がわかれば良いです。
ここで、$k-1$種類持っているということは、持っていない景品は$n-k+1$種類であり、
これは、$\frac{n-k+1}{n}$の確率で当たり(=まだ持ってない景品)を引けるクジと考えることができます。
そして、そのあたりを引くまでの回数は、$p=\frac{n-k+1}{n}$の幾何分布に従うので、
前回の記事で見た通り、その期待値は$\frac{1}{p}=\frac{n}{n-k+1}$となります。

よって、
$$
E(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}) = \frac{n}{n-k+1}
$$
ですから、
$$
E(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}\frac{n}{n-k+1} = \frac{n}{n} + \frac{n}{n-1} + \cdots + \frac{n}{1}
$$
となり、和の順番を入れ替えて$n$で括ると、
$$
\begin{align}
E(\text{全ての景品が揃うまでの回数}) &= n\left(1+\frac{1}{2}+\frac{1}{3}+\cdots \frac{1}{n}\right)\\
&= n\sum_{k=1}^{n}\frac{1}{k}
\end{align}
$$
となります。
シンプルで美しい結果になりましたね。

ついでにですが、分散も求めておきましょう。
$X$と$Y$が独立の時、期待値同様に分散も$V(X+Y)=V(X)+V(Y)$と分解できることに注意すると以下のようになります。
$$
V(\text{全ての景品が揃うまでの回数}) = \sum_{k=1}^{n}V(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}).
$$
パラメーターが$p$の幾何分布の分散は、$\frac{1-p}{p^2}$ですから、$p=\frac{n-k+1}{n}$を代入すると、
$$
\begin{align}
V(\text{k-1種類持っている状態からk種類目の景品を引くまでの回数}) &= \frac{1-\frac{n-k+1}{n}}{\left(\frac{n-k+1}{n}\right)^2}\\
&=\frac{n(k-1)}{(n-k+1)^2}
\end{align}
$$
となります。
よって、求めたい分散は、
$$
\begin{align}
V(\text{全ての景品が揃うまでの回数}) &= \sum_{k=1}^{n}\frac{n(k-1)}{(n-k+1)^2}\\
&= n\sum_{k=1}^{n-1}\frac{k}{(n-k)^2}
\end{align}
$$
となります。

この次の記事で、幾何分布の性質(期待値)を使いたいのでおさらいしておきます。

おさらい:
成功確率が$p$である独立なベルヌーイ試行を繰り返す時、初めて成功するまでの試行回数$X$が従う確率分布を幾何分布と言います。
(本やサイトによっては、初めて成功する回数ではなく、初めて成功するまでの失敗回数で定義することもあります。
負の二項分布との関係が明確になったり、台が0始まりになったりするので実は個人的にはそちらの方が好みです。)

確率質量関数は次のようになります。
$$
P(X=k) = (1-p)^{k-1}p \quad (k = 1, 2, 3, \cdots).
$$
確率$1-p$で発生する失敗を$k-1$回続けた後に、確率$p$で発生する成功を1回、と考えれば自明ですね。

期待値は$\frac{1}{p}$、分散は$\frac{(1-p)}{p^2}$です。

モーメント母関数は次の式になります。
$$
M_X(t) = \frac{pe^t}{1-(1-p)e^t} \qquad(t<-\log{(1-p)}). $$ 期待値と分散はこのモーメント母関数から算出することもできるのですが、見ての通り、この関数の微分、2回微分を計算していくのは結構手間です。 なので、幾何分布に関しては、期待値と分散は直接計算するのがおすすめです。 (といってもこれもそこそこトリッキーなことをするのですが。) では、期待値から導出していきましょう。まず期待値の定義です。 $$ E(X) = \sum_{k=1}^{\infty}k(1-p)^{k-1}p. $$ ここで、$\frac{1}{1-y}$のマクローリン展開を考えます。 $$ \frac{1}{1-y} = 1+y+y^2+\cdots = \sum_{k=0}^{\infty}y^k. $$ これを両辺$y$で微分すると次の式になります。(しれっと微分と極限の順序交換をしています。数学科の学生さんなどはこの辺りも厳密に議論することをお勧めします。) $$ \frac{1}{(1-y)^2} = \sum_{k=1}^{\infty}ky^{k-1}. $$ この式に、$y=1-p$を代入すると次のようになります。 $$ \sum_{k=1}^{\infty}k(1-p)^{k-1} = \frac{1}{(1-1+p)^2} = \frac{1}{p^2}. $$ 両辺に$p$をかけることで、 $$ E(X) = \frac{1}{p} $$ が導かれました。 成功率が$p=\frac{1}{n}$のベルヌーイ試行は、平均$\frac{1}{p}=n$回で成功する、と考えると直感とよくあいますね。 続いて、分散$V(X)$を導出するために、$E(X^2)$を計算していきましょう。 先ほどのマクローリン展開の微分の式から始めます。 $$ \frac{1}{(1-y)^2} = \sum_{k=1}^{\infty}ky^{k-1}. $$ この両辺に、$y$を掛けます。 $$ \sum_{k=1}^{\infty}ky^{k} = \frac{y}{(1-y)^2} $$ そして、この両辺をもう一回$y$で微分します。 $$ \sum_{k=1}^{\infty}k^2y^{k-1} = \frac{1+y}{(1-y)^3}. $$ $y=1-p$とすると、 $$ \sum_{k=1}^{\infty}k^2(1-p)^{k-1} = \frac{2-p}{p^3}. $$ 両辺に$p$をかけて、 $$ \sum_{k=1}^{\infty}k^2(1-p)^{k-1}p = \frac{2-p}{p^2}. $$ この左辺は$E(X^2)$ですね。 よって、分散$V(X)$は、次のように求まります。 $$ \begin{align} V(X) &=E(X^2) - E(X)^2\\ &=\frac{2-p}{p^2} - \left(\frac{1}{p}\right)^2\\ &=\frac{1-p}{p^2}. \end{align} $$

前回の記事で二項分布の期待値と分散を直接計算したわけですが、記事中でも述べている通り、
二項分布の期待値や分散を導出するのはモーメント母関数を使った方が楽です。
マイナーな方法だけ紹介しているというのも変なので、この記事で二項分布のモーメント母関数について紹介します。

早速ですが、確率関数は
$$P[X=k] = {}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}$$
なので、モーメント母関数は次のようになります。
$$
\begin{align}
M_X(t) &= E(e^{tX})\\
&= \sum_{k=0}^{n} e^{tk}{}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}\\
&= \sum_{k=0}^{n} {}_{n}\mathrm{C}_{k} (pe^{t})^k(1-p)^{n-k}\\
&= (pe^t + 1 -p)^n.
\end{align}
$$
最後の行の式変形は二項定理を使いました。

モーメント母関数を$t$で1回微分すると次式になります。
$$
\frac{d}{dt} M_X(t) = npe^t(pe^t + 1 -p)^{n-1}.
$$
$t=0$を代入することで期待値が得られます。
$$
\begin{align}
E(X) &= \left.\frac{d}{dt} M_X(t)\right|_{t=0}\\
&=np(p+1-p)\\
&=np.
\end{align}
$$
前回の記事の直接計算するのに比べて若干楽なのが感じられると思います。

続いて分散です。モーメント母関数の2回微分は次のようになります。
$$
\frac{d^2}{dt^2} M_X(t) = npe^t(pe^t + 1 -p)^{n-1} + n(n-1)p^2e^{2t}(pe^t + 1 -p)^{n-2}.
$$
正確には$n=1$の場合と$n\geq2$の場合でそれぞれ計算しないといけないのですが、結局どちらの場合も上の式で表されることが証明できます。

ここから二項分布の2次のモーメントが次のように計算できます。
$$
\begin{align}
E(X^2) &= \left.\frac{d^2}{dt^2} M_X(t)\right|_{t=0}\\
&= np(p+1-p)^{n-1} + n(n-1)p^2(p+1-p)^{n-2}\\
&= np +n^2p^2-np^2
\end{align}
$$

よって、二項分布の分散は次のように導出されます。
$$
\begin{align}
V(X) &= E(X^2)-E(X)^2\\
&= np +n^2p^2-np^2 – (np)^2\\
&= np(1-p).
\end{align}
$$

分散に関しては、直接計算するに比べてモーメント母関数を使った方がはるかに楽に導出できましたね。

おさらい:
成功確率が$p$のベルヌーイ試行を独立に$n$回行った時の成功回数を確率変数とする分布を二項分布と呼び、$B(n, p)$と書きます。
確率関数は次の式になります。
$$
P[X=k] = {}_{n}\mathrm{C}_{k}p^k(1-p)^{n-k}.
$$

期待値$E(X)$と分散$V(X)$は次の式で表されることが知られています。
$$
\begin{align}
E(X) &= np.\\
V(X) &= np(1-p).
\end{align}
$$

色々なテキストを見ると期待値や分散の導出はモーメント母関数を使われているのをよく見かけます。
最近復習と計算の練習を兼ねて、これらをモーメント母関数を使わずに定義から直接算出してみたところ、思ったより手こずったので記事に残すことにしました。

式変形の途中で二項係数の次の性質を使いますので注意してみてください。
$x\geq1$の時、$x\cdot{}_{n}\mathrm{C}_{n} = n\cdot {}_{n-1}\mathrm{C}_{x-1}$です。

証明 $x\geq1$とすると、
$$
\begin{align}
x\cdot{}_{n}\mathrm{C}_{x} &= x\frac{n!}{x!(n-x)!}\\
& = n\frac{(n-1)!}{(x-1)!((n-1) – (x-1))!}\\
& = n \cdot {}_{n-1}\mathrm{C}_{x-1}.
\end{align}
$$

それでは、本題に戻って期待値$E(X)$から算出していきます。
$$
\begin{align}
E(X) &= \sum_{x=0}^{n}x\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}x\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x} \qquad& (\text{冒頭の二項係数の性質から})\\
&= np\sum_{x=1}^{n}{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}\\
&= np\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x} \qquad&(\text{x-1をxに置き換え})\\
&= np\{p+(1-p)\}^{n-1}\\
&= np.
\end{align}
$$

以上で、$B(n, p)$の期待値が$np$であることが証明できました。
つぎは分散$V(X)$ですが、$V(X)=E(X^2)-E(X)^2$を利用して算出するので、$E(X^2)$を計算していきます。
$$
\begin{align}
E(X^2) &= \sum_{x=0}^{n}x^2\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}x^2\cdot{}_{n}\mathrm{C}_{x}p^{x}(1-p)^{n-x}\\
&= \sum_{x=1}^{n}xn\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}\\
\end{align}
$$
ここで、$\sum$の中の最初の$x$を、$x=(x-1)+1$と変形して、2項にわけます。
$$
\begin{align}
E(X^2) &= \sum_{x=1}^{n}(x-1)n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}+\sum_{x=1}^{n}n\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x}(1-p)^{n-x}\\
&= np\sum_{x=1}^{n}(x-1)\cdot{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}+np\sum_{x=1}^{n}{}_{n-1}\mathrm{C}_{x-1}p^{x-1}(1-p)^{n-x}\\
&= np\sum_{x=0}^{n-1}x\cdot{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}+np\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}
\end{align}
$$
ここで、$\sum_{x=0}^{n-1}x\cdot{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}$は$B(n-1, p)$の期待値なので、$(n-1)p$です。
さらに、$\sum_{x=0}^{n-1}{}_{n-1}\mathrm{C}_{x}p^{x}(1-p)^{n-1-x}$は$B(n-1, p)$の確率関数の全体の和なので$1$になります。
よって、
$$
E(X^2) = n(n-1)p^2 + np
$$
となります。

あとはこれを使って、
$$
\begin{align}
V(X) &= E(X^2) – E(X)^2\\
&= n^2p^2-np^2+np-(np)^2\\
&= np(1-p)
\end{align}
$$
が導出されました。

確率密度関数が$f(x)$の同一の確率分布に従う$n$個の確率変数$X_1, \dots, X_n$について、これらの最大値が従う分布を考える機会がありました。
初めは少々苦戦したのですが、綺麗に定式化できたので記録として残しておこうと思います。
元々は最大値が従う確率密度関数を直接求めようとしてちまちまと場合分けなど考えていたのですが、
確率密度関数ではなく、累積分布関数を先に求めて、それを微分して確率密度関数を得るようにするとスムーズに算出できました。

最初に記号を導入しておきます。
まず、$X_i$たちが従う確率分布の分布関数を$F(x)$とします。
そして、$Y=\max(X_1, \dots, X_n)$が従う確率分布の確率密度関数を$g(y)$,累積分布関数を$G(y)$とします。

最終的に知りたいのは$g(y)$なのですが、まず$G(y)$の方を算出していきます。
$$
\begin{align}
G(y) &= \text{Yがy以下になる確率}\\
&= X_1, \cdots, X_n \text{が全てy以下になる確率}\\
&= (X_1\text{がy以下になる確率}) \times \cdots \times (X_n\text{がy以下になる確率})\\
&= F(y)^n
\end{align}
$$

こうして、最大値$Y$の累積分布関数が$F(y)^n$であることがわかりました。
確率密度関数は累積分布関数を1回微分することで得られるので次のようになります。
$$
\begin{align}
g(y) &= \frac{d}{dy}G(y)\\
&= \frac{d}{dy}F(y)^n\\
\therefore g(y) &= nF(y)^{n-1}f(y)
\end{align}
$$

ついでに、最小値$Z=\min(X_1, \dots, X_n)$が従う分布の確率密度関数$h(z)$と累積分布関数$H(z)$についても同様に算出できるのでやっておきます。
最大値の場合と同じように$H(z)$の方を求めます。
$$
\begin{align}
H(z) &= \text{Zがz以下になる確率}\\
&= 1-(\text{Zがz以上になる確率})\\
&= 1-(X_1, \cdots, X_n \text{が全てz以上になる確率})\\
&= 1-(X_1\text{がz以上になる確率}) \times \cdots \times (X_n\text{z以上になる確率})\\
&= 1-(1-(X_1\text{がz以下になる確率})) \times \cdots \times (1-(X_n\text{z以下になる確率}))\\
&= 1-(1-F(z))^n
\end{align}
$$
これで、最小値が従う分布の累積分布関数が求まりました。あとはこれを微分して、確率密度関数にします。
$$
\begin{align}
h(z) &= \frac{d}{dz}H(z)\\
&= -n(1-F(z))^{n-1}(-F'(z))\\
\therefore h(z) &= n\{1-F(z)\}^{n-1}f(z)
\end{align}
$$
最大値より若干複雑に見えますが、これで最小値が従う分布も得られました。

先日、データを分析している中でシンプソンのパラドックスが発生しているのを見かけました。
興味深い現象なので、紹介したいと思います。
ただし、業務的な情報は書けないので記事中の用語も設定もデータも全部架空の物です。

２種類のアプリがあったとします。それぞれ旧アプリと、新アプリとします。
そしてそれらのアプリを使っているユーザーがとある属性によってグループA,グループBに分かれていたとします。

ユーザー数の内訳が次のようになっていたとします。(単位:人)

	旧アプリ	新アプリ
グループA	40000	1000
グループB	60000	9000

これらのユーザーのコンバージョン数が次の通りだったとします。

	旧アプリ	新アプリ
グループA	3200	100
グループB	1800	360

コンバージョン率を見ると次のようになりますね。

	旧アプリ	新アプリ
グループA	8%	10%
グループB	3%	4%

どちらのグループのユーザーに対しても、新アプリの方がコンバージョン率が高いことがわかりました。

しかしここで、グループごとに分けて集計することをやめて、新アプリと旧アプリを単純に比較してみます。

	旧アプリ	新アプリ
ユーザー数	100000	10000
コンバージョン数	5000	460
コンバージョン率	5%	4.6%

なんと、新アプリより旧アプリの方がコンバージョン率が高いことになりました。

このように、
集団全体を複数の集団に分けてそれぞれの集団で同じ仮説(今回の例では新アプリの方がコンバージョン率が高い)が成り立っても、
集団全体に対してはそれが成り立たないことがあることをシンプソンのパラドックスと呼びます。

これは$\frac{a}{A}>\frac{b}{B}$ かつ $\frac{c}{C}>\frac{d}{D}$ が成り立ったとしても
$$\frac{a+c}{A+C}>\frac{b+d}{B+D}$$
が成り立つわけじゃないと言う単純な数学的な事実から発生する物です。

今回の例で言えば、新アプリの方がどちらのユーザー群に対しても良い効果をもたらしているので良さそうなのに、
全体の集計だけで旧アプリの方が良いと結論づけてしまうと誤った分析をしてしまうことになります。
注意する必要がありますね。

以前も紹介した、SciPyのカーネル密度推定のメソッド、gaussian_kdeの話です。
参考: SciPyによるカーネル密度推定

最近、多次元(と言っても2次元のデータですが)に対して、カーネル密度推定を行いたいことがあり、
どうせ1次元の場合と同じように使えるのだろうと適当に書いたら思うような動きにならず苦戦しました。

何をやろうとしたかというと、
[[x0, y0], [x1, y1], [x2, y2], …, [xn, yn]]
みたいなデータをそのままgaussian_kdeに渡してしまっていました。

ドキュメント: scipy.stats.gaussian_kde
をよく読むと、
Datapoints to estimate from. In case of univariate data this is a 1-D array, otherwise a 2-D array with shape (# of dims, # of data).
と書いてあります。

僕が渡そうとしたデータは shapeが[データ件数, 2(=次元数)]になっていたのですが、実際は
[2(=次元数), データ件数]の型で渡さないといけなかったわけです。

丁寧なことに、Examplesに取り上げられている例も1Dではなく2Dの例で、np.vstackとか使って書かれているので、以前の記事を書いた時にもっとしっかり読んでおけばよかったです。

使い方がわかったので、2次元のデータに対してやってみます。
サンプルデータは今回はscikit-learnのmake_moonsを使いました。

from sklearn.datasets import make_moons
import matplotlib.pyplot as plt

# データ生成。 今回はラベルは不要なのでデータだけ取得する。
data, _ = make_moons(
    n_samples=200,
    noise=0.1,
    random_state=0
)

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111, title="サンプルデータ")
ax.scatter(data[:, 0], data[:, 1])
plt.show()

サンプルデータはこんな感じです。

生成したデータを、gaussian_kdeにそのまま渡すとうまくいきません。
出来上がるモデルは、200次元のデータ2個を学習したものになっているからです。

from scipy.stats import gaussian_kde
# これはエラーは出ないが誤り
kde = gaussian_kde(data)

# 密度推定した結果を得ようとするとエラーになる。
print(kde.pdf([0, 0]))
# ValueError: points have dimension 1, dataset has dimension 200

きちんと、2次元のデータ200個を学習させるには、データを転置させて渡します。

from scipy.stats import gaussian_kde
# (データを、(次元数, データ件数) の型で渡す)
kde = gaussian_kde(X.T)

kde.evaluate を使って、推定した結果の値を得る時も、(次元数, データ件数)の形で、データを渡す必要があります。
(一点のみなら長さが次元数分の配列を渡せば良いです。)

# 点(1, 1)での値
print(kde.evaluate([1, 1]))
# [0.03921808]
# 4点(-0.5, 0), (0, 1), (0.5, 1), (1, 0) での値
print(kde.evaluate([[-0.5, 0, 0.5, 1], [0, 1, 1, 0]]))
# [0.07139624 0.2690079  0.2134083  0.16500181]

等高線を引いて図示する場合は次のように行います。(公式ドキュメントでは等高線ではなく、imshowで可視化していますね。)
こちらの記事も参考にしてください。
参考: matplotlibで等高線

# 等高線を引く領域のx座標とy座標のリストを用意する
x = np.linspace(-1.5, 2.5, 41)
y = np.linspace(-0.8, 1.3, 22)
# メッシュに変換
xx, yy = np.meshgrid(x, y)
# kdeが受け取れる形に整形
meshdata = np.vstack([xx.ravel(), yy.ravel()])
# 高さのデータ計算
z = kde.evaluate(meshdata)

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111, title="カーネル密度推定")
ax.scatter(data[:, 0], data[:, 1], c="b")
ax.contourf(xx, yy, z.reshape(len(y), len(x)), cmap="Blues", alpha=0.5)
plt.show()

結果がこちらです。

最後に、xもしくはy座標を固定して、断面をみる方法を紹介しておきます。
これはシンプルに、固定したい方は定数値でもう一方のデータと同じ長さの配列を作って、
固定しない方のデータを動かしてプロットするだけです。
x=-1.0, 0.5, 2.0 の3つの直線で切ってみた断面を可視化すると次のようなコードになります。

fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
for x_ in [-1.0, 0.5, 2.0]:
    ax.plot(y, kde.pdf([[x_]*len(y), y]), label=f"x={x_}")
ax.set_ylabel("z")
ax.set_xlabel("y")
ax.legend()
plt.show()

結果がこちら。

だいぶ使い方の感覚が掴めてきました。