scikit-learnで将来の日本の人口を予測する。[python][非線形][機械学習]

python

はじめに

今回はscikit-learn(サイキット・ラーン)を用いて、時系列データの機械学習を行う。
scikit-learnは機械学習のライブラリで機械学習の様々なアルゴリズムを持っている。

会社の研修でデータ分析の研修を実施したため、自分でpythonを用いて何か分析したいなと思い、
ネットに無料で落ちているデータ(※)を用いて、機械学習を行うことにした。
※全国の市町村の人口推移データを今回は利用した。

冒頭の写真のように(ネット上にある)様々なデータに対して非線形回帰曲線のモデリングが可能であるため、最後までぜひ見てください。

利用するデータ

データは何でも良かったので、ネットで「無料 データ 取得」とググると様々なサイトがヒットした。

今回はその中でRESAS(リサース)というサイトからデータを取得することにした。
以下、RESASサイトからの抜粋。

地域経済分析システム(RESAS:リーサス)は、地方創生の様々な取り組みを情報面から支援するために、経済産業省と内閣官房(まち・ひと・しごと創生本部事務局)が提供しています。

自治体職員の方や、地域の活性化に関心を持つ様々な分野の方によって、効果的な施策の立案・実行・検証のためなどに広く利用されています。
お役に立つご参考情報もこちらに掲載しております。

時系列データは昔から分析してみたいと思っていたので、ここのサイトにある全国の市町村人口推移データを用いて、将来の人口分布を予測していく。

データの取得方法

下記のようにメニューより、「人口マップ > 人口構成」を選択して下さい。

すると、下記のように人口分布のマップが出力されますので、「人口推移」ボタンをクリック。

すると、下記の画面に飛びますので、データをダウンロードして下さい。

データの解凍方法

Windowsの方

そのまま、解凍して下さい。

Macの方

元データの文字コードがShift-JISとなっており、単純に解凍すると文字化けします。
アプリ等のインストール不要の手順を紹介しますので、下記の手順で解凍して下さい。

macで文字化けせずに、WindowsのZIPファイルを回答する方法。(アプリのインストール不要)

線形で機械学習を行う

さて、データをダウンロードした方は実際にscikit-learnを利用して、線形回帰曲線を求める。
以下がコードの全量となる。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression

#ローカル保存したcsvデータ
PATH =  "../../../unzip-folder/population-composition_20200827/02_人口推移/人口_人口構成_人口推移_市区町村.csv"

#SHIFT-JIS形式で読み込む。
df = pd.read_csv(PATH,encoding="SHIFT-JIS")
#列名を定義
df = df.rename(columns={"集計年":"year","都道府県コード":"t_code","都道府県名":"t_name","市区町村コード":"city_code","市区町村名":"city_name","総人口(人)":"total","年少人口(人)":"low","生産年齢人口(人)":"middle","老年人口(人)":"high"})
#2020年より大きいのデータを削除。
df = df[(df["year"] <= 2020)]

#市町村のデータでサンプル数が少ないものは削る。
#市町村データのMAXは9。9未満のデータは削除したい。
#9未満のデータのcity_codeを取得する。
vc = df['city_code'].value_counts()
vals_to_remove = vc[vc < 9].index.values
vals_to_remove.sort()

#9未満の市町村データを持つ行のインデックスを格納するリスト。
IDX = []

for index,item in df.iterrows():
    if(item.city_code in vals_to_remove):
        # vals_to_remove リストにitem.city_codeが含まれていれば、IDXにitem.city_codeを格納する。
        IDX.append(index)

#9未満の市町村データを削除。drop関数を利用。
df_drop = df.drop(index = IDX)

#都道府県(t_code)と市町村(city_code)の重複なしリストの作成。
t_code = df_drop['t_code'].unique()
city_code = df_drop['city_code'].unique()

# date_list =  [[t_code],[city_code],[coef_],[intercept_]]
date_list = [[],[],[],[]]

def create_fit(df):
    '''sklearnでfit関数を用いて切片を導出'''
    lr = LinearRegression()
    X = df[['year']].values
    Y = df[['total']].values
    lr.fit(X,Y)
    
    return [lr.coef_[0][0] , lr.intercept_[0]]

#data_listに値を代入する。
for city in city_code:
    
    df = df_drop[(df_drop.loc[:,'city_code'] == city)]
    
    date_list[0].append(df['t_code'].unique()[0])
    date_list[1].append(city)
    a,b= create_fit(df)
    date_list[2].append(int(a))
    date_list[3].append(int(b))

#データフレームのindexを指定。
N = 1

def plot_date(N):
    '''生データとモデリングした関数を比較。'''
    X_FIT = np.array(range(1980,2025,5))
    Y_FIT = X_FIT * date_list[2][N] + date_list[3][N]
    plt.plot(X_FIT,Y_FIT,'--')

    #df_drop
    df_scop = df_drop[(df_drop['t_code'] == date_list[0][N]) & (df_drop['city_code'] == date_list[1][N])]
    X = df_scop[['year']].values
    Y = df_scop[['total']].values
    plt.plot(X,Y,'o')

    plt.show()
    
    return df_scop

