気温と湿度による降雪現象の分類

気温と湿度による降雪現象の分類

北海道室蘭市の降雪現象の分類

北海道室蘭市の気象データ

 2017年11月1日から2018年1月31日にかけての、北海道室蘭市の最低気温 $\boldsymbol{\mathrm{X}}$ と平均湿度 $\boldsymbol{\mathrm{Y}}$ , および「降雪現象なし」を表す $\boldsymbol{\mathrm{T}}$ のデータを用意しました(気象庁提供)。「降雪現象なし」は、「降雪があった場合」を 0, 「降雪がなかった場合」を 1 で表します。言い換えると、雨の中に少しでも雪が混じっていれば 0, 1日中雨だったら 1 となります。

# リストM7-A-1

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from scipy.optimize import minimize

# 2017年11月1日~2018年1月31日
# 北海道室蘭市における気象観測データ

# 最低気温[℃]
X = np.array([-10.1, -9.9, -7.8, -7.4, -6.9, -6.8, -6.5, -6.4, -5.5,
              -5.4, -5.2, -5.0, -4.9, -4.9, -4.7, -4.7, -4.5, -4.5,
              -4.3, -4.3, -4.3, -4.1, -4.1, -4.1, -4.0, -3.8, -3.8,
              -3.6, -3.5, -3.4, -3.4, -3.3, -3.2, -3.1, -3.1, -3.0,
              -3.0, -2.9, -2.9, -2.9, -2.6, -2.5, -2.4, -2.3, -2.3,
              -2.2, -2.1, -2.1, -2.1, -2.0, -1.8, -1.6, -1.6, -1.5,
              -1.5, -1.4, -1.3, -1.0, -0.8, -0.4, -0.3, -0.2, -0.2,
              -0.1,  0.0,  0.3,  0.6,  0.7,  0.7,  0.7,  0.9,  1.0,
               1.1,  1.3,  1.7,  2.0,  2.5,  2.8,  3.6,  3.9,  4.0,
               4.1,  4.3,  4.4,  4.7,  7.4,  7.9,  9.7,  9.8,  11.1])

# 平均湿度[%]
Y = np.array([74, 67, 60, 63, 65, 61, 72, 70, 73,
              74, 69, 74, 67, 81, 72, 68, 74, 66,
              73, 70, 71, 65, 58, 71, 76, 66, 81,
              67, 85, 62, 70, 67, 70, 69, 75, 74,
              71, 79, 79, 69, 68, 77, 75, 79, 68,
              72, 72, 72, 75, 58, 83, 64, 73, 74,
              68, 67, 78, 67, 69, 79, 84, 67, 72,
              69, 77, 78, 69, 69, 74, 66, 73, 74,
              82, 61, 74, 70, 71, 81, 89, 69, 63,
              63, 63, 77, 58, 72, 87, 83, 80, 75])

Y2 = Y / 100

# 降雪現象なし
T = np.array([0, 0, 1, 1, 0, 1, 1, 0, 1,
              0, 0, 0, 1, 0, 0, 1, 0, 1,
              0, 0, 1, 0, 1, 1, 0, 1, 1,
              1, 0, 1, 1, 1, 0, 0, 0, 0,
              0, 0, 0, 1, 1, 1, 0, 0, 0,
              1, 0, 1, 0, 1, 0, 1, 1, 0,
              0, 0, 1, 0, 1, 0, 1, 1, 1,
              1, 1, 1, 1, 1, 1, 1, 1, 1,
              1, 1, 1, 1, 1, 1, 1, 1, 1,
              1, 1, 1, 1, 1, 1, 1, 1, 1])

 Y を 100 で割ったデータ Y2 は、scipy.optimize() でパラメータサーチを実行するときに使用するものです(百分率で入力すると少し厄介な問題が生じるので普通の比率に直してあります)。今回は気温と湿度の 2 種類のデータを入力して「降雪現象あり」と「降雪現象なし」に分類することを試みます。実際の観測データがどのような分布になっているかを Matplotlib を使って図示してみましょう。

# リストM7-A-1-1

# 降雪現象が観測された日の最低気温と平均湿度
X_0 = X[T == 0]
Y_0 = Y[T == 0]

# 降雪現象が観測されなかった日の最低気温と平均湿度
X_1 = X[T == 1]
Y_1 = Y[T == 1]

# FigureとAxesの設定
fig, ax = plt.subplots(figsize = (8, 5))
ax.grid()
ax.set_xlim(-10, 10)
ax.set_ylim(50, 90)
ax.set_xlabel("Temperature [${}^\circ$C]", fontsize = 14)
ax.set_ylabel("Daily mean humidity [%]", fontsize = 14)

# 降雪現象が観測された日のデータを青色で表示
ax.scatter(X_0, Y_0, color = "blue",
           s = 16, label = "snow", zorder = 2)

# 降雪現象が観測されなかった日のデータを赤色で表示
ax.scatter(X_1, Y_1, color = "red",
           s = 16, label = "rain", zorder = 2)

ax.legend()
plt.show()

北海道室蘭市の気温・湿度・降雪現象なしデータ

 お互いのクラスが入り混じった分布になっていてわかりにくいのですが、なんとなく次のようなラインで分類できそうな気がします。

Python 北海道室蘭市の気温・湿度・降雪データ02

 これは、あくまで直感で適当に引いた線です。

平均交差エントロピー誤差関数

 前回と同じようにロジスティック回帰によって分類を試みますが、今回は入力データが 2 変数となっています。気温 $x$ と湿度 $y$ が与えられたときに「降雪現象なし」$(t=1)$ となる確率 $P(t=1|x,y)$ を 2 変数ロジスティック関数
 
