40歳PMがAI初学者としてデータ分析をやってみた話

はじめに

初めまして。PM（プロジェクトマネジメント）歴10年以上のAI初学者です。
縁あってAIのベンチャー企業に入社しPMの役割でお仕事させていただいています。クライアント業務なこともありお客様とAI技術、サービスについてお話しする機会はたくさんありますので、勉強兼ねてAIの勉強を始めました。

本記事の概要

国分寺に住み続けたいが心配がある。その心配とは「人口減少で市が衰退していかないか」この問いに答えるべく国分寺市の人口推移を予測します。

どんな人に読んで欲しいか

国分寺が好きな人この記事に書くこと、

わかること

国分寺の人口の予測モデルをSARIMA、LSTMの２つのモデルを構築した結果の精度比較

この記事で扱わないこと

上記以外の全て

やったことの全体の流れ

• データの取得
• データの整理
• データの分析と可視化
• SARIMAモデル構築と推論結果
• LSTMモデル構築と推論結果

データ取得

データの取得元は国分寺市のHPです。e-statのデータが扱いやすいと思ったのですが、市区町村単位のデータはなかったので諦めました。
PDFだったのでcsvにするのが面倒でした。数日に分けて10時間弱かかりました。PDFを見ながら手動でcsvデータに変換していきました。
ちなみに1986年から2023年までのデータです。正直イヤになりました。

データの整理

csvを読み込みデータを整形します。必要データを追加、あるいはそのまま残し、不要データを削除します。今回の人口予測なので、世帯数や男性人数、女性人数など不要なものは削除して、総人口を新たにデータ列として加えます。

import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

df_kokubunji = pd.read_csv('/content/drive/MyDrive/PopulationOfKokubunji.csv')

index = pd.date_range('1986-01-01', '2023-12-31', freq='MS')
df_kokubunji.index = index

del df_kokubunji["era"]
del df_kokubunji["month"]
del df_kokubunji["family_fo"]
del df_kokubunji["family_mix"]
del df_kokubunji["male_fo"]
del df_kokubunji["female_fo"]

df_kokubunji["total"] = df_kokubunji["male_jp"] + df_kokubunji["female_jp"]


df_kokubunji = df_kokubunji
del df_kokubunji["family_jp"]
del df_kokubunji["male_jp"]
del df_kokubunji["female_jp"]

データの分析と可視化

データの傾向を掴む

季節調整を行い原系列をトレンド、季節変動、残差に分けて出力します。
下の図でTrendをみると上昇しているのは一目瞭然です。

fig = sm.tsa.seasonal_decompose(df_kokubunji, period = 30).plot()
plt.show()

自己相関係数と偏自己相関係数

コレドグラムで過去のデータと現在のデータに関係性があるかを確認しました。
年間のソフトクリーム売り上げのような季節と売上に相関があるような特徴は見られませんでした。一方的に上昇しているだけの人口なので季節はあまり関係ないか。

#自己相関係数
fig = plt.figure(figsize =(16,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df_kokubunji, lags=100, ax=ax1)

#偏自己相関係数
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df_kokubunji, lags=65, ax=ax2)
plt.show()

季節性は関係なさそうとはいえ、時系列データを扱う上ではSARIMAモデルを使うのは王道らしいのでいったんSARIMAモデルで推論します。

SARIMAモデル構築と推論結果

SARIMAモデルの適切なパラメータを求めます。

import itertools

# 1年間分のデータにしています
df_data = df_kokubunji.tail(12)

# 関数の定義
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(df_data,12)
結果、[(0, 1, 0), (0, 1, 0, 12), nan]と出ましたのでモデル構築でこの値を使います。以下のコードで推論結果を見てみます。

#訓練データとテストデータの分割
def split_train_test(data, train_ratio):

    # 学習データのサイズ
    train_size = int(len(data) * train_ratio)

    # 順番に分割
    train = data[:train_size]
    test = data[train_size:]
    return train, test

train_data, test_data = split_train_test(df_kokubunji, 0.8)

# SARIMAモデルの構築と学習
SARIMA_kokubunji = sm.tsa.statespace.SARIMAX(
    train_data,
    order=(0, 1, 0),
    seasonal_order=(0, 1, 0, 12)
).fit()

# BICの出力
print("BIC:", SARIMA_kokubunji.bic)

# 予測
pred = SARIMA_kokubunji.predict(start=test_data.index[0], end=test_data.index[-1])

# 予測データと実際のデータの可視化して評価
plt.figure(figsize=(12, 6))
plt.plot(train_data, label="Train")
plt.plot(test_data, label="Test")
plt.plot(pred, color="r", label="Prediction")
plt.legend()
plt.show()

# 精度評価 - 定量的に予実の当てはまりの良さを評価
mae = mean_absolute_error(test_data, pred)
mse = mean_squared_error(test_data, pred)
rmse = np.sqrt(mse)
r2 = r2_score(test_data, pred)

print("Mean Absolute Error (MAE):", mae)
print("Mean Squared Error (MSE):", mse)
print("Root Mean Squared Error (RMSE):", rmse)
print("R-squared (R2):", r2)

