コロナ禍での犯罪者数を予測してみた!!

目次

1.はじめに
2.本題
3.時系列解析(SARIMAモデル)
4.教師ありの機械学習(回帰分析)
5.まとめ
6.おわりに

1. はじめに

私の自己紹介をしたいと思います。
機械学習については、仕事で少し知識はありますが、プログラミングについては初心者です。

2. 本題

今回は月別の犯罪者数を予測してみることにしました。コロナによって、犯罪者数にどんな影響があるかも調査してみました。時系列モデルのsarimaモデルで予測を行います。
　使用するデータはe-Statという政府の統計を管理しているページからダウンロードしました。この中の「罪名別被疑事件の処理人員(2007年1月～2023年2月)」の数字を扱います。
※「犯罪認知件数」という言葉の方が一般的かもしれませんが、
　 犯罪を件数単位ではなく、人数単位でカウントしたかったのでこの統計を選びました。
　時系列モデルのsarimaで予測した結果、犯罪者数がコロナが原因で減少したのを予測出来ていましたが、極端に減少した箇所は予測出来ていなかったので、時系列データ以外の下記説明変数も使用し、犯罪者数が予測出来るかも調査しました。
失業者数、自殺者数、来訪者数、日経平均株価、コロナ患者数、物価指数

・実行環境
Jupyter Notebook
Python 3.7.3

機械学習の手法
●時系列解析(SARIMAモデル)
●教師あり機械学習(線形回帰)

3. 時系列解析(SARIMAモデル)

時系列データは 時間の経過とともに変化するデータ のことをいいます。具体的には会社の売り上げや商品の売り上げ予測、さらには来店者数の予測などで、時系列分析はビジネスにおいて非常に重要な分析技術となります。時系列分析の第一歩は、時系列データの視覚化です。

自己相関係数
過去の自分の値との相関を現した値。一方の時系列データをずらして相関係数を求めます。

偏自己相関
7日間差の自己相関係数を求めたとします。しかし、ある日のデータが前日のデータと相関があった場合、7日前→6日前→5日前......1日前→今日とデータを通して相関している可能性があります。この間の影響を取り除いて相関を求めたもの

実際の時系列データを観測すると、下記の3つのパターンが組み合わさっている。
●トレンド
データの長期的な傾向を意味します。時間の経過とともに、値が上昇している場合には 正のトレンド 下降している場合には 負のトレンド があると言います。

●周期変動
データは時間の経過に伴ってデータの値が上昇と下降を繰り返します。特に1年間での周期変動を 季節変動 といいます。

●不規則変動(残差)
時間の経過と関係なくデータの値が変動することをいいます。

時系列分析のフロー
①データを取得
②データをグラフ化
③データの整理
④モデルの決定
⑤SARIMAモデルの実行

①データを取得
e-Statからページからダウンロードした罪名別被疑事件の処理人員(2007年1月～2023年2月)」のデータを読み込む

import glob
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
import datetime
import itertools

年別にファイルが分かれているので、1つのフォルダにファイルを格納して、glob関数でファイルを取得

file_paths = glob.glob('/content/drive/MyDrive/成果物/v5909/*.xls*')
file_paths

取得したファイルを読み込む

yearDate_sumClim_list = []
for file_path in file_paths:
  df = pd.read_excel(file_path,usecols=[3], skiprows=3)
  df = df.dropna()
  犯罪者数 = df.iloc[0,0]
  file_name = file_path.split('/')[-1]
  yearDate_sumCrime_list.append([file_name, 犯罪者数])

merge_df = pd.DataFrame(yearDate_sumCrime_list,columns=['年月', '犯罪者数'])
merge_df

merge_df=merge_df.replace({'.xls': '', '.xlsx':''},regex=True)
merge_df=merge_df.replace('x', '',regex=True)
merge_df

