Pythonによるデータ処理16 ～カルマンフィルタ3

今回はオープンデータの時系列データに対して前回までのカルマンフィルタを適用した例を紹介する。

atatat.hatenablog.com

1. オープンデータの入手と加工
2. csvデータの読み込み
3. カルマンフィルタの適用
4. 日付の追加
5. まとめ

1. オープンデータの入手と加工

以下の政府統計データから1953年から今までの完全失業者の統計データをcsvで入手する。

データが大量であればpandasなどを用いて加工する必要があるが、今回は手作業で十分に対応できるので時点(年月日)と完全失業者のみの2列のデータとして、さらに1か月ごとのデータの下に半期ごとなどのデータが含まれているため、手動で削除しファイル名をtest_employee.csvとした。

2. csvデータの読み込み

　まずGoogle Colabにcsvファイルをアップロードするには、サイドにあるフォルダアイコンをクリックする。

つぎに矢印のアップロードのアイコンをクリックして、ファイルをGoogle Colabにアップロードする。

アップロード後は以下のようになる。

アップロードしたファイルはcontent/に格納されている。以下xは番号、pointは日付(文字列)、dataは完全失業者数(整数値)を格納する。

import numpy as np
from matplotlib import pyplot as plt
import csv

csv_file = open('/content/test_employee.csv')
f = csv.reader(csv_file)
header = next(f)
x=[]
data=[]
point=[]
i=1
for row in f:
  x.append(i)
  point.append(row[0])
  data.append(int(row[1]))
  i=i+1

以下のようにpointは日付の文字列、dataは数値が格納されていることが確認できる。

point[0]
>'1953/1/1'
data[0]
>79

3. カルマンフィルタの適用

csvから読み込んだdataについてカルマンフィルタを適用する。前回/前々回と基本的には同じコードを使用するが、dataのサイズのみ異なり、結果を格納するarrayサイズを変更するため例えば以下のようにする。ちなみに今回のdataのサイズは831となっている。

T = np.size(data)
m0 = np.array([[0]])
C0 = np.array([[1e7]])

# 結果を格納するarray
m = np.zeros((T, 1))
C = np.zeros((T, 1, 1))
s = np.zeros((T, 1))
S = np.zeros((T, 1, 1))

完全失業者のデータ(オリジナル)とカルマンフィルタ、カルマン平滑化のプロットが以下のように得られる。オリジナルデータに対してカルマンフィルタは滑らかになっていることが確認できる。

"完全失業者統計データのカルマンフィルタ適用" — 完全失業者統計データのカルマンフィルタ適用

カルマンフィルタとカルマン平滑化がわかりにくいので、後半一部を拡大した。前回と同様にカルマンフィルタよりもカルマン平滑化の方が滑らかになっている。

4. 日付の追加

上のプロットの横軸は番号としたため、いつの失業者数か把握しにくい。そこで横軸を日付に変更するが、なるべく節目で表示数もプロット上から見やすいように数を表示数を制限する。

流れとしてはx軸ラベルを付与するプロットの番号を指定し、指定した番号を日付の文字列に置き換える。ただし、元の日付は2022/06/01と冗長なため年のみの形式2022のようにする。

まずx軸ラベルを付与するための配列new_xticksを用意する。データファイルは月ごとの数値なので、例えば48番ごとにラベルを付与する場合は4年ごとにラベルすることになる。また最後のデータの日付を表示させたいので831番だけnp.appendを用いて最後に追加した。

new_xticks=np.arange(0,830,48)
new_xticks=np.append(new_xticks,831)

次にnew_xticksの番号に対応した日付の文字列を含むを生成する。このときに年月日から年のみの形式に変換した文字列になるようにする。

これを実行するにはdatetimeモジュールのstrptime(文字列をdatetime型に変換)を用いてnew_xticksに含まれる番号の日付をpointから読み出しdttに格納する。このときpointに格納されている日付形式を%Y/%m/%dとして指定する。次にstrftime(datetime型を文字列に変換)関数を用いて、形式を%Yとすることで年のみの日付の文字列に変換し、new_labelに格納する。

import datetime
new_label=[]
for i in new_xticks:
  dtt=datetime.datetime.strptime(point[i],'%Y/%m/%d')
  new_label.append(datetime.datetime.strftime(dtt,'%Y'))

上記で生成した、ラベル表示する番号を指定しているnew_xticksと番号に対応した日付文字列new_labelを用いて以下のようなコードでプロットの横軸を日付にすることができる。

参考元と詳細な説明は qiita.com

from matplotlib.dates import DateFormatter
import matplotlib.ticker as ticker
figure_ = plt.figure(1, figsize=(14,6))
axes_ = figure_.add_subplot(111)
xaxis_ = axes_.xaxis
xaxis_.set_major_locator(ticker.FixedLocator(new_xticks))
xaxis_.set_ticklabels(new_label)
axes_.plot(x, data, '--')
axes_.plot(x,m,'r')
axes_.plot(x,s,'m')
axes_.set_ylabel('Umemployed ( x$10^{4}$ )')