機械学習による日本の人口推移予想

みなさん、初めまして!Aidemy研修生の佐野です。

近年、よく「機械学習」という言葉も耳にする機会も多いのではないでしょうか?

そんな流行りの「機械学習」を用いて、これからの日本の人口推移を予想していこうと思います。

さて、なぜ日本に人口推移を予想しようかと考えたのかというと、今の日本社会のままだと多くの課題を解決せぬまま、少子高齢化社会の波に飲まれてしまう未来が頭に浮かぶからです。今、機械学習に多くの人の関心が向いているため、それを利用した社会の問題(人口問題)を可視化したいからです!

ああ、熱くなりすぎました…

では、気を取り直して早速、日本の人口推移を予想していこうと思います。次のような流れで説明していきます。

目次

  1. 実行環境
  2. データについて
  3. データの収集
  4. import内容
  5. データの取得
  6. データの整理
  7. SARIMAモデル
  8. SARIMAモデル(データ標準化後)
  9. 状態空間モデル
  10. 考察
  11. まとめ

1. 実行環境

まずは、実行環境ですね。私のPCはこんな感じです。

  • OS : windows8.1
  • CPU : Intel(R) Core(TM) 2.20GHz 
  • メモリ: 4GB

お世辞にも実行環境は良いとは言えませんが、頑張っていきます!

2. データについて

今回、日本の総人口の推移予想を行う予定です。なので必要なデータは、日本の総人口の時系列データです。できるだけ長期間分のデータを集めれるだけ集めたいと思います。

3. データの収集

データの収集方法としては、ネット上に厚生労働省や統計データが大量に存在しているため、それらを利用していこうと思います。以下に利用したネット上のURLを添付しておきます。

4. import内容

上記のURLから、1920年~2000年のデータと2000年~2015年までの日本の人口推移のデータがあるので探してダウンロードします。

今回の全プログラムを通して利用するimport内容は以下に示すとおりです。

import xlrd
import pprint
import numpy as np
import pandas as pd
from scipy import stats
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm
from pylab import rcParams
from datetime import datetime
from statsmodels.tsa.statespace.sarimax import SARIMAX

5. データの取得

データの取得プログラムを以下に示します。今回はcsvファイルではなく、Excelファイルからデータを取得しました。理由は、Exelが世の中では多く使われているためです。

wb = xlrd.open_workbook('男女.xlsx')
sheet = wb.sheet_by_name('第1表')
col_valuesE = sheet.col_values(0)#年代
EF = pd.DataFrame(col_valuesE)
col_valuesB = sheet.col_values(1)#日本の人口
BF = pd.DataFrame(col_valuesB)

PythonでのExcelファイルの利用方法については以下のURLを参照しました。

6. データの整理

今回使用するデータは特にかけている部分もなく1920年~2015年までの95年分のデータがあります。90年分のデータで学習し、2010年 ~2015年までの5年分を機械学習を用いて予測してみようと思います。

自己相関係数

さっそく自己相関係数から調べていきます。

自己相関係数を調べることで、何個前のデータの影響が大きいのかを調べることができます。以下に、自己相関係数を調べるプログラムと結果を示します。

fig = plt.figure(figsize =(12,8))
ax = fig.add_subplot()
sm.graphics.tsa.plot_acf(col_valuesE, lags=80, ax=ax)
rcParams['figure.figsize'] = 20,7
fig = plt.figure(figsize=(10,10),dpi=20)
plt.show()

薄く青い色が付いている範囲が95%信頼区間で、大体7~8個前までのデータに大きく影響を受けていることがこのグラフからわかります。

偏自己相関係数

次に偏自己相関係数を調べます。偏自己相関係数からは、ある点とある点の間にある関係を調べることができます。これを用いることで、一個前のデータの影響を取り除いてデータの関係を調べることができます。以下に、偏自己相関係数を調べるプログラムと、結果を示します。

fig = plt.figure(figsize =(12,8))
ax = fig.add_subplot()
fig = sm.graphics.tsa.plot_pacf(col_valuesB, lags=80, ax=ax)
rcParams['figure.figsize'] = 20,7
fig = plt.figure(figsize=(10,10),dpi=20)
plt.show()

この結果を見て、左から1, 27, 42, 43, 44, 45, 46, 47, 48, 49, 57番目のデータが影響を与えていることがわかります。

特に1, 45, 46, 47のデータが大きな影響を与えており、残りの27, 42, 43, 44, 48, 49, 57のデータは95%信頼区間に近いことから、何らかの影響を与えている可能性があります。

ADF検定

ADF検定では、単位根過程でないかどうかを判定することができます。以下に、ADF検定のプログラムと結果を示します。