#生データとモデルをplotする。
plot_date(N)

 

結果は以下となる。

北海道の札幌市中央区の人口推移を線形回帰曲線を用いて導出した。
なんとなく直線上に生データはあるが、ズレている部分もある。
※非線形回帰を用いると、ほぼ一致する。非線形回帰はこのあと説明するので、少々お待ちを。

ズレているのは、一旦置いておいて、コードの重要な部分を解説しよう。

2020年より大きいのデータを削除。

以下の部分の解説。

df = df[(df["year"] <= 2020)]

このコードを実施する前のdfは以下の形になっている。

そう、ここのサイトは2020年より大きいデータも予測しているのだ。
ただ、今回は過去のデータを元にモデリングしたいため、2020年より大きいデータを削除したい。
下記のように”year”列で2020以下の部分はTrue、2020より大きい値はFalseが返ってくる。

df["year"] <= 2020

さらに、下記のようにdf[…]で囲むと、Trueの列のみ抽出される。
結果、2020年以下のデータフレームを取得できる。

市町村のデータでサンプル数が少ないものは削る。

次にcity_codeを要素数を確認する。
value_counts()を利用すると、ある列の要素数を確認できる。

city_codeを見ると、1個のものもある。
これは市町村合併などでデータが欠損していることを表す。
市町村のデータ個数は揃えた方が学習しやすいので、今回は9個未満のデータは全て削除する。

#city_codeの要素数をカウント。
vc = df['city_code'].value_counts()

#city_codeの要素数が9未満のindexを取得。
vals_to_remove = vc[vc < 9].index.values

#vals_to_removeを昇順にソート。
vals_to_remove.sort()

#vals_to_removeを表示。
vals_to_remove

これがdfの中にあるcity_codeが9未満のindexの配列である。
これを元のdfから削除(drop)すると、city_codeが9未満のdfを取得することができる。

最終的に16074行となり、前節のデータフレームの16818行より少なくなっている事がわかる。

sklearnでfit関数を用いて傾きと切片を導出

今回はある市町村の線形回帰を求めるため、大量のfit関数を実行する。
何度も実行するので、その線形回帰の傾きとy切片を求めるcreate_fit関数を作成する。
その時の引数は指定の市町村が入っているデータフレーム型とする。

その時の関数が以下である。

def create_fit(df):
    '''sklearnでfit関数を用いて傾きと切片を導出'''
    #インスタンスを生成。
    lr = LinearRegression()
    #X,Yの列を作成。
    X = df[['year']].values
    Y = df[['total']].values
    #fit関数で線形回帰を作成。
    #fit関数を実行後、lrにconf_(傾き)とintercept_(切片)を求める事ができる。
    lr.fit(X,Y)
    
    #戻り値ではrのconf_(傾き)とintercept_(切片)を設定。
    return [lr.coef_[0][0] , lr.intercept_[0]]

生データとモデリングした関数を描写

最後にplot_data関数を用いて、モデリングした線形回帰のグラフと生データを描写したい。
引数はデータフレームのインデックスをもつ。

plotのオプションに’–‘をつけると、点線でグラフを描写し、’o’とすると散布図を描写できる。
最後にdf_scopを返しているのは、描写した都道府県と市町村を知るためだ。

def plot_date(N):
    '''生データとモデリングした関数を描写。Nはデータフレームのインデックスを表す。'''
    #傾きとy切片をfit関数で作成したモデル線形回帰のグラムを描写。
    X_FIT = np.array(range(1980,2025,5))
    Y_FIT = X_FIT * date_list[2][N] + date_list[3][N]
    plt.plot(X_FIT,Y_FIT,'--')

    #df_scopはある都道府県かつ市町村を持っているデータフレーム。
    #その生データをscatter plotする。
    df_scop = df_drop[(df_drop['t_code'] == date_list[0][N]) & (df_drop['city_code'] == date_list[1][N])]
    X = df_scop[['year']].values
    Y = df_scop[['total']].values
    plt.plot(X,Y,'o')

    plt.show()
    
    #戻り値にdf_scopを返す。
    return df_scop

 

非線形で機械学習を行う(2乗以上)

線形回帰だと直線となり、モデルと生データが合致しないことが多い。
sklearnには非線形回帰もfit関数で実装できるため、非線形回帰の曲線も作成してみる。

ソースは以下。

from sklearn.kernel_ridge import KernelRidge
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

#ローカル保存したcsvデータ
PATH =  "../../../unzip-folder/population-composition_20200827/02_人口推移/人口_人口構成_人口推移_市区町村.csv"

#SHIFT-JIS形式で読み込む。
df = pd.read_csv(PATH,encoding="SHIFT-JIS")
#列名を定義
df = df.rename(columns={"集計年":"year","都道府県コード":"t_code","都道府県名":"t_name","市区町村コード":"city_code","市区町村名":"city_name","総人口(人)":"total","年少人口(人)":"low","生産年齢人口(人)":"middle","老年人口(人)":"high"})
#2020年より大きいのデータを削除。
df = df[(df["year"] <= 2020)]

