統計関係が続いているのですが、統計検定に受かるまで続くんじゃないかと思っています笑

基礎は理解している（と思っている）ものの、「幾何分布」がいまいち自分の中で定着しなかったので書きます。

急に「平均は $\frac{1}{p}$ やねん」とだけ書かれているテキストも多いのが原因かなと思い、数学的な導出とpython での実装をミックスした自分専用の「幾何分布大全」でも作ろうかと。
自分のテーマとして「読むだけでもなんとなく理解できる」ことをいつも目的としているため、文章多めかも。

前提

幾何分布は「ある事象（確率p）が初めて観測されるまでの期待値とかその確率はどのくらい？」ということを知りたいことが背景にあります。

統計検定では「営業している人が初めてちゃんと受け答えしてくれる確率を求めよ」みたいなのがありました。

それ以外にも初めてコインで裏が出る期待値とか、色々実生活の思考実験などしても面白そうですね
（ちなみに、幾何分布は離散ですが、連続型を考慮すると指数分布というものになります。ここでは省略）

幾何分布の式

では、実際に具体例をみながらイメージをさらに固めていきます。

例えば、ある一軒家にめちゃくちゃ怖い犬がいて、人が通り過ぎるとpの確率でめちゃくちゃ吠えるとします（だいぶ特異な例を持ち出していることは自覚している）。

ここで、その犬が朝起きてから初めてx人目に向かって初めて吠える確率関数をP(x)とすると、P(x)の式は
$$P(x) = (1 - p)^{x-1}p$$
となります。

ここまではOKかなと思います
（x回目に初めて事象が発生するため、その前x-1回目までは全て(1-p)）

別にこの辺はわかるんだよ！って方、いよいよかもです

幾何分布の期待値と分散

確率変数を求める部分までは覚えなくても大体わかるよっていう人も、急に
期待値 = $\frac{1}{p}$ 、分散 = $\frac{1-p}{p^2}$
というものは暗記しましょう！というのはだいぶ酷な気がしました。

なので、ちゃんと導出しましょう。

期待値

まず期待値の定義式に則り、今回は離散であることも踏まえると

\begin{align}
E[x] &= \sum_{x=1}^{∞}x(1-p)^{x-1}p \\
&=　p\sum_{x=1}^{∞}x(1-p)^{x-1}
\end{align}

実はなぜにテキストが急に期待値とか分散になると下を向いて「そーなりますねん」っていうスタンスになるのかといえば、
マクローリン展開という割とガチガチの数学的手法が必要とされるからだと調べていくうちに知りました。
（統計学は割と文理問わず最近は学ぶ人も多いので、マクローリン展開って言われてもね〜みたいな人への配慮かもしれませんね）

マクローリン展開は以下。
$$\frac{1}{1-x} = 1 + x + x^2 + ・・・ = \sum_{k=0}^{∞}x^k$$

（ここではマクローリン展開がなぜ成り立つのか？みたいな導出の過程は省略します）

これをそのまま使うわけでもなく、両辺xで微分すると、
$$\frac{1}{(1-x)^2} = \sum_{k=0}^{∞} kx^{k-1}$$
なので、x = 1 - p, k=xとすれば、

\begin{align}
E[x] &= p\sum_{x=1}^{∞}x(1-p)^{x-1} \\
&=　p * \frac{1}{(1 - (1-p))^2} \\
&=　p * \frac{1}{p^2} \\
&=　\frac{1}{p}
\end{align}

導出完了。
さすがに、この導出方法を毎度覚えておくことは無謀であり、最終的には「覚える」というのが最短であることは否めませんが、それでも一度ちゃんと導出しておけば、理解も深まると個人的に思います。

分散

次に分散。
定義式は以下。
$$ V[x] = E[x^2] - (E[x])^2$$

先ほどのマクローリン展開を微分したものをさらに微分したものを使います。
微分しちゃいましょう

$$\frac{2}{(1-x)^3} = \sum_{k=2}^{∞} k(k-1) x^{k-2}$$

ここに、先ほどと同じくx = 1-p, k=xを入れると

$$\sum_{x=2}^{∞} x(x-1) (1-p)^{x-2} = \frac{2}{p^3}$$

このままの導出は難しいので、両辺に(1-p)をかけます。
（こういう数学独特の発想みたいなのは今回は追求しません笑）

$$\sum_{x=2}^{∞} x^2 (1-p)^{x-1} = \frac{2(1-p)}{p^3}$$

そこで、以下の天才しか思いつかないであろう期待値を用意します。
$$ E[x(x-1)] = p\sum_{x=1}^{∞} x(x-1)(1-p)^{x-1} = \frac{2(1-p)}{p^2}$$

では、分散を求めます

\begin{align}
V[x] &= E[x^2] - (E[x])^2 \\
&=　E[x(x-1)] + E[x] - (E[x])^2 \\
&=　 \frac{2(1-p)}{p^2} + \frac{1}{p} - \frac{1}{p^2}\\
&=　\frac{1-p}{p^2}
\end{align}

言いたいことはわかります。グッと堪えてください笑

どう考えてもこの導き方は天才以外思いつきませんw
ですが、導出を問題にすることは実生活上意味がないので、あくまで「導出してみた」というくらいで十分だと思います

pythonで幾何分布を扱う

では、最後にpythonで少しみてみます
いつもながら、こんな感じ、くらいで。

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats


# 確率（今回サイコロを想定）
p=1/6
x = np.arange(1, 20)
y = stats.geom.pmf(k=x, p=p)

# それぞれの統計量
x_stats = stats.geom.stats(p=p, loc=0, moments='mv')

# 公式からの導出
mean = 1/p
var = (1-p)/p**2
print(f'公式で計算\n Mean: {mean}, Variance: {var:.0f}')

plt.figure(figsize=(8, 5))
plt.plot(x, y, 'bo')
plt.vlines(x, 0, y) # vlines(x, ymin, ymax)
plt.ylabel('Probability')
plt.title(f'Mean: {int(x_stats[0])}, Variance: {int(x_stats[1])}', fontsize=20)
plt.show()

実際に公式で計算した場合と、pythonで自動で導出した時と一致しました！
そして、幾何分布のグラフイメージもこれでバッチリですな。

今回参考にさせていただいた記事

幾何分布の確率関数からの期待値と分散の導出
Matplotlib で水平線と垂直線をプロットする方法
scipy.stats.geom