RMSEがよくわからない大きな値になっています。推論結果のグラフもPredictionとTestが全然重ならず・・・結果がよさそうに見えません。

LSTMモデル構築と推論結果

LSTMで学習するには、

• データを標準化する
• 次元を増やす
• データセットを時間ステップ（シーケンス長）で作成

など面倒があり以下のように実装し学習します。

#LSTMを使った予測
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import StandardScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense

#訓練データとテストデータに分割
train, test = split_train_test(df_kokubunji, 0.8)

# データの標準化
scaler = StandardScaler()
scaled_train = scaler.fit_transform(train)
scaled_test = scaler.transform(test)

# データセットの作成
def create_dataset(dataset, time_step=1):
    X, y = [], []
    for i in range(len(dataset)-time_step-1):
        a = dataset[i:(i+time_step), 0]
        X.append(a)
        y.append(dataset[i + time_step, 0])
    return np.array(X), np.array(y)

# 時間ステップ（シーケンス長）の設定
# 過去の20つのデータポイントを使って次の値を予測する
time_step = 20
X_train, y_train = create_dataset(scaled_train, time_step)
X_test, y_test = create_dataset(scaled_test, time_step)

# データの整形 [サンプル数, 時間ステップ, 特徴数]
X_train = X_train.reshape(X_train.shape[0], X_train.shape[1], 1)
X_test = X_test.reshape(X_test.shape[0], X_test.shape[1], 1)


# LSTMモデルの構築
model = Sequential()
# LSTMモデルのサイズの指定 データサイズの指定にはinput_shape= (行数, 列数) を用います
model.add(LSTM(input_shape=(time_step, 1), units=50, return_sequences=True))
model.add(LSTM(units=50, return_sequences=False))
model.add(Dense(1))

#model.compile(optimizer='adam', loss='mean_squared_error', metrics=['accuracy'])
model.compile(optimizer='adam', loss='mean_squared_error')

# モデルの訓練
model.fit(X_train, y_train, epochs=30, batch_size=1, verbose=1)

モデルが作成できたので、訓練データを使って予測してみます。
「訓練データ」で訓練したのだからもちろん良い結果だろという疑問がありますが、これは正しく学習されていることを確認したいのでやっています。精度をみたいわけではないです。

その後、テストデータで推論し、結果を表示します。

 訓練データで予測の作成
train_predict = model.predict(X_train)

# inverse_transformで、同じ平均・分散を使って調整したデータを元に戻せます。
train_predict = scaler.inverse_transform(train_predict)
y_true = scaler.inverse_transform(y_train.reshape(-1, 1))


# 結果のプロット
fig = plt.figure()
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2)

#plt.figure(figsize=(8,4))
ax1.plot(train.index[time_step+1:], y_true, label='True Population')
ax1.plot(train.index[time_step+1:], train_predict, label='Predicted Population')
ax1.legend()

# テストデータ予測の作成
test_predict = model.predict(X_test)

# inverse_transformで、同じ平均・分散を使って調整したデータを元に戻せます。
test_predict = scaler.inverse_transform(test_predict)
X_test_true = scaler.inverse_transform(y_test.reshape(-1, 1))

# 結果のプロット
#plt.figure(figsize=(8,4))
ax2.plot(test.index[time_step+1:], X_test_true, label='True Population')
ax2.plot(test.index[time_step+1:], test_predict, label='Predicted Population')
ax2.legend()
plt.show()

score = model.evaluate(X_test, y_test, verbose=0)

# 評価結果の出力
if isinstance(score, (list, tuple)):
    print('Validation loss:', score[0])
    print('Validation accuracy:', score[1])
else:
    print('Validation loss:', score)

compileでパラメーターとしてmetricsを設定しないのは、これは２値分類の問題ではないからです。人口が増加または減少の２値分類であればmetricsでAccuracy、RecallやPrecisionなど見るべきですが、この問題では未来の数値予測なので実際の値と推論した値がどれだけ離れているかを示す数値で精度をみます。
実行すると以下が得られます。

上：訓練データを推論
下：テストデータを推論

上の折れ線グラフで確かに学習はできてそうだとわかります。
下の折れ線グラフは、なんかシェイプは追従していていい感じです。Lossの値も小さいので良さげです。

今後の活用

エンジニアばかりに頼らずに自身でもデータ分析をしてみたい、簡易的にモデルを作ってみたいと思いました。

おわりに

実際にモデルを構築してみると、モデル開発の全体の規模や流れがより理解できました。
基礎は大事ですが、あまり正直に真面目に100%理解しようとせずわかったふりで早く何かを生み出すことにパワーを注いだ方が結果的に効率よい学びが得ることができます。
これである程度、言葉の意味や実作業がわかったのでAIエンジニアと少し会話がスムーズに、相互理解が深まるコミュニケションができそうだと思いました。これからはエンジニアばかりに頼らずに自身でもデータ分析をしてみたいと思います。

補足
なおこのブログはお世話になったAidemy Premiumのカリキュラムの一環で、受講修了条件を満たすために公開しています