#市町村のデータでサンプル数が少ないものは削る。
#市町村データのMAXは9。9未満のデータは削除したい。
#9未満のデータのcity_codeを取得する。
vc = df['city_code'].value_counts()
vals_to_remove = vc[vc < 9].index.values
vals_to_remove.sort()

#9未満の市町村データを持つ行のインデックスを格納するリスト。
IDX = []

for index,item in df.iterrows():
    if(item.city_code in vals_to_remove):
        # vals_to_remove リストにitem.city_codeが含まれていれば、IDXにitem.city_codeを格納する。
        IDX.append(index)

#9未満の市町村データを削除。drop関数を利用。
df_drop = df.drop(index = IDX)

#都道府県(t_code)と市町村(city_code)の重複なしリストの作成。
t_code = df_drop['t_code'].unique()
city_code = df_drop['city_code'].unique()
# date_list =  [[t_code],[city_code],[Y]]
date_list = [[],[],[]]

def create_fit(df):
    
    clr = KernelRidge(alpha=0.1, degree=3, kernel='poly')
    X = df[['year']].values - df[['year']].values.min()
    Y = df[['total']].values - df[['total']].values.min()
    clr.fit(X,Y)
    
    return clr.predict(X) + df[['total']].values.min()
    
for city in city_code:
    
    df = df_drop[(df_drop.loc[:,'city_code'] == city)]
    
    date_list[0].append(df['t_code'].unique()[0])
    date_list[1].append(city)
    date_list[2].append(create_fit(df))

#print(date_list)
import matplotlib.pyplot as plt

N =1

X = df_scop[['year']].values
plt.plot(X,date_list[2][N],'--')

df_drop
df_scop = df_drop[(df_drop['t_code'] == date_list[0][N]) & (df_drop['city_code'] == date_list[1][N])]
Y = df_scop[['total']].values
plt.plot(X,Y,'o')
plt.show()

df_scop

決定係数と呼ばれるモデリングの関数がどれだけ生データと近い値かを示す指標がある。
1に近ければ、近いほど生データに合致しているこの指標は、
線形回帰の時は0.3位だったが、非線形回帰モデルの時は0.8-0.9が平均となっていた。

よって、今回の人口推移モデルは非線形回帰の方がより良いモデルである事がわかる。

KernelRidgeのライブラリを用いることで、非線形回帰の曲線を作成できる。
以下、線形回帰と差分となっている箇所を重点的に解説する。

非線形回帰曲線を描写

線形回帰と同じく、非線形回帰を求めたい。
線形回帰はreturn文に係数と切片を設定したが、
非線形回帰の場合はYの値のリストを戻り値とする。

KernelRidgeのインスタンス作成時のkernelは関数の種類を表ており、
線形・多項・シグモイド・ガウス・RBFなど様々ある。
今回は非線形回帰でn乗の和の形を表ているため、kernel=’poly’としている。
また、次数(degree)は3次の多項式を表す。

def create_fit(df):
    
    clr = KernelRidge(alpha=0.1, degree=3, kernel='poly')
    X = df[['year']].values - df[['year']].values.min()
    Y = df[['total']].values - df[['total']].values.min()
    clr.fit(X,Y)
    
    return clr.predict(X) + df[['total']].values.min()

ポイントはX,Yを求める時に自身の列の最小値を引いていることだ。
こうすることで生データが原点近くに推移し、KernelRidgeでfit関数を実行がしやすくなる。(※)
※正規化・標準化の類。この工夫をするだけで、結果が全く違う。

戻り値はpredict関数を使って、モデル曲線のYの値のリストを返す。
その時に自身の列の最小値を引いたため、足し直すことを忘れずに!

以降は線形代数の同様のやり方でグラフ描写でいるため、割愛する。

最後に

以上でscikit-learnを用いたpythonで線形・非線形のモデリングの作成は終了である。
今回はN(インデックス)を変更することで、各都道府県・市町村の人口推移モデルを作成できる。
また、このモデルを用いることで、将来の人口推移も予測することが可能である。

データのサンプル数が9つと少なかったため、若干過学習になりがちだが、
そこはデータの前処理をうまく実施することで、回避も可能だ。
※例えば、市町村のモデルに都道府県の人口推移のパラメータを加味するとか。

pythonの良いところは、このように機会学習のライブラリが多くあるため、
これらを活用することで、初心者でも簡単にデータ分析できることが強みだ。

このブログで紹介したLinearRegressionとKernelRidgeを利用すれば、線形・非線形モデルは楽勝だと思う。(fit関数を実行するだけなので)

他にも決定木やk-means法など様々な機械学習があるので、どんどん紹介していきたい。

質問などありましたら、コメント下さい。

以上です。ありがとうございました。

コメント

タイトルとURLをコピーしました