adf_result = sm.tsa.stattools.adfuller(col_valuesB, autolag ='AIC')
adf = pd.Series(adf_result[0:4], index = ['Test Statistic', 'p-  value', '#Lags Used', 'Number of Observations Used'])
print(adf)

pの値が、0.671959となっています。p>0.05であるため、単位根過程ということになります。
単位根過程とは、データが「上か下かどちらかの一方に増えていく」ような時系列データのことです。

今回、単位根過程ということで良い学習データとは言えませんが、1920年~2010年までの学習を行い、2010年~2015年までの予測を行わせてみようと思います。

7. SARIMAモデル

モデルの決定と学習

今回使用するデータが時系列データなので、研修時に履修したSARIMAモデルを使用して学習していきます。SARIMAモデルとは、時系列分析の一手法です。複数の時系列モデルを複合した手法と考えることができます。以下に示す図のような関係性があります。

SARIMAモデルについて詳しく示してあるURLを以下に添付します。

学習プログラムに関しては以下に示します。

N=90 #90点までを学習とし、その先5点を予測する
test=BF[:]
learn=BF[:N] #グラフよりs=1,45,46,47のいずれかが良い可能性が高い
SARIMA_BF=sm.tsa.statespace.SARIMAX(learn,order=(1,0,1),seasonal_order = (1,0,0,47), enforce_stationarity = False,
enforce_invertibility = False, trend = "n").fit(trend='nc', disp = False)
print(SARIMA_BF.summary())
pred = SARIMA_BF.predict()
pred2 = SARIMA_BF.forecast(5)

上記のプログラムでは、1920年~2010年までの90点分のデータで学習させ、学習データをもとに2015年前での5年間の人口推移を予想させるプログラムになります。

SARIMA(1,0,1)(1,0,0,47)としています。

モデルの学習結果を以下に示します。

学習結果は、AIC、BICがそれぞれ600を超えています。モデルとしてあまり良くない結果となりました。一般的に、AIC、BICは値が小さいほうが良いモデルとされる傾向があります。

次に、学習結果のデータをグラフにしてみます。グラフ化のプログラムは以下に示します。

plt.title("JAPAN population ")
plt.xlabel("data")
plt.ylabel("population*1000")
plt.ylim([50000,135000])
plt.xlim([0,100])
plt.plot(BF,color="b", label = "Actual value")
plt.plot(pred,color="r", label = "learn value")
plt.plot(pred2,color="y", label = "Predicted value")
plt.legend(loc = 2)
plt.show()

以下に、グラフを示します。

青のグラフが元のデータ、赤のグラフが学習させたデータ、黄色が予測させたデータになります。

黄色の予測データの部分が従来のグラフから外れており、改善が必要となっています。

8. SARIMAモデル(データ標準化後)

データの標準化を行う

先ほどの予測が従来のグラフから外れてしまった一つの要因として、AIC、BICのどちらの値も600を超えるモデルとなっていることから、モデルとして良くなく、学習結果も従来のグラフから外れてしまう結果となってしまったことが考えられます。

そこで、上記の理由より、学習がうまくいかなかった要因の一つとして考えられる、データの標準化を行おうと思います。

標準化については以下に添付するURLを参照してみてください。

追記するプログラムは数行ほどです。

付け加えたプログラムは、上から7行目の標準化のプログラムと、下から3~4行目の出力時の標準化から戻すプログラムのみです。

それでは、さっそくAICとBICから確認していきます。

wb = xlrd.open_workbook('男女.xlsx')
sheet = wb.sheet_by_name('第1表')
col_valuesE = sheet.col_values(0)#年代
EF = pd.DataFrame(col_valuesE)
col_valuesB = sheet.col_values(1)#日本の人口
BF = pd.DataFrame(col_valuesB)
#標準化を行う
BF_mean = np.mean(BF.values) #平均
BF_std = np.std(BF.values) #標準偏差
BF_b = (BF - BF_mean)/BF_std #標準化
N=90#90点までを学習とし、その先5点を予測する
test=BF[:]
learn=BF_b[:N]#グラフよりs=1,45,46,47のいずれかが良い可能性が高い
SARIMA_BF=sm.tsa.statespace.SARIMAX(learn,order=(1,0,1),seasonal_order = (1,0,0,47), enforce_stationarity = False,
enforce_invertibility = False, trend = "n").fit(trend='nc', disp = False)
print(SARIMA_BF.summary())
pred = SARIMA_BF.predict()
pred2 = SARIMA_BF.forecast(5)
plt.title("JAPAN population ")
plt.xlabel("data")
plt.ylabel("population*1000")
plt.ylim([50000,135000])
plt.xlim([0,100])
del pred[0]
plt.plot(BF,color="b", label = "Actual value")
#標準化していたので出力時に元の形に戻す。
plt.plot(pred * BF_std + BF_mean,color="r", label = "learn value")
plt.plot(pred2 * BF_std + BF_mean,color="y", label = "Predicted value")
plt.legend(loc = 2)
plt.show()

