Examples¶

世界の年平均気温の変化¶

1991年から2020年の平均からの偏差

気象庁のページ https://www.data.jma.go.jp/cpdinfo/temp/list/an_wld.html より取得したもの

In [52]:
import pandas as pd
url = "https://toyoki-lab.ee.yamanashi.ac.jp/~toyoki/lectures/PracDataSci/data/annual_temp_world.csv"
df = pd.read_csv(url)
df.head()
Out[52]:
年 世界全体 北半球 南半球
0 1891 -0.78 -0.88 -0.68
1 1892 -0.89 -1.00 -0.74
2 1893 -0.94 -1.06 -0.79
3 1894 -0.86 -0.93 -0.77
4 1895 -0.82 -0.95 -0.67
In [53]:
import matplotlib
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'IPAPGothic'  # 全体のフォントを設定
df.plot(x="年")
Out[53]:
<Axes: xlabel='年'>

前処理：機械学習用にデータ加工¶

  • 説明変数(年)の原点を最初の年に移動。べき近似のため

  • べき関数を基底関数にとることとして、説明変数のデータ（計画行列）を作成

    • カーネル法(SVR, RVRでは必要ない)

  • 正規化 (経験上、ニューラルネットモデルの場合には必要）

In [54]:
import numpy as np

# べき関数近似を行うために説明変数(年)の原点移動をして学習用データを作成
first_year = 1891

# 目的変数
target_col = "世界全体"

# M次多項式近似の場合は、べき乗の値を列とするM+1列の行列(計画行列)
deg = 5
poly_feature = np.vander((df['年'] - first_year).to_list(), deg+1)

# SVR, RVRはカーネルを指定するだけなので、1列だけでよい
single_feature = np.array([(df['年'] - first_year).to_list()]).T

# 説明変数の正規化 (Neural Netについては、正規化データを与えないとよい結果が得られない
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler().fit(poly_feature) #正規化
scaled_poly_feature = scaler.transform(poly_feature)

#single_feature
target = df[target_col].to_list()

# テストデータ
n=100
test_single = np.linspace(0, 2022-first_year, n)
test_single = test_single[:, np.newaxis]
test_poly = np.vander(np.linspace(0, 2022 - first_year, n)[:], deg+1)
scaled_test_poly = scaler.transform(test_poly)

モデルの設定・学習・テスト¶

In [109]:
# scikit-learnの各種モデル
# {model, title(name), 説明変数, 目的変数, test(predict)用データ}の配列を用意する

models = {}

# 非線形回帰:多項式近似 (vander行列をデータとして入力 ⇒ 線形重回帰へ入力)
from sklearn import linear_model as lm
models['lm'] = {"model": lm.LinearRegression(),
               "name": "Simple Polynomial Model" + str(deg) +"次",
              "feature": poly_feature, "target": target, "test_data": test_poly}

# SVR
from sklearn import svm
models['svr'] = {"model": svm.SVR(kernel='rbf'),
               "name": "Support Vector Regression",
              "feature": single_feature, "target": target, "test_data": test_single}

# RVR (pipでインストールする必要あり)
from sklearn_rvm import EMRVR
models['rvr'] = {"model": EMRVR(kernel='rbf', gamma="scale"),
               "name": "Relevant Vector Regression",
              "feature": single_feature, "target": target, "test_data": test_single}

# Neural Network (よいパラメータが見つけられなかった)
from sklearn.neural_network import MLPRegressor as MLPR
models['nn'] = {"model": MLPR((100,), activation="identity", solver="lbfgs",
                            max_iter=2000, tol=1e-8),
               "name": "Neural Network",
               "feature": scaled_poly_feature, "target": target, "test_data": scaled_test_poly}


# Random Forest
from sklearn.ensemble import RandomForestRegressor as RFR
models['rfr'] = {"model": RFR(max_depth=3, n_estimators=100), # 深さをいろいろ変えてみる
               "name": "Random Forest Regression",
               "feature": poly_feature, "target": target, "test_data": test_poly}

# Mixture model

結果のプロット¶

In [110]:
import matplotlib.pyplot as plt
import matplotlib
plt.rcParams['font.family'] = 'IPAPGothic'  # 全体のフォントを設定

fig, ax = plt.subplots()
for m in models.values():
    m['model'].fit(m['feature'], m['target'])
    result = m['model'].predict(m['test_data'])
    ax.plot(test_single+first_year, result, label=m['name'])
    
ax.legend()
df.plot(kind='scatter',x="年", y=target_col, ax=ax)
ax.set_ylabel("気温の偏差 (基準：1991年～2020年の平均)")
Out[110]:
Text(0, 0.5, '気温の偏差 (基準：1991年～2020年の平均)')

線形重回帰とニューラルネットはほとんど重なっている。

表示するモデルをインタラクティブに選択するためのipywidgetsを使った例を下に示す。（チェックボックスのほうがわかりやすいのだが、うまい方法が見つからなかった。）

複数選択するときは、ctrlキーを押しながら選択すればよい。（html版では、表示されない。）

In [126]:
import ipywidgets

def plot_results(plot_m):
    fig, ax = plt.subplots()
    for k in plot_m:
        m = models[k]
        m['model'].fit(m['feature'], m['target'])
        result = m['model'].predict(m['test_data'])
        ax.plot(test_single+first_year, result, label=m['name'])
    ax.legend()
    df.plot(kind='scatter',x="年", y=target_col, ax=ax)
    ax.set_ylabel("気温の偏差 (基準：1991年～2020年の平均)")
        
select_plots = ipywidgets.SelectMultiple(
    options=models.keys(),
    value=["lm", "svr", "rfr"],
    #rows=10,
    description='表示モデル',
    disabled=False
)

ipywidgets.interact(plot_results, plot_m = select_plots)

interactive(children=(SelectMultiple(description='表示モデル', index=(0, 1, 4), options=('lm', 'svr', 'rvr', 'nn', …
Out[126]:
<function __main__.plot_results(plot_m)>

Relevance Vector Regression(RVR)では誤差範囲が得られる

In [57]:
# 学習を行う
m = models['rvr']
m['model'].fit(m['feature'], m['target'])
  
# テスト
y_test, y_std = m['model'].predict(m['test_data'], return_std=True) # standard dev. can be also obtained if "return_std" is set to be True
    
# plot prediction and the the width of std. dev.
fig, ax = plt.subplots()
ax.plot(m['test_data']+first_year, y_test)
ax.fill_between(list(m['test_data'][:,0]+first_year),
                y_test - y_std, y_test + y_std, color="darkorange", alpha=0.2) #xデータは1次元リストに変換
# 生データのプロット
df.plot(kind='scatter',x="年", y=target_col, ax=ax)
ax.set_ylabel("気温の偏差 (基準：1991年～2020年の平均)")
Out[57]:
Text(0, 0.5, '気温の偏差 (基準：1991年～2020年の平均)')

Mixture Linea Model

In [112]:
%reload_ext autoreload
import mixture_lm as M_LM
feature = M_LM.PolynomialFeatures(degree=1)
X_train = feature.transform(single_feature)
In [113]:
# 条件付き混合モデルの用意（モデル数:n_models) 
n_models =2
y_scaler = 100.0 # targetの値を拡大する(ガウシアンの計算エラー回避)

model = M_LM.ConditionalMixtureModel(n_models=n_models)
# 学習（パラメータの最尤推定）
model.fit(X_train, np.array(target)*100,n_iter=1000)
# 推定
Y = model.predict(X_train)

[-94.34951014   0.73908574]
In [108]:
import math

# Figureの初期化
fig = plt.figure(figsize=(7, 5))
# Figure内にAxesを追加()
ax = fig.add_subplot(111) #...2

for i, c in enumerate(model.getResponsibility()): # データ点ごとに負担率(responsibility)で色付け描画
    if len(c) == 2:
        color = [c[0], 0, c[1]]
    else:
        color = c
    ax.scatter([df['年'][i]],[df[target_col][i]], color=color)

# 回帰線の描画
ax.plot(single_feature+first_year, Y[:, 0]/y_scaler, "r-")
ax.plot(single_feature+first_year, Y[:, 1]/y_scaler, "b-")

#ax.plot(x, Y[:, 2])
std_dev = math.sqrt(model.var)
# 標準偏差の幅で色付けしてみる
#ax.fill([x[0],x[0],x[-1],x[-1]],
#        [Y[0,0]-std_dev, Y[0,0]+std_dev, Y[-1,0]+std_dev, Y[-1,0]-std_dev], alpha=0.2)
#ax.fill([x[0],x[0],x[-1],x[-1]],
#        [Y[0,1]-std_dev, Y[0,1]+std_dev, Y[-1,1]+std_dev, Y[-1,1]-std_dev], alpha=0.2)
#ax.fill([x[0],x[0],x[-1],x[-1]],
#        [Y[0,2]-std_dev, Y[0,2]+std_dev, Y[-1,2]+std_dev, Y[-1,2]-std_dev], alpha=0.2)
plt.show()
print("傾き(℃/100年) %f  %f" % (model.coef[1][0], model.coef[1][1]))
#print(model.coef/y_scaler) # 回帰係数
#print(model.var/y_scaler)  # 分散

傾き(℃/100年) 0.350644  0.880092
In [51]:
%%html
<link rel="stylesheet" type="text/css" href="custom.css">
In [ ]: