コンプガチャのシミュレーション

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

おさらいしておくと、$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}
$$
最大値より若干複雑に見えますが、これで最小値が従う分布も得られました。

シンプソンのパラドックス

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

2種類のアプリがあったとします。それぞれ旧アプリと、新アプリとします。
そしてそれらのアプリを使っているユーザーがとある属性によってグループ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で多次元のカーネル密度推定

以前も紹介した、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()

結果がこちら。

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

Scipyで対数正規分布

対数正規分布の話の続きです。
参考: 前回の記事

Pythonで対数正規分布を扱う場合、(もちろんスクラッチで書いてもいいのですが普通は)Scipyに実装されている、
scipy.stats.lognormを使います。
これに少し癖があり、最初は少してこずりました。

まず、対数正規分布の確率密度関数はパラメーターを二つ持ちます。
次の式の$\mu$と$\sigma$です。
$$
f(x) =\frac{1}{\sqrt{2\pi}\sigma x}\exp\left(-\frac{(\log{x}-\mu)^2}{2\sigma^2}\right).
$$

それに対して、scipy.stats.lognorm の確率密度関数(pdf)は、
s, loc, scale という3つのパラメーターを持ちます。 ややこしいですね。

正規分布(scipy.stats.norm)の場合は、 loc が $\mu$に対応し、scaleが$\sigma$に対応するので、とてもわかりやすいのですが、lognormはそうはなっていなっておらず、わかりにくくなっています。
(とはいえ、ドキュメントには正確に書いてあります。)

Scipyの対数正規分布においては、確率密度関数は
$$
f(x, s) = \frac{1}{sx\sqrt{2\pi}}\exp\left(-\frac{\log^2{x}}{2s^2}\right)
$$
と定義されています。
そして、$y=(x-loc)/scale$と変換したとき、lognorm.pdf(x, s, loc, scale) は、 lognorm.pdf(y, s) / scale と等しくなります。(とドキュメントに書いてあります。)
わかりにくいですね。

$\mu$と$\sigma$とはどのように対応しているのか気になるのですが、上の2式を見比べてわかるとおり、(そしてドキュメントにも書いてある通り、)
まず、引数のsが、パラメーターの$\sigma$に対応します。 (scaleじゃないんだ。)
そして、引数のscaleは$e^{\mu}$になります。(正規分布の流れで、locと$\mu$が関係してると予想していたのでこれも意外です。)
逆にいうと、$\mu=\log{scale}$です。

locはまだ登場していませんが、一旦ここまでの内容をプログラムで確認しておきます。
lognorm(s=0.5, scale=100) とすると、 $\sigma=0.5$で、$\mu=\log{100}\fallingdotseq 4.605$の対数正規分布になります。
乱数をたくさん取って、その対数が、期待値$4.605…$、標準偏差$0.5$の正規分布に従っていそうかどうかみてみます。


from scipy.stats import lognorm
import matplotlib.pyplot as plt
import numpy as np

data = lognorm(s=0.5, scale=100).rvs(size=10000)
log_data = np.log(data)
print("対数の平均:", log_data.mean())
# 対数の平均: 4.607122120944554
print("対数の標準偏差:", log_data.std())
# 対数の標準偏差: 0.49570170767616645

fig = plt.figure(figsize=(12, 6), facecolor="w")
ax = fig.add_subplot(1, 2, 1, title="元のデータ")
ax.hist(data, bins=100, density=True)
ax = fig.add_subplot(1, 2, 2, title="対数")
ax.hist(log_data, bins=100, density=True)
# 期待値のところに縦線引く
ax.vlines(np.log(100), 0, 0.8)

plt.show()

出力された図がこちらです。

想定通り動いてくれていますね。

残る loc ですが、 これは分布の値を左右に平行移動させるものになります。
先出の$y=(x-loc)/scale$ を $x = y*scale + loc$ と変形するとわかりやすいかもしれません、

わかりやすくするために、locに極端な値($\pm 1000$と$0$)を設定して、乱数をとってみました。
(sとscaleは先ほどと同じです。)


locs = [-1000, 0, 1000]

fig = plt.figure(figsize=(18, 6), facecolor="w")
for i, loc in enumerate(locs, start=1):
    data = lognorm(s=0.5, scale=100, loc=loc).rvs(size=1000)
    ax = fig.add_subplot(1, 3, i, title=f"loc={loc}")
    ax.hist(data, bins=100, density=True)

plt.show()

出力がこちらです。

分布の左端が それぞれ、locで指定した、 -1000, 0, 1000 付近になっているのがみて取れると思います。

対数正規分布

最近とあるデータのモデリングのために対数正規分布を使うことがありました。
そのとき使ったScipyのlognormには少し癖があったのでそれを紹介しようと思ったのですが、その前に対数正規分布そのものについて紹介させてもらいます。

対数正規分布とは一言で言ってしまうと、「それに従う確率変数の対数を取ったら正規分布に従うような確率分布」です。
正規分布の対数ではなく、対数を取ったら正規分布なので、注意が必要です。

もう少し正確に書きます。
確率変数$Y$が、正規分布$N(\mu, \sigma^2)$に従うとします。
このとき、$X=e^{Y}$は対数正規分布に従います。

定義に従って確率密度関数$f(x)$を導出しておきましょう。
まず、正規分布 $N(\mu, \sigma^2)$ の確率密度間数$g(y)$は、
$$
g(y) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right)
$$
です。
あとは$y=\log{x}$と$\frac{dy}{dx}=\frac{1}{x}$ に気をつけて、次のように計算できます。

$$
\begin{align}
f(x) &= g(y)\frac{dy}{dx}\\
&=g(\log{x})\frac{1}{x}\\
\therefore f(x) &=\frac{1}{\sqrt{2\pi}\sigma x}\exp\left(-\frac{(\log{x}-\mu)^2}{2\sigma^2}\right).
\end{align}
$$

ここで、$x$の値域は、 $0 <x < \infty$ です。
$\mu$ と $\sigma$ はパラメーターになります。
この導出の過程から明らかですが、それぞれ、対数正規分布に従う確率変数の「対数の」期待値と標準偏差です。
この二つがそのまま対数正規分布の期待値や標準偏差になるわけではないので気をつける必要があります。

対数ではない、確率変数$X$自体の期待値と分散は次の値になります。
参考: 対数正規分布 – Wikipedia
$$
\begin{align}
E[X] &= e^{\mu+\sigma^2/2}\\
V[X] &= e^{2\mu+\sigma^2}(e^{\sigma^2}-1)
\end{align}
$$