2023年12月7日Python,データサイエンス,統計

こちらの記事で正規分布の和、差、積、商の分布を可視化しましたが、それを通してコーシー分布はやっぱりおかしいなと思ったので、コーシー分布の振る舞いをもっと詳しくみていこうと思います。

コーシー分布とは？

連続な確立分布で、基本的な情報をまとめると以下のようになります。

なお、\(x_0 \in \mathbb{R}, \gamma > 0\)はそれぞれ位置母数、尺度母数とよばれる、分布の形や位置を決めるパタメータです。また、コーシー分布は自由度1のt分布です。

これらの導出や積率母関数、特製関数の計算は数学の景色さんが詳しいのでご参照ください。

また、各種証明に加え、ある確率変数に対してどんな変換をすればコーシー分布になるかなど、@Mark-Nさんの「初めてのコーシー分布」が非常に参考になります。

この記事ではコーシー分布に従う確率変数の特徴や振る舞いの異常さをいろんな可視化を通して確認し、理解しようとする（コーシー分布と仲良くなる）ことを試みます。

なお、本記事では\(x_0 = 0, \gamma = 1\)の標準コーシー分布についてみていき、健常な分布の代表として標準正規分布、自由度5のt分布も一緒に可視化していきます。

コーシー分布のここが異常

概形でまともな分布だと思わせてくるが、ヒストグラムを描いた時点で他と様子が異なる猫被り分布

標準コーシー分布、自由度5のt分布、標準正規分布の密度関数を重ねて図示してみます。なお、以降では基本的にそれぞれコーシー分布、t分布、正規分布と書きます。

コーシー分布、t分布、標準正規分布の順に裾が重たくなっています。t分布の自由度を無限に大きくすると標準正規分布になり、1にするとコーシー分布になるということからも、この大小関係は直感的にわかりやすいです。

コーシー分布は左右対称、釣鐘型であり、平均、期待値、中央値が一致しそうで、数学的にかなり扱いやすい分布に見えます。

次に、実際にそれぞれの分布に従う乱数を生成してヒストグラムを密度で描いてみます。

乱数を発生させて相対頻度にしてヒストグラムを表示し、その上にそれぞれの真の密度関数を重ねています。t分布、正規分布はいい感じですが、コーシー分布は早速おかしいですね。x軸を狭めて拡大すると次の通りです。

このように、ばらつきが大きすぎてx軸の範囲やbinの数をうまく設定しないと姿さえ現してくれません。

期待値が発散する→大数の法則やデルタ法が使えない

\(X_1, X_2, i.i.d. \sim (\mu, \sigma)\)（互いに独立に同じ分布に従う）を見ると「Nが大きければ中心極限定理で平均はは正規分布に従う！」とやりたくなりますが、コーシー分布の場合期待値が存在しないため、正規分布に従うと言わせてくれません。

期待値が発散するのは積分で期待値計算をしてみればわかりますが、分布の形が正規分布に似てるのもあって、どうも納得いきません。ということで、30万サンプル生成し、1サンプルづつ増やして計算した平均値の推移を可視化してみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import cauchy, t, norm

def moving_average(x):
    ma = []
    for i in range(len(x)):
        ma.append(np.mean(x[:(i+1)]))
    return ma

# それぞれの分布からサンプルを30万点生成
np.random.seed(0)
x_cauchy = cauchy.rvs(size=300_000)
x_t = t.rvs(size=300_000, df=5)
x_norm = norm.rvs(size=300_000)

# 1サンプルづつ増やして計算した平均をプロット
fig, axes = plt.subplots(3, 1, tight_layout=True, figsize=(8,8))
axes[0].plot(moving_average(x_cauchy))
axes[1].plot(moving_average(x_t))
axes[2].plot(moving_average(x_norm))
axes[0].set_ylabel("コーシー分布")
axes[1].set_ylabel("t分布")
axes[2].set_ylabel("正規分布")
plt.show()

t分布、正規分布が1万サンプルくらいで早々に収束しているところ、コーシー分布はいつまでもフラフラしています。

もっと極端な場合を見てみます。サンプルサイズ別に平均値を計算することを10回繰り返し、平均値の平均値を棒グラフとして描いてその上に平均値の標準誤差を重ねてみます。ちょっとややこしいですが、要するにサンプルサイズによって標本平均がどの程度ばらつくのかを確認します。

