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

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

おさらいしておくと、$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} $$

勤務先のテックブログの宣伝

現在、週2回のペースでこの(私用の)ブログを更新していますが、実はここ以外にも勤務先であるオープンワーク社のテックブログにも記事を投稿しています。
投稿頻度はこのブログに比べて非常に低く、僕はまだ4記事しか投稿していないのですが、
その分、1記事1記事は丁寧に時間をかけて書いてきました。

こちらのブログには書いてこなかった、実際に仕事でやっている内容と密接した記事も書いていますので、
もしご興味のあるかたがいらっしゃいましたらこれらの記事も読んでいただけると嬉しいです。

僕が投稿した記事は以下の4記事になります(新しい順。)

OpenWorkの年齢別年収機能の裏側
ABテストの目的と分析時にアナリストが考えていること
企業の”採用力”を指標化しようとして失敗した話
オープンワークのアナリストが分析していること

次の更新は未定ですが、また新しい記事を更新したらこちらのブログでも紹介させていただこうと思います。

二項分布のモーメント母関数とそれを用いた期待値と分散の導出

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

早速ですが、確率関数は
$$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}
$$
が導出されました。

2021年上半期(1月~6月)によく読まれた記事

2021年もあっという間に半分が終わってしまいました。
ここで恒例(?)のよく読まれた記事ランキングを掲載したいと思います。

参考ですが、昨年1年間のよく読まれた記事ランキングはこちらです。
参考: 2020年のまとめ

では早速発表していきます。
集計期間は2021年1月から6月まで。pvでソートしています。

  1. matplotlibのグラフを高解像度で保存する
  2. ネットワークグラフの中心性
  3. pyenvで作成した環境を消す
  4. TensorflowやKerasでJupyterカーネルが落ちるようになってしまった場合の対応
  5. Pythonで連続した日付のリストを作る
  6. numpyのpercentile関数の仕様を確認する
  7. matplotlibのデフォルトのフォントを変更する
  8. INSERT文でWITH句を使う
  9. kerasのto_categoricalを使ってみる
  10. scipyで階層的クラスタリング

相変わらず、プログラミングのちょっとした小ネタのような記事が人気を集めていますね。
このブログ自体そういう記事が多いので、やむを得ないことですが。

データサイエンティストのブログらしい記事としては、
今回ネットワークグラフの中心性の記事が2位にランクインしました。
これを書いたのは昨年なのですが、当時ネットワーク解析について色々勉強して書いた記事だったので、
ニーズがあって嬉しいです。