\[P(t=1|x,y)=L(x,y)=\frac{1}{1+\exp\{-(ax+by+c)\}}\]
で表すことにして、パラメータベクトル $\boldsymbol{q}=(a,\ b,\ c)$ を最適化するために、平均交差エントロピー誤差
 
\[E(\boldsymbol{q})=\frac{1}{N}\sum_{k=0}^{N-1}\{-t_k\log L_k-(1-t_k)\log (1-L_k)\}\]
の最小値を求めます。こちらの記事 と同じように計算すると、$E$ の勾配 (gradient) は
 
\[\nabla E(\boldsymbol{q})=\begin{bmatrix}\cfrac{\partial E}{\partial a}\\
\cfrac{\partial E}{\partial b}\\
\cfrac{\partial E}{\partial c}\end{bmatrix}
=\frac{1}{N}\begin{bmatrix}\displaystyle\sum_{k=0}^{N-1}(L_k-t_k)x_k\\
\displaystyle\sum_{k=0}^{N-1}(L_k-t_k)y_k\\
\displaystyle\sum_{k=0}^{N-1}(L_k-t_k)\end{bmatrix}\]
で与えられることがわかります。
 

いちばんやさしいPythonの教本 人気講師が教える基礎からサーバサイド開発まで 「いちばんやさしい教本」シリーズ

関数の実装

 2 変数ロジスティック関数とその勾配、平均交差エントロピー誤差関数を実装しておきます。

# リストM7-A-2

# 2変数ロジスティック関数
def logistic_2(x, y, q):
    f = 1 / (1 + np.exp(-(q[0] * x + q[1] * y + q[2])))
    return f

# 平均交差エントロピー誤差関数
def cee_logistic_2(q, x, y, t):
    f = logistic_2(x, y, q)
    c = - np.sum((t * np.log(f) + (1 - t) * np.log(1 - f)))
    return c / len(x)

# 平均交差エントロピー誤差関数の勾配
def grad_cee_logistic_2(q, x, y, t):
    f = logistic_2(x, y, q)
    grad = np.zeros(3)
    grad[0] = np.sum((f - t) * x)
    grad[1] = np.sum((f - t) * y)
    grad[2] = np.sum(f - t)
    return grad / len(x)

 
 scipy.minimize() を使ってパラメータを最適化します。

# リストM7-A-3

# パラメータの初期値
q0 = [2, -5, 5]

# 平均交差エントロピー誤差関数が最小値をとるパラメータを探索
fit = minimize(cee_logistic_2, q0, args = (X, Y2, T),
               jac = grad_cee_logistic_2, method = "CG")

# パラメータを取得
q = fit.x

print("(a, b, c) = ({0:.3f}, {1:.3f}, {2:.3f})".format(q[0], q[1], q[2]))
(a, b, c) = (0.467, -15.486, 12.765)

 
 mpl_toolkits.mplot3d を使って、2変数ロジスティック関数を 3D 表示してみましょう。

# リストM7-A-3-1

# 格子点生成関数を定義
def mesh_2d(xr, yr, n):
    x = np.linspace(xr[0], xr[1], n)
    y = np.linspace(yr[0], yr[1], n)
    X, Y = np.meshgrid(x, y)
    return X, Y

# FigureとAxes
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection="3d")
ax.set_title("Logistic function of 2 Variables", fontsize = 16,
             position = (0.5, 1.1))
ax.set_xlabel("x", size = 15)
ax.set_ylabel("y", size = 15)
ax.set_zlabel("P (t=1)", size = 15)
ax.set_zlim([0, 1])

# 格子点を作成
XX, YY = mesh_2d([-10, 10], [0.5, 1.0], 129)

# 高度を計算
TT = logistic_2(XX, YY, q)

# 視点の設定
ax.view_init(45, 45)

# 曲面を描画
ax.plot_surface(XX, YY, TT, color = "green")

plt.show()

2変数ロジスティック関数の曲面
 
 最適化されたロジスティック関数が 0.5 となる $x$ と $y$、すなわち $L(x,y)=0.5$ を表す等高線が決定境界となります。

# リストM7-A-4
def mesh_2d(xr, yr, n):
    x = np.linspace(xr[0], xr[1], n)
    y = np.linspace(yr[0], yr[1], n)
    X, Y = np.meshgrid(x, y)
    return X, Y

# 等高線描画データ
xr = [-10, 10]
yr = [0.5, 1.0]
n = 257
XX, YY = mesh_2d(xr, yr, n)
TT = logistic_2(XX, YY, q)

# Figureと3DAxesの設定
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111)
ax.set_xlabel("Temperature [${}^\circ$C]", fontsize = 14)
ax.set_ylabel("Daily mean humidity [%]", fontsize = 14)
ax.set_xlim([xr[0], xr[1]])
ax.set_ylim([yr[0], yr[1]])

# 等高線をプロット
ct = ax.contour(XX, YY, TT, colors = "black", levels = [0.5])
ax.clabel(ct, fmt = "%0.1f", fontsize = 10, zorder = 1)

# 降雪現象が観測された日の最低気温と平均湿度
X_0 = X[T == 0]
Y_0 = Y2[T == 0]

# 降雪現象が観測されなかった日の最低気温と平均湿度
X_1 = X[T == 1]
Y_1 = Y2[T == 1]

# 降雪現象が観測された日のデータを青色で表示
ax.scatter(X_0, Y_0, color = "blue",
           s = 16, label = "snow", zorder = 2)

# 降雪現象が観測されなかった日のデータを赤色で表示
ax.scatter(X_1, Y_1, color = "red",
           s = 16, label = "rain", zorder = 2)

ax.set_yticklabels(np.arange(50, 110, 10))

ax.legend()

Python 等高線によるクラス分類

 今回の記事で分類の基本問題は終わりです。
 次回からいよいよニューラルネットワークについて学びます。