年月と犯罪者数のDataFrame作成
年月 犯罪者数
0 07-02-02 183207.0
1 07-01-02 169416.0
2 07-03-02 190646.0
3 07-04-02 167558.0
4 07-05-02 188723.0
... ... ...
187 22-09-02 67621.0
188 22-08-02 67716.0
189 22-10-02 75177.0
190 22-11-02 79348.0
191 22-12-02 91676.0
192 rows × 2 columns

ファイル名に-2と付いていたので削除し、年月列の文字型を日付型に変換

merge_df['年月'] = merge_df['年月'].map(lambda x: x[:-3])
merge_df = merge_df.sort_values('年月')
merge_df['年月'] = pd.to_datetime(merge_df['年月'], format='%y-%m')
merge_df=merge_df.reset_index()
merge_df=merge_df.drop('index',axis=1)
merge_df

yearDate列をindexに設定し、コロナ期間のデータを削除したDFを作成

merge_df=merge_df.set_index("年月")
merge_df.index=pd.DatetimeIndex(merge_df.index.values,freq=merge_df.index.inferred_freq)
merge_df_nondel=merge_df
merge_df=merge_df.drop(merge_df.iloc[156:].index)
merge_df

②データをグラフ化
2020年はコロナの影響で減少が大きいですが、それ以降はコロナ前と同様に減少しています。

③データの整理
・自己相関係数

# 自己相関を求める
df_acf = sm.tsa.stattools.acf(merge_df['犯罪者数'], nlags=40)

#  自己相関のグラフ
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_acf(merge_df['犯罪者数'], lags=40)
plt.show()

青色の影は95%信頼区間です。
この点線の中におさまっていれば自己相関を持たない、つまり現在の自分と過去の自分は関係ないということです。12ヶ月前までは自己相関がある、大きく影響を受けていることがこのグラフからわかります。

偏自己相関

# 偏自己相関を求める
df_pacf = sm.tsa.stattools.pacf(merge_df['犯罪者数'], nlags=20, method='ols')

# 偏自己相関を可視化する
fig = plt.figure(figsize=(12,8))
fig = sm.graphics.tsa.plot_pacf(merge_df['犯罪者数'], lags=20)

この結果を見て、左から0～3,9,10,12,13,17番目のデータが影響を与えていることがわかります。

・データパターンの確認
sm.tsa.seasonal_decompose()を使用すると原系、トレンド、 季節変動、 残差に簡単に分けることができます。

res = sm.tsa.seasonal_decompose(merge_df["犯罪者数"], period=12, extrapolate_trend='freq')
plt.show()
fig = res.plot()

④モデルの決定
どのモデルで分析、予測していくかを決める必要があります。ARMAモデルは定常データ、ARIMAモデルは非定常データ、SARIMAモデルは非定常で季節性があるデータの解析に適しています。
上記のsm.tsa.seasonal_decompose()で取得したデータは上から原系、トレンド、季節性、残差です。
トレンドと季節性があるのが分かったので、今回はSARIMAモデルを使用します。また、季節性があるのでs=12とします。
(p, d, q)(sp, sd, sq, s)の７つのパラメーターを決めますが、自動で最も適切にしてくれる機能はないので、情報量規準 (今回の場合は BIC(ベイズ情報量基準) )によってどの値が最も適切なのか調べるプログラムを書かなければなりません。BICの場合は この値が低ければ低いほどパラメーターの値は適切です。

import itertools
import statsmodels.api as sm
def selectparameter(DATA,s):
   p = d = q = range(0, 2)
   pdq = list(itertools.product(p, d, q))
   seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
   parameters = []
   BICs = np.array([])
   for param in pdq:
       for param_seasonal in seasonal_pdq:
           try:
               mod = sm.tsa.statespace.SARIMAX(DATA,
                                           order=param,
                                           seasonal_order=param_seasonal)
               results = mod.fit()
               parameters.append([param, param_seasonal, results.bic])
               BICs = np.append(BICs,results.bic)
           except:
               continue
   return parameters[np.argmin(BICs)]
# 周期を埋めて
selectparameter(merge_df,12)

⑤SARIMAモデルの実行