下に示すのが、標準化後のデータから学習した際のモデルの結果です。AICとBICの値が、先ほどの600から-200まで下がっています。

この時、負号が反転してしまっていますが、AICとBICの値は、0に近いほど良いモデルとなっていることに変わりはありませんので、モデルとして先ほどのモデルより、良いモデルとなっていると考えられます。

そして、こちらが標準化したデータを用いた予測結果です。

青が実際値グラフ、赤がデータ学習結果のグラフ、黄色が予測結果のグラフになります。 このままでは、予測精度が向上しているかわかりにくいので、下に、元のモデルとの比較のグラフを添付します。

こちらが、予測結果の比較のグラフになります。

青が実際値のグラフ、赤がデータ処理なしの予測結果のグラフ、黄色が標準化後の予測結果のグラフになります。

データの標準化を行った結果により、データ処理を行わないときに比べ少しだけ従来のグラフに近づいています。

しかし、このままでは、人口推移の予測など、とてもじゃありませんができそうにありません。

そこで、次から新しいモデルを用いて日本の人口推移を予測しようと思います。

9. 状態空間モデル

次に、状態空間モデルを用いた日本の人口推移を予想していこうと思います。

状態空間モデルとは、多くの統計モデルを統一的に表すことが可能な統計モデルになります。 今回、状態空間モデルを利用する理由としては、 分散分析や回帰分析、正規線形モデルだけでなく、一般化線形モデル、ARIMAなどを状態空間モデルというたった一つのモデルで表すことが可能であるからです。このモデルであれば、先ほどの日本人口推移を予想することができるかもしれません!

状態空間モデルについて簡単に説明がしてあるURLを以下に添付します。

では、さっそく、やっていこうと思います。

SARIMAモデルのプログラムに使用したimport内容を利用します。

下に先ほどと同様にExcelファイルを読み込んだ後に、90点までの学習データを、分けておきます。

wb = xlrd.open_workbook('男女.xlsx')
sheet = wb.sheet_by_name('第1表')
col_valuesE = sheet.col_values(0)#年代
EF = pd.DataFrame(col_valuesE)
col_valuesB = sheet.col_values(1)#日本の人口
BF = pd.DataFrame(col_valuesB)
N = 90 #90点まで学習させる。
BF_t = BF[:N]

今から、さっそく学習を進めていこうと思います。今回は、2つのモデルを用いて、精度の良い方を用いて予想を行うこととします。

ローカルレベルモデル

ローカルレベルモデル とは、実際のデータに単にノイズを追加したデータを加えることで、次の状態を推定するモデルになります。

簡単なモデルですが、簡単な割に精度は高い方だと思います。以下に、そのプログラムを添付します。

mod_local_level = sm.tsa.UnobservedComponents(BF_t, 'local level')
res_local_level = mod_local_level.fit()

ローカル線形トレンドモデル

ローカル線形トレンドモデルとは、先ほどのローカルレベルモデルに傾きの要素が加わったものになります。日本の人口推移のグラフを見た時、右肩上がりのグラフとなっています。つまり、右肩上がりの傾きがあるということがわかります。

これより、ローカル線形トレンドモデル は、今回の日本の人口推移の予想にはうってつけのモデルになります。日本の人口推移は右肩上がり一方のデータになっており、上昇トレンドであるため、このモデルを用いることでモデル化がうまくいくと予想されます。こちらも以下にそのプログラムを添付します。

mod_trend = sm.tsa.UnobservedComponents(BF_t,'local linear trend')
res_trend = mod_trend.fit(method='bfgs')

今回、2つのモデルを利用します。その際、2つのモデルを互いに比較し、AICの値が良い方を利用することとします。AICは値が小さい方が精度が良いので、値の小さい方を使用することとします。下にプログラムを添付します。

AIC_list = pd.DataFrame(index=range(2), columns=["MODEL", "AIC"])
AIC_list.ix[0]["MODEL"] = "res_local_level"
AIC_list.ix[0]["AIC"] = res_local_level.aic
AIC_list.ix[1]["MODEL"] = "res_trend"
AIC_list.ix[1]["AIC"] = res_trend.aic
print(AIC_list)

こちらが、プログラムの実行結果となります。

MODEL AIC
0   res_local_level 1492.15
1   res_trend 1399.25

上がローカルレベルモデルで、下がローカル線形トレンドモデルになります。ローカル線形トレンドモデルの方が、AICの値が100程低いので、ローカル線形トレンドモデルを使用することとします。

下に示すのが、日本の人口推移の予測を行うプログラム部分となります。学習させた後の値の最初に0が格納されてしまうため、最初の要素だけを2行目で削除しています。

