はじめに
統計を勉強するにあたり、
教科書の演習問題を解くだけではモチベーションが上がらないので、実際のデータに当てはめて考察してみました。
テーマは「指数分布、中心極限定理、ガンマ分布」です😀
データセットの入手
指数分布とガンマ分布は、いずれもランダムに起きる事象を扱う確率分布です。
なので、ランダムに起きていそうな現象を探しました。
Google Cloud 一般公開データセットを見ていたところ、
シカゴの犯罪に関するデータがあったので、それを使ってみることにしました。
犯罪のカテゴリーとしては、ランダムに起こりそうなものを選びたかったので、「HOMICIDE（殺人）」を選びました。
以下のリンクからGoogle Couldにログインします。
（データセットのソース ↓
https://data.cityofchicago.org/）
「chicago」で検索すると、以下のデータセットが表示されます。
「Chicago Crime Data」 をクリックします（2020/5/21現在の表示）
次に「データセットを表示」をクリックします。
画面上部の空白にSQLのクエリを書いて「実行」を押します。
すると、所望のデータだけを抽出することができます。
今回は、以下のクエリにより
2001年1月1日 ～ 最終更新（2020年5月9日15時過ぎ）までに起きた殺人事件の発生時刻を抜き出しました。トータルで10223件です。
SELECT
date
FROM
`bigquery-public-data.chicago_crime.crime`
WHERE
primary_type = "HOMICIDE"
ORDER BY
date
抽出し終わったら、「クエリを保存」を押し、csvに保存します。
ちなみに、上記の操作にはGoogle Cloud Platformというシステムを使っています。
有料ですが、12 か月間 $300 分の無料トライアルが用意されているので、この程度の操作であれば無料で使うことができます😃
指数分布
単位時間にλ回起きるランダムな現象の発生間隔を表します。
$$f(x)=λe^{-λx} (x：時間）・・・(1)$$
殺人事件がランダムに起こっているのであれば、その発生間隔が指数分布に従うはずです。先ほどのcsvを使って、殺人事件の発生間隔を算出してプロットします。
発生間隔のプロット
Pythonを使って行います。
import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import expon
from scipy.stats import gamma
import japanize_matplotlib
import seaborn as sns
sns.set(font='MS Gothic')
# データの読み込み
data = pd.read_csv("データセットから抽出したcsvのパス")
# date列から"UTC"を削除する
UTC_del = lambda x: x.replace(" UTC","")
data["date"] = data["date"].map(UTC_del)
#date列をdatetime型に変換する
data["date"] = pd.to_datetime(data["date"])
# date列でソートする（不要かも）
data = data.sort_values('date')
# date列の1階差分を取る
data_diff = data.copy()
data_diff["delta time"] = data_diff["date"].diff()
# 0行目は欠損するので除く
data_diff = data_diff.iloc[1:,:]
# delta time列を分に変換する
trans_min = lambda x: x.total_seconds()/60
data_diff["delta time (min)"] = data_diff["delta time"].map(trans_min)
# 発生間隔の最短、最長、平均を求める
print("殺人事件が起こらなかった最短時間：{}分".format(round(data_diff["delta time (min)"].min(),1)))
print("殺人事件が起こらなかった最長時間：{}分".format(round(data_diff["delta time (min)"].max(),1)))
print("殺人事件が起こらなかった平均時間：{}分".format(round(data_diff["delta time (min)"].mean(),1)))
殺人事件が起こらなかった最短時間：0.0分
殺人事件が起こらなかった最長時間：12631.0分
殺人事件が起こらなかった平均時間：995.7分
平均すると16.6時間に1件の殺人が発生しているようです。
# 殺人件数が発生する平均的な時間間隔
inv_lambda = data_diff["delta time (min)"].mean()
# 指数分布を発生させるための時間軸
x = np.arange(0, 6000, 1)
# ヒストグラムを指数分布の理論曲線を描く
plt.figure(figsize=(20, 10))
# 実データのヒストグラム
plt.hist(data_diff["delta time (min)"],
bins=int(round(np.log2(len(data_diff))+1)),
density=True,label="殺人事件の発生間隔の頻度")
# 指数分布の理論曲線（λ = 1分当たりの殺人事件発生回数の頻度）
plt.plot(x,expon.pdf(x,scale=inv_lambda),
color="red",label="λ＝[1分あたりの発生回数]とした指数分布")
plt.xlabel("発生間隔(分)",fontsize=18)
plt.ylabel("確率密度",fontsize=18)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.xlim(data_diff["delta time (min)"].min(),6000)
plt.legend(fontsize=18,loc="best")
plt.grid(b=True, which='major', axis='both')
指数分布ぽくなっていますね。
ヒストグラムを描く際の階級数には注意が必要です。
階級数の取り方で、ヒストグラムの見え方が如何様にも変わるからです💦
今回は以下のスタージェスの公式を使いました。
$$k≒1+log_2N(k：階級数、N：観測数)$$
さて、上のグラフの理論曲線（赤線）を見て何か疑問に思われないでしょう？
「なぜ、発生間隔0（分）で確率密度が最大になるのか？」と。 これは、ある殺人事件が起こった次の瞬間に殺人事件が起こる確率が最も高いという
ことです。感覚的には、殺人事件が1件が起こったら、しばらく平和な時間が続きそうじゃないですか？😅（というか、続いて欲しい）
同じような疑問を感じられた方がいましたら、是非以下の動画をご覧ください。
まさにその点が解説されています。しっくりくる説明はたくみ先生にお任せするとして、
数式の観点からこの疑問を解決します。ある殺人事件が起こってから、ｓの間殺人事件が発生しなかったという前提を置きます。
その時にさらにｔ経過しても殺人事件が発生しない確率を考えます。
この確率は以下のように書けます。
$$P(X≧s+t|X≧s)$$
$$=\frac{P(X≧s+t,X≧s)}{P(X≧s)}$$
$$=\frac{P(X≧s+t)}{P(X≧s)}$$
$$=\frac{e^{-λ(s+t)}}{e^{-λs}}$$
$$=e^{-λt}$$
$$=P(X≧t)$$
ある殺人事件が発生してから、tの間殺人事件が発生しない確率と等しくなりました。
つまり、ランダムに起こる事象に対しては「sの間事象が発生しなかったから、次のtの間も発生しにくいとか、発生しやすいということは全く言えない」ということです。
これを「無記憶性」と言います。この性質があるので、ある殺人事件が起こった次の瞬間に別の殺人事件が起こる可能性が最も高いわけです。
統計学的には、
「天災は忘れたころに忘れる前にやってくる」と言うのが正しそうです。
分布関数を求めてみます。
plt.plot(x, expon.cdf(x,scale=inv_lambda)*100, color="red",alpha=0.5, label="殺人事件の発生確率".format(round(inv_lambda,1)))
plt.xlabel("発生間隔(分)",fontsize=24)
plt.ylabel("分布関数×100",fontsize=24)
plt.xticks(fontsize=22)
plt.yticks(fontsize=22)
plt.xlim(data_diff["delta time (min)"].min(),6000)
plt.legend(fontsize=28,loc="best")
plt.grid(b=True, which='major', axis='both')
print("1時間以内に殺人が起きる確率：{} %".format(round(expon.cdf(60,scale=inv_lambda)*100,1)))
print("24時間以内に殺人が起きる確率：{} %".format(round(expon.cdf(60*24,scale=inv_lambda)*100,1)))
1時間以内に殺人が起きる確率：5.8 %
24時間以内に殺人が起きる確率：76.5 %
中心極限定理
せっかくなので、中心極限定理が成り立っているか確認します。
中心極限定理とは、
母集団がどんな分布でも、サンプルサイズが十分大きければ、
標本平均の確率分布が正規分布に従うと見なしてよいという定理です。
式で書くと以下のようになります。
$$標本平均\frac{(X_1+X_2+・・・+X_n)}{n} ～ N(μ,\frac{σ^2}{n})・・・(2)$$
殺人事件の発生間隔の母平均、母分散は不明なので、
不偏推定量である標本平均と不偏分散を使います。
標本数を100に固定して、
サンプルサイズn = 1、10、100、500についてヒストグラムと(2)式から求めた正規分布を
書くと以下のようになります。
# 殺人事件の間隔をランダムサンプリングして標本平均を求める関数
# sample_size：サンプルサイズ、number_of_samples：標本数
def cal_mean(sample_size , number_of_samples):
# 標本平均を入れるリストを作成
mu_list = []
for i in range(number_of_samples):
# indexをランダムに生成
random_index = random.sample(list(data.index), k=sample_size)
# 平均を算出
mu = data_diff["delta time (min)"][random_index].mean()
# リストに格納する
mu_list.append(mu)
return mu_list
# 標本数は100に固定する
number_of_samples = 100
# サンプルサイズ1の時の標本平均値
n1 = cal_mean(1 , number_of_samples)
# サンプルサイズ10の標本平均値
n10 = cal_mean(10 , number_of_samples)
# サンプルサイズ100の標本平均値
n100 = cal_mean(100 , number_of_samples)
# サンプルサイズ500の標本平均値
n500 = cal_mean(500 , number_of_samples)
# 時間軸
X = np.arange(0,3000,0.1)
#標本平均（＝母平均の不偏推定量）
mu = data_diff["delta time (min)"].mean()
# 不偏標準偏差（＝母標準偏差の不偏推定量）
std = data_diff["delta time (min)"].std()
# サンプルサイズ1の時の描画
Y1 = norm.pdf(X,mu,std/np.sqrt(1))
plt.plot(X,Y1,color='#123456',linestyle="-",label="N(μ,$σ^{2}$/1)")
plt.hist(n1,
bins=int(round(np.log2(len(n1))+1)),
color='#123456',
density=True,label="サンプルサイズ = 1",
alpha=0.3)
# サンプルサイズ10の時の描画
Y10 = norm.pdf(X,mu,std/np.sqrt(10))
plt.plot(X,Y10,color='#4169e1',linestyle="-",label="N(μ,$σ^{2}$/10)")
plt.hist(n10,
bins=int(round(np.log2(len(n10))+1)),
color='#4169e1',
density=True,label="サンプルサイズ = 10",
alpha=0.3)
# サンプルサイズ100の時の描画
Y100 = norm.pdf(X,mu,std/np.sqrt(100))
plt.plot(X,Y100,color='#1300e6',linestyle="-",label="N(μ,$σ^{2}$/100)")
plt.hist(n100,
bins=int(round(np.log2(len(n10))+1)),
color='#1300e6',
density=True,label="サンプルサイズ = 100",
alpha=0.3)
# サンプルサイズ500の時の描画
Y500 = norm.pdf(X,mu,std/np.sqrt(500))
plt.plot(X,Y500,color='#e6bf00',linestyle="-",label="N(μ,$σ^{2}$/500)")
plt.hist(n500,
bins=int(round(np.log2(len(n10))+1)),
color='#e6bf00',
density=True,label="サンプルサイズ = 500",
alpha=0.3)
plt.xlabel("発生間隔の平均値(分)",fontsize=24)
plt.ylabel("確率密度",fontsize=24)
plt.xticks(fontsize=22)
plt.yticks(fontsize=22)
plt.xlim(data_diff["delta time (min)"].min(),3000)
plt.legend(fontsize=28,loc="best")
plt.grid(b=True, which='major', axis='both')
中心極限定理
確かに、サンプルサイズが大きくなるにつれてヒストグラムが(2)式の正規分布に近づいていることが分かります。サンプルサイズ10ですでに峰が見え始めています。そして、サンプルサイズが増加するに従って、全データを使って算出した標本平均995.7分の辺りにデータが集中してくることが分かります。
最後に10件の殺人事件が起こるのに要する時間がガンマ分布にあてはまるかどうか見てみます。
ガンマ分布
ガンマ分布は平均的に期間\(β\)ごとに1回起きるランダムな事象が\(α\)回起こるまでに要する
時間の分布を表します。
$$Ga(α,β)$$
$$=\frac{1}{Γ(α)}\frac{1}{β}(\frac{x}{β})^{α-1}exp(\frac{-x}{β}) , x>0$$
$$ただし、Γ(α)は以下で定義される。$$
$$Γ(α)=\int^∞_0 y^{α-1}exp(-y)dy$$
ここで、\(α\)=10、\(β\)=995.7分（発生間隔の平均）です。
# dateの10階差分を取る
data_diff10 = data.copy()
data_diff10["delta time"] = data_diff10["date"].diff(periods=10)
# 10行目までは欠損するので削除する
data_diff10 = data_diff10.iloc[10:,:]
# delta time列をトータル分に変換する
trans_min = lambda x: x.total_seconds()/60
data_diff10["delta time (min)"] = data_diff10["delta time"].map(trans_min)
# 時間軸
x = np.arange(0, 40000, 1)
plt.figure(figsize=(20, 10))
# 殺人事件が10件発生するのに要する時間
plt.hist(data_diff10["delta time (min)"],
bins=int(round(np.log2(len(data_diff10))+1)),
density=True,label="殺人事件10件の発生間隔の頻度")
# ガンマ分布の理論曲線
plt.plot(x, gamma.pdf(x,a=10,scale=inv_lambda),
color="red",label="Ga(α,β)")
plt.xlabel("殺人事件10件の発生間隔(分)",fontsize=24)
plt.ylabel("確率密度",fontsize=24)
plt.xticks(fontsize=22)
plt.yticks(fontsize=22)
plt.xlim(data_diff10["delta time (min)"].min(),30000)
plt.legend(fontsize=28,loc="best")
plt.grid(b=True, which='major', axis='both')
ガンマ分布とヒストグラムの比較
うーん、結構ガンマ分布から外れています。
（ということは、指数分布からもそれなりに外れていたということか。）
おそらく、\(β\)= 一定としてしまっていることに無理があるのではないかと思います🤔
は2001年1月1日 ～ 2020年5月9日の間の平均的な殺人事件の発生間隔であり、ランダムにこの頻度で殺人事件が発生していれば、ヒストグラムとGa(α,β)は一致するはずです。
しかし、実際には上記の期間には様々なことが起こっています。
9.11テロ、リーマンショック、大統領の交代 etc.
殺人事件の発生に影響を及ぼす経済状況の変化や法律の改正があったかもしれません。
したがって、完全にランダムに殺人事件が発生していたとは言えないように思います。
つまり、殺人事件が起きやすい期間と起きにくい期間があったのではないかと。
数式で表すと、
$$Ga(α,β(t))$$
$$=\frac{1}{Γ(α)}\frac{1}{β(t)}(\frac{x}{β(t)})^{α-1}exp(\frac{-x}{β(t)}) , x>0$$
βが刻々と変わる時間の関数なのではないかと。
この\( β(t) \)の出し方は勉強が進んだら検討してみたいと思います。
（そもそも数学的に関数を導出できるのか否かすら分かりません😅）
そもそも、ランダムであることって数学的に厳密に定義するとどうなるんだろう🙄？
最後に
統計を勉強をするにあたっては、世の中の現象がどんな確率分布に従っているか見てみる
と面白いです。教科書どおりの数式にピッタリあてはまることは少ないかもしれせん。でも、逆にそれが勉強意欲をそそってくれますｗ
オープンデータが豊富にあり、PythonやRでサクッと可視化できる時代になったのでこれらを勉強に活用しない手はないですね🤗