SARIMA_numberof_crime = sm.tsa.statespace.SARIMAX(merge_df, order=(0, 0, 0),seasonal_order=(1, 1, 1, 12)).fit() 
#構築したSARIMAモデルのBICを出力します
print(SARIMA_numberof_crime.bic)

2008年から2025年の犯罪者数を予測します

pred = SARIMA_numberof_crime.predict('2008-01-01','2025-12-01')
pd.plotting.register_matplotlib_converters()
plt.plot(pred)
plt.show()

2020年までは減少しますが、2020年以降はほぼ横ばいです。

plt.plot(merge_df)
plt.plot(pred, color="r")
plt.show()

元の時系列データと予測データを同時に出力して予測が正しいか確認してみる。（赤色が予測)
image.png

コロナ期間(2020/1～2022/12)を拡大してプロットしました

comp = pd.merge(merge_df_nondel,pred, left_index=True, right_index=True, how="outer")
comp.columns= ["コロナ禍","予測（コロナがなかったら）"]
print(comp)

predict = pd.DataFrame(data = pred[145:193], columns =["犯罪者数"])
print(predict)
comp = pd.merge(merge_df_nondel, pred, left_index=True, right_index=True, how="outer")
comp.columns= ["コロナ禍","予測（コロナがなかったら）"]
comp_f=comp[144:195]

x = comp_f.index
y1 = comp_f["コロナ禍"]
y2 = comp_f["予測（コロナがなかったら）"]

plt.figure(figsize=(20,10))   
plt.title("Number of foreigners to Japan",fontsize=30)
plt.xlabel("Date",fontsize=20)
plt.ylabel("Numbers(In Millions)",fontsize=30)
plt.grid(True)
plt.plot(x,y1, color="b",label="REAL")
plt.plot(x,y2, color="r",label="NO COVIT19")
plt.legend()
plt.show()

2020年以降はコロナが原因で実際の犯罪者数は減少しているが、2020/1-2022/12のコロナ期間を削除したデータで予測した結果では横ばいなので、コロナが起きなければ2020年以降は横ばいだったのではと思われる。

このグラフは2020/1-2022/12のコロナ期間を削除していないデータで予測した結果で、外出制限で極端に外出が減少した際の犯罪者数以外は結構予測出来ている。

SARIMAモデルでは予測出来ていない箇所があるため、他に予測出来るモデルがないか考えました。SARIMAモデルは、外因性変数を考慮することできないが、回帰モデルは、外因性変数を考慮することができるため、外部の要因が時系列データに影響を与える場合には回帰モデルが適していると考え、回帰分析を行いました。

4. 教師ありの機械学習(回帰分析)

下記の項目が犯罪件数と相関ありそうなので、説明変数として使用
a失業者数　多いと犯罪件数も増える
b来日者数　直接は関係ないが、景気が良いと増える
c自殺者数　直接は関係ないが、景気が悪いと増える
d物価指数　大きいと犯罪件数も増える
e日経平均　小さいと犯罪件数も増える

分析のフロー
①犯罪者のデータを取得
②失業者数、来日者数、自殺者数、物価指数、日経平均データを取得
④それぞれのデータを合体
⑤学習データとテストデータに分割
⑥特徴量の作成とラベリング
⑦予測モデルを構築し、予測精度を計測
⑧実測値と予測値を同じグラフにプロット

①データの取得
a失業者数

#データの読み込み、日付の削除
df_shitugyou = pd.read_excel('/content/drive/MyDrive/成果物/shitugyou.xlsx', sheet_name=0, usecols=[0, 1, 19])

df_shitugyou=df_shitugyou.iloc[657:851] 
df_shitugyou = df_shitugyou.dropna()

df_shitugyou.columns = ['年', '月', '失業者数']
df_shitugyou['月'] = df_shitugyou['月'].map(lambda x: x[:-1])
df_shitugyou['年月'] = df_shitugyou['年'].astype(str)+ '-' + df_shitugyou['月'].astype(str)