t分布、正規分布はこれくらい大きなサンプルサイズだと0に収束していますね。一方コーシー分布はというと、10万のサンプルを集めても平均値が大きくばらつきます。ただ、一応平均値はサンプルサイズが増えるにつれて小さくなっているように見えます、もっとサンプルサイズを増やせば収束するのではないでしょうか？コーシー分布だけもっと極端な場合も見てみます。

収束の希望が打ち砕かれました。1億サンプル集めても収束しない標本平均ほど厄介なものも珍しいですね。これ以上極端な状況は実行時間の都合で諦めました。

分散も発散する

こちらも期待値と同様、積分計算により確かめられますが、やはりあの分布形状を見ると納得いきません。期待値の場合と同様に、サンプルサイズごとの標本分散の推移とより極端なサンプルサイズでの分散の平均と標準誤差を見てみます。

やはりコーシー分布だけ全く収束せず、全体的に見ると徐々に増加しているように見えますね。時折発生する莫大な外れ値により平均、分散ともにこのようなカクカクしたグラフになっているようです。

次に極端なサンプルサイズの場合を見てみます。

案の定、t分布と正規分布が安定している（コーシーが大きすぎて潰れてますが、ほぼ真値です）のに対し、コーシーはどんどん分散が大きくなっています。より極端な場合は以下の通りです。

サンプルサイズに比例して標本分散が大きくなっています。観測すればするほど未知の値が出てくる、恐ろしい分布ですね。

\(k\)次モーメントも発散

2,3,4次モーメントに関しても標本平均、標本分散と同様に、サンプルサイズに応じてどう変化していくのかをみてみます。他の分布はどうせ収束するので今回はコーシー分布だけ可視化します。他と同様、\(k\)次モーメントは\(k\)乗の平均値で推定します。

やはりダメですね。省略しますが、t分布と正規分布に関しては収束しそうな動きをしていました。

外れ値の外れ具合が半端じゃない

これまでの結果で分かる通りですが、コーシー分布に従う確率変数は外れ値の程度が尋常じゃないです。本当にあの密度関数からあれほど異常な外れ値が出るのでしょうか。3つの分布に対して、極端な値を取る確率（\(P(|X| > x)\)）を計算してみます。

thresh = 100
out_cauchy = 1 - cauchy.cdf(thresh) + cauchy.cdf(-thresh)
out_t = 1 - t.cdf(thresh, df=5) + t.cdf(-thresh, df=5)
out_norm = 1 - norm.cdf(thresh) + norm.cdf(-thresh)

コーシー分布に従う確率変数は、絶対値が1万を超える確率が0.01%もあります。したがって、例えば、1万サンプル集めて少なくとも1サンプルは\(|X| > 90\)となる確率は\(1 – (1 – 0.0001)^{1万} \approx 63\text{%}\)であり、十分あり得るということになります。t分布は裾が重いので\(|X|>10|\)の場合は0.02%ですが、他の場合では正規分布含め生起確率はほとんどゼロです。

おわりに

以上、いろんな観点からの可視化でコーシー分布の異常さを確認してみました。最後に、気になったことを箇条書きでまとめます。

• コーシー分布の活用について…莫大な外れ値がある、という特徴は地震などの自然現象ではよく見かけるのではないでしょうか。実際、（ハーフコーシー分布を含め）裾の重さを利用して階層モデルの事前分布などにもよく使用されています。そういった状況で使えそうとは思いつつ、あまりにも裾が広いためやっぱりちょっと使いにくい、ということにもなりそうで、なかなか難しそうだなと感じました。
• コーシー分布に従う確率変数の標本平均の収束について…本記事での実験では「収束しそうにない」と結論づけましたが、全体的に見ると平均値は揺れながらも0に向かっていると見えなくもないです。1兆以上の規模ではどうなるのか、ある程度分散が逓減するなどの特徴はあるのか、など、気になりました。

参考

• 数学の景色さん、「コーシー分布の定義と性質とその証明」
• Mark-Nさん、「初めてのコーシー分布」

Follow me!

@data_everyday