3行目で、2010年から2015年までの予測を行っています。

pred_t = res_trend.predict()
del pred_t[0]
pred2 = res_trend.forecast(5)

下には、先ほどとほとんど同じグラフ表示用のプログラムを添付します。

rcParams['figure.figsize'] = 8, 6
plt.title("JAPAN population ")
plt.xlabel("data")
plt.ylabel("population*1000")
plt.plot(BF,"b",label = "Actual value")
plt.plot(pred_t, "r", label = "learn value")
plt.plot(pred2,"y",label = "Predicted value")
plt.legend(loc = 2)
plt.show()

そして、こちらが学習結果をグラフ表示させたものになります。

青のグラフが、もともとのデータを表しています。赤のグラフが、学習後のデータを表しています。黄色のグラフが、予測結果を表しています。

最初の方こそ、少しの誤差が生じていますが、その後はほとんど元のデータのグラフと同じ曲線を描いています。そのことから、SARIMAモデルを用いたときに比べ学習は各段に向上していることがわかります。

また、黄色のグラフを見てください。先ほどのSARIMAモデルを用いた予測結果と比べても格段に、正規のデータに近いグラフとなっています。

同じグラフ上にこの部分を拡大したグラフを表示させようと思います。

下に示すのが、表示範囲を拡大したグラフになります。

青が元のデータ、黄色がローカル線形トレンドモデル、緑がSARIMAモデルのグラフになります。

SARIMAモデルに比べ、ローカル線形トレンドモデル精度が格段に良いことがわかります。

最終的に、かなり良い精度の予測ができたのではないでしょうか。これらに関しての考察を行っていこうと思います。

10. 考察

今回、日本の人口推移をSARIMAモデルを用いて予想しようとしたが、実際のグラフから大きく外れる結果となりました。日本に限らず、人口推移の時系列データはGDPなどと同様に、単位根過程のものであり、今回のデータは特に出力に大きな影響を与えたのだと考えられます。

それらの要因を踏まえたうえで、ローカル線形トレンドモデルを利用したところ、SARIMAモデルに比べ格段に精度が向上したことから、それぞれのモデルの特徴から原因を考えようと思います。

まず、SARIMAモデルはブログの途中に添付した図からもわかるように、ARIMAモデルの拡張として、季節成分とトレンド性を加えたものになります。しかし、SARIMAモデルのトレンド性は、グラフが上下に動きながら徐々に上昇していくようなトレンド性に対しての学習は強いですが、今回のようにただ上昇していくデータの場合、トレンドの傾きがかなり大きくなってしまいます。そのため、予測結果がそのトレンド成分の影響を受けることで、実際の値より大きい予測結果が出力されたのだと考えられます。

また、SARIMAモデルのトレンド成分は、学習データの最初のデータから最後のデータまでのデータから得たトレンド成分である可能性が考えられます。その理由は、本来のデータがdata 60 ~ data 80の段階で曲線が緩やかになっているため、トレンド成分も緩やかになるべきであるのに対して、予測結果は、data 0 ~ data 90までの直線の延長線上にグラフが出力されていることから推測されます。

それに比べ、ローカル線形トレンドモデルの場合は、時間(データ)ごとに切片とトレンド成分の変化をさせていくことができるモデルであるため、data 60 ~ data 80までの曲線が緩やかになる部分のデータから予測データが予測されていると考えられます。

今回の日本の人口推移のモデルにはローカル線形トレンドモデルが向いていると考えられます。また、ローカル線形トレンドモデルは、人口だけでなく、緩やかなトレンド性をもつデータの予測に対して大きな力を発揮してくれると考えられます。

11. まとめ

今回、日本に人口推移のみについて機械学習を用いたが社会背景と照らし合わせて考えると、現在の日本は少子高齢化社会を迎える段階に至っており、人口推移は昔に比べ緩やかな曲線となっています。そのため、今回のように昔のデータで学習させたSARIMAモデルで現在の人口推移を予想させると、予想データのグラフが従来データのグラフの値を上回ってしまう結果となります。

しかし、人口は未だ減少傾向にあるデータがない(地方の人口データは例外)ため、SARIMAモデルが向いていないと完全に断言することはできないです。

ただ、緩やかなトレンド性を持ち、尚且つ前のデータの影響が大きく影響するようなデータの場合は、ローカル線形トレンドモデルの使用がお勧めできます。

1 COMMENT

Amy

I could not resist commenting. Perfectly written! I could not resist commenting.
Very well written! I’ve been browsing online more than 4 hours today, yet I never found any interesting article like yours.
It is pretty worth enough for me. In my opinion, if all website owners
and bloggers made good content as you did, the internet will be
much more useful than ever before. http://Nestle.com/

返信する

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です