df_shitugyou['年月'] = pd.to_datetime(df_shitugyou['年月'], format='%Y-%m')
df_shitugyou = df_shitugyou.drop(['年',"月"], axis=1)

df_shitugyou=df_shitugyou.set_index("年月")
print(df_shitugyou)

b来日者数

sumVisitor_list = []
for j in range(4):
  df = pd.read_excel('/content/drive/MyDrive/成果物/since2003_visitor_arrivals.xlsx',sheet_name=[j], skiprows=[0,1,2], usecols=[2, 4, 6, 8, 10, 12,14,16,18,20,22,24])
  sum_visitor = df[j].iloc[0]
  for j in sum_visitor:
    sumVisitor_list.append([j])
for i in range(4,17,1):
  df = pd.read_excel('/content/drive/MyDrive/成果物/since2003_visitor_arrivals.xlsx',sheet_name=[i], skiprows=[0,1,2], usecols=[1, 3, 5, 7, 9, 11,13,15,17,19,21,23])
  sum_visitor = df[i].iloc[0]
  for i in sum_visitor:
    sumVisitor_list.append([i])

Year=range(2023, 2006, -1)
Month=range(1, 13, 1)

yearDate_v_list=[]
for Y in Year:
  for M in Month:
    yearDate=(str(Y)+ '-' + str(M))
    yearDate=yearDate
    yearDate_v_list.append([yearDate])

sumVisitor_df = pd.DataFrame(sumVisitor_list)
yearDate_v_df = pd.DataFrame(yearDate_v_list)
merge_df_v = pd.merge(yearDate_v_df, sumVisitor_df, left_index=True, right_index=True, how="outer")
merge_df_v.columns=['年月','来日者数']
merge_df_v['年月'] = pd.to_datetime(merge_df_v['年月'], format='%Y-%m')
merge_df_v = merge_df_v.sort_values('年月')
merge_df_v=merge_df_v.dropna(subset=['来日者数'])
merge_df_v=merge_df_v.set_index("年月")
print(merge_df_v)

c自殺者数

# ライブラリのインポート
import datetime as dt
# 日付のリストを作成
df_jisatu = pd.read_excel('/content/drive/MyDrive/成果物/jisatusha.csv.xlsx', sheet_name=0, header=None, usecols=[0], skiprows=[1, 3, 5, 7,9,11,13,15,17,19])
df_jisatu.iloc[0,0]

Year=range(2022, 2012, -1)
Month=range(1, 13, 1)

yearDate_list=[]
for Y in Year:
  for M in Month:
    yearDate=(str(Y)+ '-' + str(M))
    yearDate_list.append([yearDate])

rateSuicide_list = []
row=range(10)#0-8
i=range(10)#0-11
for r in row:
    df_jisatu_1 = df_jisatu.iloc[r,0]
    rateSuicide= df_jisatu_1.split(' ')
    for n in rateSuicide [0:12]:
      n= n.replace(',', '')
      rateSuicide_list.append([n])

yearDate_df = pd.DataFrame(yearDate_list)
yearDate_rateSuicide_df = pd.DataFrame(rateSuicide_list)
#df_stock['月'] = df_shitugyou['月'].map(lambda x: x[:-1])

yearDate_rateSuicide_df = yearDate_rateSuicide_df.astype(float)
merge_df_jisatu = pd.merge(yearDate_df, yearDate_rateSuicide_df, left_index=True, right_index=True, how="outer")
merge_df_jisatu.columns=['年月','自殺者数']
merge_df_jisatu['年月'] = pd.to_datetime(merge_df_jisatu['yearDate'], format='%Y-%m')
merge_df_jisatu = merge_df_jisatu.sort_values('年月')
merge_df_jisatu =merge_df_jisatu.set_index("年月")
print(merge_df_jisatu)

d物価指数

#データの読み込み、日付の削除
df_bukka = pd.read_excel('/content/drive/MyDrive/成果物/bukkashisuu.xlsx', sheet_name=0, usecols=[8, 12])
df_bukka=df_bukka.iloc[457:650] 
df_bukka = df_bukka.dropna()

