久しぶりの機械学習。
データセット
Boston house-prices (ボストン住宅価格データセット)
本題に入る前にデータセットについて。
毎回題材に合いそうなデータセットを用意するのも難しい場合があるので、公開されているデータセットを活用してみます。
scikit-learnのBoston house-prices (ボストン住宅価格データセット)と言われているものを使用します。
この他にもいくつか公開されているようです。
7. Dataset loading utilities — scikit-learn 0.22.1 documentation
データの中身は以下のような内容になっているようです。
■データ
| カラム名 | 説明 |
|---|---|
| CRIM | 人口 1 人当たりの犯罪発生数 |
| ZN | 25,000 平方フィート以上の住居区画の占める割合 |
| INDUS | 小売業以外の商業が占める面積の割合 |
| CHAS | チャールズ川によるダミー変数 (1: 川の周辺, 0: それ以外) |
| NOX | NOx の濃度 |
| RM | 住居の平均部屋数 |
| AGE | 1940 年より前に建てられた物件の割合 |
| DIS | 5つのボストン市の雇用施設からの距離 |
| RAD | 環状高速道路へのアクセスしやすさ |
| TAX | $10,000 ドルあたりの不動産税率の総計 |
| PTRATIO | 町毎の児童と教師の比率 |
| B | 町毎の黒人 (Bk) の比率を次の式で表したもの。 1000(Bk – 0.63)^2 |
| LSTAT | 給与の低い職業に従事する人口の割合 (%) |
| MEDV | 所有者が占有している家屋の$ 1000単位の中央値 |
全14カラムのうち、最初の13カラムが説明変数、最後のMEDVが目的変数となります。
実際にデータを取得してプロットしてみます。
・例
import numpy as npy import matplotlib.pyplot as plt %matplotlib inline #ボストン住宅価格データセットの読み込み from sklearn.datasets import load_boston boston = load_boston() # 説明変数 X_array = boston.data # 1カラム目（CRIM）を取得 X=X_array[:,0] # 目的変数 y_array = boston.target Y=y_array plt.figure(figsize=(4,4)) plt.plot(X,Y,marker='o',linestyle='None',markeredgecolor='black',color='cornflowerblue')
今回は説明変数として「CRIM」のみを使用しますのであえて取り出しています。
本来はデータセットまるごと使用するものでしょう。
最小二乗法
まずは上記の「CRIM」と「MEDV」の相関関係をこれまでも行った最小二乗法によって導出します。
・使用データを絞る
import numpy as npy import matplotlib.pyplot as plt %matplotlib inline #ボストン住宅価格データセットの読み込み from sklearn.datasets import load_boston boston = load_boston() #説明変数 X_array = boston.data # 1カラム目（CRIM）を取得 X=X_array[:,0] #目的変数 y_array = boston.target Y=y_array plt.figure(figsize=(4,4)) # 0番目から15番目のデータをプロット plt.plot(X[0:15],Y[0:15],marker='o',linestyle='None',markeredgecolor='black',color='cornflowerblue')
・出力結果
これらのデータから線形回帰モデルとなる直線を以前と同様の方法で求めます。
・ライブラリによる最小二乗法
# sklearn.linear_model.LinearRegression クラスを読み込み from sklearn import linear_model import pandas as pd import numpy as npy import matplotlib.pyplot as plt #ボストン住宅価格データセットの読み込み from sklearn.datasets import load_boston clf = linear_model.LinearRegression() boston = load_boston() #説明変数 X_array = boston.data #目的変数 y_array = boston.target X = npy.array([X_array[:,0][0:15]]).T Y = y_array[0:15] # 予測モデルを作成 clf.fit(X, Y) # 回帰係数 print(clf.coef_) # 切片 (誤差) print(clf.intercept_) # 決定係数 print(clf.score(X, Y)) # 散布図 plt.scatter(X, Y) # 回帰直線 plt.plot(X, clf.predict(X))
・出力結果
[-16.03069468]
26.562651256893634
0.23087204323047003
プロットの位置からも予想がつく通り、これらのデータから求められる直線は誤差が大きいです。
そして、犯罪率が高くなったとしても価格がマイナス値に振れるというのは現実的ではないでしょう。
線形基底関数モデル
線形基底関数モデル
これまでは線形回帰モデルは全て直線で考えてきました。
しかし現実のデータは当然ながら直線モデルで十分に表現できるものばかりではありません。
そこで、今回は曲線モデルについて考えてみます。
数学の曲線というと二次関数、三次関数、対数関数、指数関数などがありますが、基底関数（φ( x)）はこれらのベースとなる関数のことを指しています。
線形基底関数モデルとは、線形回帰モデルのxを基底関数φ( x)で置き換えたものになります。
ガウス基底関数
今回はガウス関数を基底関数としたガウス基底関数モデルを扱ってみます。
ここで先にガウス関数の性質について軽く確認しておきましょう。
・ガウス関数
# ライブラリのインポート
import numpy as npy
import matplotlib.pyplot as plt
# Jupiter Notebookで結果を表示するためのおまじない
%matplotlib inline
# 関数の定義
def f1(x,mu,s):
return npy.exp(-(x-mu)**2/(2*s**2))
# X軸に0から9までの配列
xline = npy.linspace(-10,10,500)
yline1 = f1(xline,0,1)
yline2 = f1(xline,1,1)
yline3 = f1(xline,0,2)
# pltに要素を設定する
plt.plot(xline, yline1, color ='black', label='$func1$')
plt.plot(xline, yline2, label='$func2$')
plt.plot(xline, yline3, label='$func3$')
# 凡例を表示する
plt.legend(loc="upper left")
# タイトルを表示する
plt.title('$gaussModel$')
# X軸に名前を付ける
plt.xlabel('$x$')
# Y軸に名前を付ける
plt.ylabel('$y$')
# グリッドを表示する
plt.grid(True)
# グラフを描写する
plt.show()・出力結果
形が分かると何となくイメージがつく気がしますが正規分布のモデルになります。
関数の性質としては「μ」が分布の位置（中央値）を表しており、「s」がグラフの幅を表していることが分かります。
ガウス関数を用いた線形基底関数モデル
ガウス関数を用いて線形基底関数モデルを考えていきます。
元々の線形回帰モデルは以下のようなものでした。
ｘを基底関数Φに置き換えることにより線形基底関数モデルとして表現することができます。
以下はM=4の場合の線形基底関数モデルです。各関数Φはガウス関数を表します。
線形基底関数モデルは重みパラメータである「w」を各基底関数に乗算し足し合わせたものになります。
「w4」については他のパラメータと異なり関数Φがついていないため「Φ4(x)=1」というダミー関数をつけることで対応します。
これは前回にも適用した考え方です。
これにより線形基底関数yは以下のように表現することができます。
このように表現できたことにより、平均二乗誤差も簡潔な表現をすることができます。
また、パラメータwについても以前求めた解析解を置き換えることで基底関数モデルの場合でもwを求めることができます。
ここでのΦは計画行列と呼ばれるもので以下を表しています。
またこれは入力パラメータが複数の場合にも対応できる多次元入力モデルのベクトルxを用いた表現にするのが普通のようです。
線形基底関数モデルの実装
全体像
最適なパラメータwを求めるためにPythonで実装を行います。
先に実装例を載せます。
・実装例
import numpy as npy
import matplotlib.pyplot as plt
%matplotlib inline
#ボストン住宅価格データセットの読み込み
from sklearn.datasets import load_boston
def gauss( x, mu, s):
#return 1+s/(x+mu)
return npy.exp(-(x-mu)**2 / (2*s**2))
def gauss_func(w,x):
m = len(w) -1
mu = npy.linspace(0,0.5,m)
s = mu[1] -mu[0]
y = npy.zeros_like(x)
for j in range(m):
y = y + w[j]*gauss(x, mu[j], s)
y = y + w[m]
return y
def mse_gauss_func(x, t, w):
y = gauss_func(w, x)
mse = npy.mean((y-t)**2)
return mse
def fit_gauss_func(x,t,m):
mu = npy.linspace(0,0.5,m)
s = mu[1] - mu[0]
n = x.shape[0]
phi = npy.ones((n,m+1))
for j in range(m):
phi[:,j] = gauss(x,mu[j],s)
phi_T = npy.transpose(phi)
b = npy.linalg.inv(phi_T.dot(phi))
c = b.dot(phi_T)
w = c.dot(t)
return w
def show_gauss_func(w):
xb = npy.linspace(X_min, X_max, 30)
y = gauss_func(w, xb)
plt.plot(xb, y, c=[.5, .5, .5], lw=4)
boston = load_boston()
#説明変数
X_array = boston.data
#目的変数
y_array = boston.target
X = npy.array([X_array[:,0][0:40]])
Y = y_array[0:40]
xline = X[0]
yline = Y
X_min = 0
X_max = max(max(X))
X_n = len(X[0])
plt.figure(figsize=(4,4))
M = 4
W = fit_gauss_func(xline,yline,M)
show_gauss_func(W)
plt.plot(xline,yline,marker='o',linestyle='None',markeredgecolor='black',color='cornflowerblue')
plt.xlim(X_min, X_max)
plt.grid(True)
mse = mse_gauss_func(X[0],Y[0],W)
print('W=' + str(npy.round(W,1)))
print("SD={0:.2f}cm".format(npy.sqrt(mse)))
plt.show()・出力結果
W=[16.8 -5.9 -1.5 5.6 15.7]
SD=6.02cm
以下では関数を一つづつ確認していきます。
ガウス関数
先ほども示したガウス関数です。
def gauss( x, mu, s):
#return 1+s/(x+mu)
return npy.exp(-(x-mu)**2 / (2*s**2))
パラメータwを求める関数
Jを最小にする最適なパラメータwを求めるための関数です。
実装中の変数「phi」は計画行列Φを表しています。
後半部分ではΦを転置させ行列の乗算を行うことで最終的にパラメータwを求めています。
ガウス関数の返却値はxの数に対応するベクトルを返却します。
「phi」は行列なので行ごとにガウス関数の返却値を格納しています。
def fit_gauss_func(x,t,m):
mu = npy.linspace(0,0.5,m)
s = mu[1] - mu[0]
n = x.shape[0]
phi = npy.ones((n,m+1))
for j in range(m):
phi[:,j] = gauss(x,mu[j],s)
phi_T = npy.transpose(phi)
b = npy.linalg.inv(phi_T.dot(phi))
c = b.dot(phi_T)
w = c.dot(t)
return w
線形基底関数モデル
xとwを入力パラメータとしてyを返却する以下の式を表しています。
実装中の「mu」はμを表し位置を0から0.5の間をm個に分割してyに加算していきます（総和を求める）。
def gauss_func(w,x):
m = len(w) -1
mu = npy.linspace(0,0.5,m)
s = mu[1] -mu[0]
y = npy.zeros_like(x)
for j in range(m):
y = y + w[j]*gauss(x, mu[j], s)
y = y + w[m]
return y
平均二乗誤差Jを求める関数
式で表されている通り先ほど作成したgauss_funcを用いてyを求めている。
def mse_gauss_func(x, t, w):
y = gauss_func(w, x)
mse = npy.mean((y-t)**2)
return mse
線形基底関数モデルを描写する関数
線形基底関数モデルを呼び出しグラフへの描写を行います。
def show_gauss_func(w):
xb = npy.linspace(X_min, X_max, 30)
y = gauss_func(w, xb)
plt.plot(xb, y, c=[.5, .5, .5], lw=4)
出力結果は・・・
概ね散布図に従ったグラフが描写できています。
正直深く分かっていない部分もあるのでパラメータの設定にかなり苦労しました。
単純な直線ではないのでかなりそれらしい形になっていると思います。これぐらいできると入力データがよければ実践的にも使用可能なレベルになってきたようにも感じます。
ほんとに難しくなってきた・・・
・今回のソース
python_dev/Supervised_learning5.ipynb at master · wantanblog/python_dev · GitHub