df_bukka.columns = ['年月', '物価指数']
#df_stock['月'] = df_shitugyou['月'].map(lambda x: x[:-1])
#df_stock['年月'] = df_stock['年'].astype(str)+ '-' + df_shitugyou['月'].astype(str)

df_bukka['年月'] = pd.to_datetime(df_bukka['年月'], format='%Y年%m月')
#df_stock = df_shitugyou.drop(['年',"月"], axis=1)

df_bukka=df_bukka.set_index("年月")
print(df_bukka)

e日経平均

#データの読み込み、日付の削除
df_stock = pd.read_excel('/content/drive/MyDrive/成果物/nikkei_stock_average_monthly_jp.xls', sheet_name=0, usecols=[0,1])
df_stock=df_stock.iloc[84:278] 
df_stock = df_stock.dropna()

df_stock.columns = ['年月', '日経平均']
#df_stock['月'] = df_shitugyou['月'].map(lambda x: x[:-1])
#df_stock['年月'] = df_stock['年'].astype(str)+ '-' + df_shitugyou['月'].astype(str)

df_stock['年月'] = pd.to_datetime(df_stock['年月'], format='%Y/%m月')
#df_stock = df_shitugyou.drop(['年',"月"], axis=1)

df_stock=df_stock.set_index("年月")
print(df_stock)

上記５つを説明変数として、犯罪者数を予測

merge_df= pd.merge(merge_df_nondel, df_shitugyou, left_index=True, right_index=True, how="outer")
merge_df= pd.merge(merge_df, merge_df_v, left_index=True, right_index=True, how="outer")
merge_df= pd.merge(merge_df, merge_df_jisatu, left_index=True, right_index=True, how="outer")
merge_df= pd.merge(merge_df, df_bukka, left_index=True, right_index=True, how="outer")
merge_df= pd.merge(merge_df, df_stock, left_index=True, right_index=True, how="outer")
merge_df= merge_df.dropna()
print(merge_df)
X = merge_df.drop('犯罪者数', axis=1)
y = merge_df['犯罪者数']

#データの分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

#訓練、評価
model =  LinearRegression()
model.fit(X_train, y_train)
R2 = model.score(X_test, y_test)

#相関係数
print("{:.5f}".format(R2))

# テスト用データを用いて学習結果をpred_yに格納
pred_y = model.predict(X_test)
pred_y = pd.DataFrame(data = pred_y, columns =["犯罪者数"], index=y_test.index)

結果R2が0.54799と精度はあまり良くない

⑧実測値と予測値を同じグラフにプロット

comp = pd.merge(y_test, pred_y, left_index=True, right_index=True, how="outer")
comp = comp.sort_index()
comp.columns= ["実際","予測"]
comp_f=comp
x = comp_f.index
y1 = comp_f["実際"]
y2 = comp_f["予測"]

plt.figure(figsize=(20,10))   
plt.title("Number of Crime in Japan",fontsize=30)
plt.xlabel("Date",fontsize=20)
plt.ylabel("Numbers",fontsize=30)
plt.grid(True)
plt.plot(x,y1, color="b",label="REAL")
plt.plot(x,y2, color="r",label="PREDICT")
plt.legend()
plt.show()

極端に減少したり増加している箇所はやはり予測出来ていない

5. まとめ

犯罪者数は年々減少していたが、2020年以降はもしコロナが起きなければそれ以上は減少せず横ばいだったと思われる。実際は、コロナの影響で外出もできなかったので、犯罪者数も減少した。
学習ありの機会学習ではR2は0.58で精度はあまり良くなかった。外出が出来ないと犯罪者数には即時に影響するが、景気に関連している指標には即時的には影響せずに、時間をおいて影響するのかもしれない。

6. おわりに

講座を受講しているだけでは、理解できていない事がたくさんあったのを実感出来た。特にデータ取得、データ加工方法は色々なパターンがありとても勉強になった。