正則化の意味～制約条件付最適化

2020-05-05 / tau / コメントする

疑　問

L2ノルムによるRidge回帰やL1ノルムによるLasso回帰に登場する罰則項。Lpノルムで表現すると次のようになる。

(1) $\begin{equation*} L(\boldsymbol{w}) = E(\boldsymbol{w}) + \alpha \sum_i |w_i|^p \end{equation*}$

ここでEは正則化がない場合の下の損失関数、w=(w₁, … , w_m)は特徴量に対する重みでLは正則化考慮後の損失関数。

このLを最小となるようなwを計算していくことになるのだが、多くの文献やサイトでL1正則化やL2正則化を説明するのに次のような図を使って、制約条件付きの最小化問題としている。

確かに式(1)の形はLagrrangeの未定乗数法に似ているのだが、ここで分からなくなるのが、αがハイパーパラメーターとして設定される点である。未定乗数法の場合、wとαを未知数として、Lを偏微分して解いていくが、ここではαは変数としては用いられない。

考え方

まず、以下のように重みwに対する制約を設定する。

(2) $\begin{equation*} \sum_i |w_i|^p \le \frac{1}{\mu} \qquad (\mu > 0) \end{equation*}$

これは、重みのノルムが高々1/μであることを意味している。これを0以下の制約条件として以下のように変形する。

(3) $\begin{equation*} g(\boldsymbol{w}) = \mu \sum_i |w_i|^p - 1 \le 0 \end{equation*}$

一般化すると $a\sum |w_i|^p + b$ の形になるが、不等式の両辺をaで割ることで、このような表現とすることは差し支えない。

これを制約条件として、元々の損失関数E(w)を最小化するため、λを導入して不等式制約条件付きのLagrangeの未定乗数法の問題を立てる。

(4) $\begin{align*} &\mathrm{minimize} \quad E(\boldsymbol{w}) \\ &\mathrm{subject~to} \quad g(\boldsymbol{w}) = \mu \sum_i |w_i|^p - 1 \le 0 \\ &\rightarrow \quad M(\boldsymbol{w}, \lambda) = E(\boldsymbol{w}) - \lambda \left( \mu \sum_i |w_i|^p - 1 \right) \end{align*}$

ここでKKT条件からλ ≤ 0となるので、ξ = −λ ≥ 0を導入して、以下のように表す。

(5) $\begin{align*} M(\boldsymbol{w}, \lambda) = E(\boldsymbol{w}) + \xi \left( \mu \sum_i |w_i|^p - 1 \right) \end{align*}$

以降、冒頭の正則化の式と上記の最小化問題の式を並べる。

(6) $\begin{align*} &L(\boldsymbol{w}) = E(\boldsymbol{w}) + \alpha \sum_i |w_i|^p \\ &M(\boldsymbol{w}, \xi) = E(\boldsymbol{w}) + \xi \left( \mu \sum_i |w_i|^p - 1 \right) \end{align*}$

これらをw_iで偏微分する。

(7) $\begin{align*} &\frac{\partial L}{\partial w_i} = \frac{\partial E}{\partial w_i} + \alpha \frac{d |w_i|^p}{d w_i} \\ &\frac{\partial M}{\partial w_i} = \frac{\partial E}{\partial w_i} + \xi \mu \frac{d |w_i|^p}{d w_i} \end{align*}$

上式から、wに関しては2つの問題は等価と言える。しかしながら、Lの最小化についてはwの解が求まるのに対して、Mの最小化で解を確定するためには、ξを未知数としてこれについても偏微分をとり、未知数と方程式を1つ加えなければならない。どうしてLの最小化だけで最適解が定まるのか。

解　釈

ところで、αはハイパーパラメーターとして外部で設定するのであった。Lに関してはw_iの個数分の未知数と方程式となり、wの解が確定する。すなわち、元の正則化の式を解くだけで「最適化された」係数のセットがw空間において1つに定まる。

ここからがミソで、wが求められるとこれを満たす制約条件の境界が定まり、同時に|w|^pに対する制約も定まるというのがどうもこの計算の仕組みらしい。

そこで、αを設定した後のながれを追ってみる。

αを設定してLを最適化する
w = (w₁, … , w_m)が定まる
点wを通る制約条件の境界g(w) = 0が定まる → Σ|w_i|^p = μ⁻¹

すなわち、αを設定してLを最適化することで制約条件の境界g = 0が決定され、重みに対する制約の強さμも決まることになる。

なお、α = 0として罰則項の効果をなくす場合のほかは必ず制約条件の境界上で解を求めることになる。一般的なLagrangeの未定乗数法のように、与えられた制約条件と目的関数の関係によって純粋な極値問題となるということはない(後述する∇Eと∇gのgradientの向きの関係からも確認できる)。

ξやμはどのように定まるのか

本来なら生真面目に未定乗数法を解いて求めるはずのξはどこへ行ったのか。

式(7)のノルム項の比較よりα = ξμ
αは事前に設定した値で、μはLの計算結果から導出されることから、ξ = μ/αも計算可能
なおMについて、点wにおいては∇E = −ξμ∇gが成り立っており、この点でEとgのgradientが平行(−ξμの符号から、両ベクトルは逆向き)

すなわち、ξはμと組み合わさって、境界条件における目的関数Eと制約条件gのgradientの比となっている。

μについては、α = ξμの関係から、αを大きくするとμも大きくなる傾向になりそうである。それにつれて式(2)から重みをより小さく制約する方向に効きそうである。ただしこの場合、ξも変動するので単純な比例関係にはならず、αの値を試行錯誤で変化させて学習効果を確認することになる。

参考サイト

本記事をまとめるにあたって、下記サイトが大変参考になったことに感謝したい。

過学習を防ぐ「正則化」とは？～○×（まるぺけ）つくろーどっとコム

Ridge回帰の理解

2020-04-26 / tau / コメントする

定義

Ridge回帰は多重回帰の損失関数に罰則項としてL2正則化項を加味する。正則化の意味についてはこちらに詳しくまとめている。

L2ノルムは原点からのユークリッド距離。

(1) $\begin{equation*} \| \boldsymbol{w} \| _2 = \sqrt{w_1 ^2 + \cdots + w_m^2} \end{equation*}$

ただしリッジ回帰では、根号の中の二乗項で計算する。

(2) $\begin{equation*} \mathrm{minimize} \quad \sum_{i=1}^n (y_i - \hat{y}_i) + \alpha \sum_{j=1}^m w_j^2 \end{equation*}$

定式化

最小化すべき関数は、

(3) $\begin{align*} L &= \sum_{i=1}^n ( \hat{y}_i - y_i )^2 + \alpha ({w_1}^2 + \cdots + {w_2}^2) \\ &= \sum ( w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i )^2 + \alpha ({w_1}^2 + \cdots + {w_m}^2) \end{align*}$

重み係数を計算するために、それぞれで偏微分してゼロとする。

(4) $\begin{align*} \frac{\partial L}{\partial w_0} &= 2 \sum (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) = 0 \\ \frac{\partial L}{\partial w_1} &= 2 \sum x_{1i} (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) + 2 \alpha w_1 = 0 \\ \vdots\\ \frac{\partial L}{\partial w_m} &= 2 \sum x_{mi} (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) + 2 \alpha w_m = 0\\ \end{align*}$

その結果得られる連立方程式は以下の通り。

(5) $\begin{align*} n w_0 + w_1 \sum x_{1i} + \cdots + w_m \sum x_{mi} &= \sum y_i \\ w_0 \sum x_{1i} + w_1 \left( \sum {x_{1i}}^2 + \alpha \right) + \cdots + w_m \sum x_{1i} x_{mi} &= \sum x_{1i} y_i \\ \vdots \\ w_0 \sum x_{mi} + w_1 \sum x_{1i} x_{mi} + \cdots+ w_m \left( \sum {x_{mi}}^2 + \alpha \right) &= \sum x_{mi} y_i \\ \end{align*}$

ここでそれぞれの和を記号Sと添字で表し、さらに行列表示すると以下の通り。

(6) $\begin{equation*} \left[ \begin{array}{cccc} n & S_1 & \cdots & S_m \\ S_1 & S_{11} + \alpha & & S_{1m} \\ \vdots & \vdots & & \vdots \\ S_m & S_{m1} & \cdots & S_{mm} + \alpha \end{array} \right] \left[ \begin{array}{c} w_0 \\ w_1 \\ \vdots \\ w_m \end{array} \right] = \left[ \begin{array}{c} S_y \\S_{1y} \\ \vdots \\ S_{my} \end{array} \right] \end{equation*}$

ここで $w_0$ を消去して、以下の連立方程式を得る。

(7) $\begin{align*} &\left[ \begin{array}{ccc} ( S_{11} + \alpha ) - \dfrac{{S_1}^2}{n} & \cdots & S_{1m} - \dfrac{S_1 S_m}{n} \\ \vdots & & \vdots \\ S_{m1} - \dfrac{S_m S_1}{n} & \cdots & ( S_{mm} + \alpha )- \dfrac{{S_2}^2}{n} \end{array} \right] \left[ \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right] \\&= \left[ \begin{array}{c} S_{1y} - \dfrac{S_1 S_y}{n} \\ \vdots \\ S_{my} - \dfrac{S_m S_y}{n} \end{array} \right] \end{align*}$

これを分散・共分散で表すと、

(8) $\begin{equation*} \left[ \begin{array}{ccc} V_{11} + \dfrac{\alpha}{n} & \cdots & V_{1m} \\ \vdots & & \vdots \\ V_{m1} & \cdots & V_{mm} + \dfrac{\alpha}{n} \end{array} \right] \left[ \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right] = \left[ \begin{array}{c} V_{1y} \\ \vdots \\ V_{my} \end{array} \right] \end{equation*}$

ここで仮に、x_jiとx_kiが完全な線形関係にある場合を考えてみる。 $x_j = a x_i + b$ とすると、分散・共分散の性質より、

(9) $\begin{equation*} V_{jj} = a^2V_{ii}, \; V_{ji} = V_{ij} = aV_{ii}, \; V_{jk} = V_{kj} = aV_{ji} = aV_{ij} \end{equation*}$

このような場合、通常の線形回帰は多重共線性により解を持たないが、式(8)に適用すると係数行列は以下のようになる。

(10) $\begin{align*} \left[ \begin{array}{ccccccc} V_{11} + \dfrac{\alpha}{n} & \cdots & V_{1i} & \cdots & aV_{1i} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} + \dfrac{\alpha}{n} & \cdots & aV_{ii} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ aV_{i1} & \cdots & aV_{ii} & \cdots & a^2V_{ii} + \dfrac{\alpha}{n} & \cdots & aV_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & aV_{mi} & \cdots & V_{mm} + \dfrac{\alpha}{n} \end{array} \right] \end{align*}$

対角要素にαが加わることで、多重共線性が強い場合でも係数行列の行列式は正則となり、方程式は解を持つ。また正則化の効果より、αを大きな値とすることによって係数の値が小さく抑えられる。

行列による表示

式(3)の損失関数を、n個のデータに対する行列で表示すると以下の通り（重回帰の行列表現はこちらを参照）。

(11) $\begin{align*} L &= \left( \boldsymbol{Xw} - \boldsymbol{y} \right)^T \left( \boldsymbol{Xw} - \boldsymbol{y} \right) + \alpha \boldsymbol{w}^T \boldsymbol{w} \\ &= \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - 2\boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} + \alpha \boldsymbol{w}^T \boldsymbol{w} \end{align*}$

これをwで微分してLを最小とする値を求める。

(12) $\begin{gather*} \frac{dL}{d\boldsymbol{w}} = 2\boldsymbol{X}^T \boldsymbol{Xw} - 2 \boldsymbol{X}^T \boldsymbol{y} + 2 \alpha \boldsymbol{w} = \boldsymbol{0} \\ \boldsymbol{w} = \left( \boldsymbol{X}^T \boldsymbol{X} + \alpha \boldsymbol{I} \right)^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{gather*}$

waveデータセット – 線形回帰

2020-04-05 / tau / コメントする

O’Reillyの”Pythonではじめる機械学習”に載っている、scikit-learnの線形回帰のwaveデータセットへの適用の再現。

waveデータセットのサンプル数を60、train_test_split()でrandom_satet=42として、書籍と同じグラフを得る。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from mglearn.datasets import make_wave

xmin, xmax = -3, 3
ymin, ymax = -3, 3

X_source, y_source = make_wave(n_samples=60)
X_train, X_test, y_train, y_test = train_test_split(X_source, y_source, random_state=42)

linreg = LinearRegression()
linreg.fit(X_train, y_train)

X_test = np.linspace(xmin, xmax, 2).reshape(-1, 1)
y_test = linreg.predict(X_test)

print(linreg.coef_[0], linreg.intercept_)

fig, ax = plt.subplots(figsize=(6.4, 6.4))

ax.scatter(X_source, y_source, s=20)
ax.plot(X_test, y_test, c="tab:orange")

ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid()

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

ax.set_aspect('equal')

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LinearRegression

from mglearn.datasets import make_wave

xmin, xmax = -3, 3

ymin, ymax = -3, 3

X_source, y_source = make_wave(n_samples=60)

X_train, X_test, y_train, y_test = train_test_split(X_source, y_source, random_state=42)

linreg = LinearRegression()

linreg.fit(X_train, y_train)

X_test = np.linspace(xmin, xmax, 2).reshape(-1, 1)

y_test = linreg.predict(X_test)

print(linreg.coef_[0], linreg.intercept_)

fig, ax = plt.subplots(figsize=(6.4, 6.4))

ax.scatter(X_source, y_source, s=20)

ax.plot(X_test, y_test, c="tab:orange")

ax.spines['bottom'].set_position('zero')

ax.spines['left'].set_position('zero')

ax.spines['top'].set_visible(False)

ax.spines['right'].set_visible(False)

ax.grid()

ax.set_xlim(xmin, xmax)

ax.set_ylim(ymin, ymax)

ax.set_aspect('equal')

plt.show()

また、訓練結果の係数、切片とスコアについても同じ結果を得ることができる。

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from mglearn.datasets import make_wave

X, y = make_wave(n_samples=60)

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)

linreg = LinearRegression()
linreg.fit(X_train, y_train)

print("coef_     : {}".format(linreg.coef_))
print("intercept_: {}".format(linreg.intercept_))

print("training score: {:.3f}".format(linreg.score(X_train, y_train)))
print("test score    : {:.3f}".format(linreg.score(X_test, y_test)))

# coef_     : [0.39390555]
# intercept_: -0.031804343026759746
# training score: 0.670
# test score    : 0.659

import numpy as np

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LinearRegression

from mglearn.datasets import make_wave

X, y = make_wave(n_samples=60)

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)

linreg = LinearRegression()

linreg.fit(X_train, y_train)

print("coef_ : {}".format(linreg.coef_))

print("intercept_: {}".format(linreg.intercept_))

print("training score: {:.3f}".format(linreg.score(X_train, y_train)))

print("test score : {:.3f}".format(linreg.score(X_test, y_test)))

# coef_ : [0.39390555]

# intercept_: -0.031804343026759746

# training score: 0.670

# test score : 0.659

Breast cancer データセット – Logistic回帰による学習率曲線

2020-04-05 / tau / コメントする

概要

breast-cancerデータセットにscikit-learnのLogisticRegressionクラスでLogistic回帰を適用した結果。

手法全般の適用の流れはLogistic回帰～cancer～Pythonではじめる機械学習よりを参照。

ここではハイパーパラメーターを変化させたときの学習率の違いをみている。

学習率曲線

scikit-learnのLogisticRegressionクラスで、正則化のパラメーターを変化させたときの学習率曲線。同クラスにはsolver引数で収束計算のいくつかの手法が選択できるが、収束手法の違いによって意外に学習率曲線に違いが出た。またtrain_test_split()のrandom_stateを変えても違いがある。569のデータセットで訓練データとテストデータを分けてもいるが、その程度では結構ばらつきが出るということかもしれない。

まず、random_state=0とした場合の、4つの収束手法における学習率曲線を示す。L-BFGSは準ニュートン法の1つらしいので、Newton-CGと同じ傾向であるのは頷ける。SAG(Stochastic Average Gradient)はまた違った計算方法のようで、他の手法と随分挙動が異なる。収束回数はmax_iter=10000で設定していて、これくらいでも計算回数オーバーの警告がいくつか出る。回数をこれより2オーダー多くしても、状況はあまり変わらない。

random_state=11としてみると、liblinearでは大きく違わないが、他の3つの手法では傾向が違っていて、特にsagを用いた場合は訓練データの学習率の方がテストデータの学習率よりも低くなっている。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_breast_cancer

ds = load_breast_cancer()

df = pd.DataFrame(ds.data, columns=ds.feature_names)

X_train, X_test, y_train, y_test = \
    train_test_split(df, ds.target, stratify=ds.target, random_state=0)

C_sup = np.linspace(5, -4, 20)
C_val = 10**C_sup

solvers = ['liblinear', 'lbfgs', 'newton-cg', 'sag']


fig, axs = plt.subplots(2, 2, figsize=(8, 8))
axs_1d = axs.reshape(-1)

for ax, solver in zip(axs_1d, solvers):
    train_scores = np.empty(0)
    test_scores = np.empty(0)
    for C in C_val:
        logreg = LogisticRegression(C=C, solver=solver, max_iter=10000)
        logreg.fit(X_train, y_train)
        train_scores = np.append(train_scores, logreg.score(X_train, y_train))
        test_scores = np.append(test_scores, logreg.score(X_test, y_test))

    ax.plot(C_val, train_scores, label="Training scores")
    ax.plot(C_val, test_scores, label="Test scores")

    ax.set_xscale('log')
    ax.set_ylim(0.9, 1)
    ax.grid(True)
    ax.legend()
    ax.set_title(solver)

plt.show()

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.model_selection import train_test_split

from sklearn.linear_model import LogisticRegression

from sklearn.datasets import load_breast_cancer

ds = load_breast_cancer()

df = pd.DataFrame(ds.data, columns=ds.feature_names)

X_train, X_test, y_train, y_test = \

train_test_split(df, ds.target, stratify=ds.target, random_state=0)

C_sup = np.linspace(5, -4, 20)

C_val = 10**C_sup

solvers = ['liblinear', 'lbfgs', 'newton-cg', 'sag']

fig, axs = plt.subplots(2, 2, figsize=(8, 8))

axs_1d = axs.reshape(-1)

for ax, solver in zip(axs_1d, solvers):

train_scores = np.empty(0)

test_scores = np.empty(0)

for C in C_val:

logreg = LogisticRegression(C=C, solver=solver, max_iter=10000)

logreg.fit(X_train, y_train)

train_scores = np.append(train_scores, logreg.score(X_train, y_train))

test_scores = np.append(test_scores, logreg.score(X_test, y_test))

ax.plot(C_val, train_scores, label="Training scores")

ax.plot(C_val, test_scores, label="Test scores")

ax.set_xscale('log')

ax.set_ylim(0.9, 1)

ax.grid(True)

ax.legend()

ax.set_title(solver)

plt.show()

Logistic回帰

2020-04-04 / tau / コメントする

概要

適用対象となる問題

ロジスティック回帰(Logistic regression)は「回帰」という名称だが、その機能は2クラス分類である。複数の特徴量に対して、2つのクラス1/0の何れになるかを判定する。

たとえば旅行会社の会員顧客20人の年齢と、温泉(SPA)とレジャーランド(LSR)のどちらを選んだかというデータが以下のように得られているとき、新たな顧客にどちらを勧めればよいか、といったような問題。

上のデータの”is SPA”列は、温泉を選んだ場合に1、レジャーランドを選んだ場合に0としている。これを散布図にすると以下の通り。

年齢が高いほど温泉を選ぶ傾向があるのは分かるが、「年齢が与えられたときに温泉とレジャーランドのどちらを選ぶか」というモデルを導出するのが目的。

(1) $\begin{equation*} (x_1, \ldots, x_m) \rightarrow \left \{ \begin{array}{cc} 1 \\ 0 \end{array} \right. \end{equation*}$

モデルを線形モデルとし、特徴量の線形和の値によってクラス1/0の何れになるかが決まるとすると

(2) $\begin{equation*} b + w_1 x_1 + \cdots + w_m x_m \rightarrow \left \{ \begin{array}{cc} 1 \\ 0 \end{array} \right. \end{equation*}$

線形回帰による場合

ここでそのままデータに対して線形回帰を適用して、たとえば回帰式の値0.5に対する年齢を境にして温泉とレジャーランドを判定してもよさそうに見える。ただ、この時の回帰式がデータに対してどのような意味を持つのかが定かではない。

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

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")
df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)
y = df['is_spa']

linreg = LinearRegression()
linreg.fit(X, y)
b = linreg.intercept_
w = linreg.coef_[0]
x_border = (0.5 - b) / w
print(x_border)

xl = 0
xr = 100
yl = b + w * xl
yr = b + w * xr

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')
ax.plot([xl, xr], [yl, yr], c='tab:orange')
ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)
ax.set_ylim(-0.1, 1.1)
ax.set_yticks(np.linspace(0, 1, 11))
ax.grid(True)

plt.show()

import os

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.linear_model import LinearRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")

df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)

y = df['is_spa']

linreg = LinearRegression()

linreg.fit(X, y)

b = linreg.intercept_

w = linreg.coef_[0]

x_border = (0.5 - b) / w

print(x_border)

xl = 0

xr = 100

yl = b + w * xl

yr = b + w * xr

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')

ax.plot([xl, xr], [yl, yr], c='tab:orange')

ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)

ax.set_ylim(-0.1, 1.1)

ax.set_yticks(np.linspace(0, 1, 11))

ax.grid(True)

plt.show()

確率分布の仮定

単純な線形回帰ではなく、特徴量の線形和に対してクラス1となる確率分布を考える。

(3) $\begin{align*} \Pr (y=1) &= f(z) = f( b + w_1 x_1 + \cdots + w_m x_m ) \\ f(0) &= 0.5 \end{align*}$

このときの確率分布の関数形として以下が必要となる。

xの(−∞, +∞)の範囲に対して、値が(0, 1)となる
xが小さい→確率が0に近くなる→クラス0と判別されやすくなる
xが大きい→確率が1に近くなる→クラス1と判別されやすくなる

これをグラフにすると、以下のような形になる。

このような関数の回帰によるクラス分類は、確率分布に用いる関数形からLogistic回帰と呼ばれている。

Logistic回帰の定式化

確率分布関数の仮定

説明変数が $x_1, \ldots, x_m$ であり、それに対してある事象が発生する／発生しないの2通りがあるとする。

いま、 $i=1, \ldots, n$ の観察対象について $x_{i1}, \ldots, x_{im}$ のデータが観測されたとする。このとき、以下のように定式化する。

(4) $\begin{align*} & z_i = b + w_1 x_{i1} + \cdots + w_m x_{im} = b + {\boldsymbol a}{\boldsymbol x} \\ & y_i = \left\{ \begin{array}{cl} 1 & \textrm{(occured)}\\ 0 & \textrm{(not occured)} \end{array} \right. \\ & \Pr (y=1) = p = \frac{1}{1 + e^{- ( b + {\boldsymbol w}{\boldsymbol x} ) }} \end{align*}$

$x_{ij} \; (j=1, \ldots, m)$ が得られたとき、その多項式を使った確率によって事象が発生することを意味している。この確率分布の式はLogistic関数と呼ばれる。

ロジスティック関数(logistic function)はシグモイド関数(sigmoid function)とも呼ばれ、その値は(0, 1)の間にあることから確率分布として採用している。

(5) $\begin{equation*} {\rm logistic} (t) = \sigma (t) = \frac{1}{1 + e^{-t}} = p \end{equation*}$

なお、ロジスティック関数の逆関数はロジット関数(logit function)と呼ばれる。

(6) $\begin{equation*} {\rm logit} (p) = \ln \frac{p}{1 - p} = t \end{equation*}$

ここでLogit関数の対数の中はオッズの形になっている。

最尤推定

$x_{ij} (i=1, \ldots , n)$ のm×n個の説明変数データと、 $y_i (1 \; \textrm{or} \; 0)$ のターゲットデータが得られたとする。

このとき、ターゲットデータが得られたパターンとなる確率は、それぞれが独立事象であるとすれば確率の積になるから、

(7) $\begin{equation*} L &= \prod_{i=1}^n \Pr (y_i = 1) ^{y_i} \Pr (y_i = 0) ^{1 - y_i} \end{equation*}$

この確率は、パラメーター $b , w_1, \ldots, w_m$ に対する尤度関数であり、その対数尤度関数は以下のようになる。

(8) $\begin{align*} \ln L &= \sum_{i=1}^n \left( y_i \ln \Pr(y_i = 1) + (1 - y_i) \ln \Pr(y_i = 0) \right) \\ &= \sum_{i=1}^n \left( y_i \ln \Pr(y_i = 1) + (1 - y_i) \ln \left(1 - \Pr(y_i = 1) \right) \right) \\ \end{align*}$

ここで、

(9) $\begin{align*} \ln \left( 1 - \Pr(y_i = 1) \right) &= \ln \left(1 - \frac{1}{1 - e^{-(b + {\boldsymbol w}{\boldsymbol x})}} \right) \\ &= \ln \frac{e^{-(b + {\boldsymbol w}{\boldsymbol x})}}{1 + e^{-(b + {\boldsymbol a}{\boldsymbol x})}} \\ &= \ln e^{-(b + {\boldsymbol w}{\boldsymbol x})} + \ln \Pr (y_i = 1 ) \\ &= -b - {\boldsymbol w}{\boldsymbol x} + \ln \Pr (y_i = 1 ) \end{align*}$

これを使って、対数尤度関数は以下のように変形される。

(10) $\begin{align*} \ln L &= \sum_{i=1}^n \left( y_i \ln \Pr(y_i = 1) + (1 - y_i) \left( -b - {\boldsymbol w}{\boldsymbol x} + \ln \Pr (y_i = 1 ) \right) \right) \\ &= \sum_{i=1}^n \left( \ln \Pr (y_i = 1 ) + (1 - y_i) (-b - {\boldsymbol w}{\boldsymbol x}) \right) \\ &= \sum_{i=1}^n \left( \ln \frac{1}{1 + e^{-(b + {\boldsymbol w}{\boldsymbol x})}} + (1 - y_i) (-b - {\boldsymbol w}{\boldsymbol x}) \right) \\ &= \sum_{i=1}^n \left( -\ln \left( 1 + e^{-(b + {\boldsymbol w}{\boldsymbol x})} \right) + (1 - y_i) (-b - {\boldsymbol w}{\boldsymbol x}) \right) \\ \end{align*}$

得られたデータに対してこの対数尤度を最大とするような ${\boldsymbol w}$ を求めるのには、通常数値計算が用いられる。

1変数の場合

定式化

線形表現は以下のようになる。

(11) $\begin{equation*} z = b + w x \end{equation*}$

ロジスティック関数による確率表現は

(12) $\begin{equation*} \Pr (y=1) = p = \frac{1}{1 + e^{-(b + w x)}} \end{equation*}$

nこのデータセット $(x_i, y_i)$ が与えられたとき、この確率分布に関する尤度関数は

(13) $\begin{equation*} L(b, w) = \prod_{i=1}^n \left( \frac{1}{1 + e^{-(b + wx_i)}} \right)^{y_i} + \left( 1 - \frac{1}{1 + e^{-(b + wx_i)}} \right)^{1 - y_i} \end{equation*}$

その対数尤度関数は

(14) $\begin{equation*} \ln L(b, w) &=& \sum_{i=1}^{n} \left( y_i \ln \frac{1}{1 + e^{-b - w x_i}} + (1 - y_i) \ln \left( 1 - \frac{1}{1 + e^{-b - w x_i}} \right) \right) \end{equation*}$

ここでデータセット(x_i, y_i)として、が与えられたとき、対数尤度を最大とするようなb, wを計算する。

LogisticRegressionモデル

温泉のデータセットに対して、scikit-learnのLogisticRegressionモデルを適用した結果が冒頭のグラフで、コードは以下の通り。

LogisticRegressionのコンストラクターの引数としてsolver='lbfgs'を指定しているが、2020年4月時点ではこれを指定しないと警告が出るため。収束アルゴリズムにはいくつかの種類があって、データサイズや複数クラス分類・1対他分類の別などによって使い分けるようだ。

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")
df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)
y = df['is_spa']

logreg = LogisticRegression(C=1e5, solver='lbfgs')
logreg.fit(X, y)
b = logreg.intercept_[0]
w = logreg.coef_[0][0]
x_border = 0.0 - b / w
print("intercept  : {}".format(b))
print("coefficient: {}".format(w))
print("x_border   : {}".format(x_border))

xl = 0
xr = 100

x_graph = np.linspace(xl, xr)
y_graph = 1 / (1 + np.exp(- b - w * x_graph))

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')
ax.plot(x_graph, y_graph, c='tab:orange')
ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)
ax.set_ylim(-0.1, 1.1)
ax.set_yticks(np.linspace(0, 1, 11))
ax.grid(True)

plt.show()

# intercept  : -13.44995797571776
# coefficient: 0.23116842117614878
# x_border   : 58.18250566961732

import os

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from sklearn.linear_model import LogisticRegression

path = os.path.join(os.path.dirname(__file__), "data_spa.csv")

df = pd.read_csv(path)

X = np.array(df['age']).reshape(-1, 1)

y = df['is_spa']

logreg = LogisticRegression(C=1e5, solver='lbfgs')

logreg.fit(X, y)

b = logreg.intercept_[0]

w = logreg.coef_[0][0]

x_border = 0.0 - b / w

print("intercept : {}".format(b))

print("coefficient: {}".format(w))

print("x_border : {}".format(x_border))

xl = 0

xr = 100

x_graph = np.linspace(xl, xr)

y_graph = 1 / (1 + np.exp(- b - w * x_graph))

fig, ax = plt.subplots()

ax.scatter(X, y, c='tab:blue')

ax.plot(x_graph, y_graph, c='tab:orange')

ax.plot([x_border, x_border], [0, 1], c='k', linestyle='dashed')

ax.set_xlim(xl, xr)

ax.set_ylim(-0.1, 1.1)

ax.set_yticks(np.linspace(0, 1, 11))

ax.grid(True)

plt.show()

# intercept : -13.44995797571776

# coefficient: 0.23116842117614878

# x_border : 58.18250566961732

また、同じデータセットについてExcelで解いた例はこちら。

正則化

LogisticeRegressionモデルのコンストラクターの引数Cは正則化の強さを正の実数で指定する。大きな値を設定するほど正則化が弱く、小さくするほど正則化が効いてくる。

温泉のデータセットに対してCを変化させてみたところ、1000あたりから上はほとんど関数形が変わらず、ほぼ正則化がない状態。デフォルトのC=1(赤い線)は結構強い正則化のように思える。

興味深いのはこれらのグラフがほぼ1点で交わっており、その点での確率が0.6近いことである。温泉かレジャーランドかを選ぶ年齢と確率が正則化の程度に寄らず一定だということになるが、これが何を意味しているか、今のところよくわからない。

一方で、確率0.5に相当する年齢は正則化を強めていくにしたがって低くなっていて、これを判定基準とするとC=0.001ではほぼ全年齢で温泉が選択されることになってしまう。

さらにCが更に小さい値になると、確率曲線はどんどんなだらかになっていき、C=1e-6では全年齢にわたって0.5、すなわち温泉、レジャーランドのいずれを選ぶかは年齢に関わらず1/2の賭けのようになっている。

2変数の場合

2変数の例として、”O’REILLYの”Pythonではじめる機械学習”にあるforgeデータの例をトレースしたものをまとめた。

Logistic回帰～forgeデータ～Pythonではじめる機械学習より

特徴量の係数

式(4)によるとz_iの値が大きいほどy = 1となる確率が高くなり、小さいほどy = 0となる確率が高くなる。また、特徴量x_i ≥ 0とすると、係数w_iが正のときにはz_iを増加させる効果、負の時にはz_iを減少させる効果がある。y = 0, 1にターゲット0, 1が割り当てられ、それらがクラス0、クラス1を表すとすると、係数が正のときにはその特徴量の増加がクラス1を選択する方向に働き、係数が負の時にはクラス0をを選択する方向に働く。

以下は特徴量数2、クラス数2の簡単な例で、Feature-0は係数が正なのでクラス0を選択する確率を高め、Feature-1の係数は係数が負なのでクラス1を選択する確率を高める。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.linear_model import LogisticRegression

X, y = make_blobs(n_samples=100, n_features=2, centers=2, random_state=9)

logreg = LogisticRegression().fit(X, y)
coef = logreg.coef_.reshape(logreg.coef_.size)
f0 = np.linspace(-15, 3)
f1 = (- logreg.intercept_ - coef[0] * f0) / coef[1]
print("intercept:{}".format(logreg.intercept_))
print("coef:{}".format(coef))

fig, ax = plt.subplots()
ax.scatter(X[y==0][:, 0], X[y==0][:, 1], label="Class-0")
ax.scatter(X[y==1][:, 0], X[y==1][:, 1], label="Class-1")
ax.plot(f0, f1, c='g')
ax.set_xlim(-14, 4)
ax.set_ylim(-12, 6)
ax.set_xlabel("Feature-0")
ax.set_ylabel("Feature-1")
ax.legend()
plt.show()

# intercept:[2.11704743]
# coef:[ 0.89424625 -0.66371147]

import numpy as np

import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs

from sklearn.linear_model import LogisticRegression

X, y = make_blobs(n_samples=100, n_features=2, centers=2, random_state=9)

logreg = LogisticRegression().fit(X, y)

coef = logreg.coef_.reshape(logreg.coef_.size)

f0 = np.linspace(-15, 3)

f1 = (- logreg.intercept_ - coef[0] * f0) / coef[1]

print("intercept:{}".format(logreg.intercept_))

print("coef:{}".format(coef))

fig, ax = plt.subplots()

ax.scatter(X[y==0][:, 0], X[y==0][:, 1], label="Class-0")

ax.scatter(X[y==1][:, 0], X[y==1][:, 1], label="Class-1")

ax.plot(f0, f1, c='g')

ax.set_xlim(-14, 4)

ax.set_ylim(-12, 6)

ax.set_xlabel("Feature-0")

ax.set_ylabel("Feature-1")

ax.legend()

plt.show()

# intercept:[2.11704743]

# coef:[ 0.89424625 -0.66371147]

Logistic回帰～1変数・Excelによる解

2020-04-04 / tau / コメントする

例題のデータ

ある旅行会社の会員顧客20人の年齢と、温泉(SPA)とレジャーランド(LSR)のどちらを選んだかというデータが以下のように得られているとき、新たな顧客にどちらを勧めればより適切か。

このようなクラス分けの問題にLogistic回帰を使うのにPythonのパッケージなどによる方法もあるが、ここではExcelを使った方法を示す。

その流れは、各観光客の選択結果のカテゴリー変数と年齢から個別の尤度と合計の尤度の計算式を定義し、切片と係数の初期値を設定しておいてから、尤度が最大となるような切片・係数を求めるためにExcelのソルバーを使う。

元となるデータは、各観光客の年齢と、行先に選んだのが温泉(SPA)かレジャーランドか(LSR)の別、それらに対して温泉を選んだ場合は1、レジャーランドを選んだ場合は0となるカテゴリー変数。

計算表の準備

このデータから以下のような表を作る。各セルの意味と内容は以下の通り。

coef：線形式の切片Aと係数Bの初期値としてそれぞれ0をセットし、収束計算の結果が入る
intercept：切片の計算のために使われるデータで、全て固定値の1
prob：coefがA, Bの値の時に各顧客の年齢に対してis_spa=1となる確率で、Logistic関数の計算値
- セルの内容は計算式で=1/(1+EXP(-$A*Y-$B*Z))
- $A, $Bは固定座標を表し、全てのデータに対してこれらのセルの内容を使う
LH：is_spaの値に対する尤度(likelihood)
- セルの内容は計算式でX*LN(C)+(1-X)*LN(1-C)
MLE：全データのLHの和で、このデータセットのパターンに対する最大尤度の結果が入る

収束計算

データタブの一番右にあるソルバーに入る(ない場合はファイル→オプション→アドイン→設定からソルバーアドインにチェックを入れる)。

ソルバーのパラメーター設定ダイアログで、

目的セルを上記のDで選択
変数セルを上記のA:Bの範囲で選択。
目標値は「最大値」を選択

「解決」ボタンを押して収束計算すると、Dの値を最大化するA:Bの内容がセットされる。

この場合の結果は以下の通り

coef：-13.6562, 0.234647
MLE：-7.45298

確率0.5(線形式の値が0)を温泉とレジャーランドの閾値とするなら、それに相当する年齢は以下のように計算される。

(1) $\begin{equation*} -13.6562 + 0.234647 \times x = 0 \quad \rightarrow \quad x = \frac{13.6562 }{0.234647 } =58.2 \end{equation*}$

Pythonのscikit-learnのLogisticRegressionモデルを同じデータに適用した結果(C=1e5)は以下の通りで、かなり近い値となっている。

intercept_ = [-13.38993211]
coefficient_ = [0.23015561]

得られた係数の値を使って、以下の関数式のグラフを描いてみたのが以下の図でLogistic曲線が現れている。

Boston house pricesデータセットの俯瞰

2020-03-28 / tau / コメントする

概要

Boston house pricesデータセットは、持家の価格とその持家が属する地域に関する指標からなるデータセットで、多変量の特徴量から属性値を予想するモデルに使われる。

各特徴量の分布

データセットからBostonにおける506の地域における13の特徴量と住宅価格の中央値が得られるが、それぞれ単独の分布を見ておく。最後のMEDVは持家価格(1000ドル単位)の中央値(Median Value)。

特徴量CHASはチャールズ川の川沿いに立地しているか否かのダミー変数で、0/1の2通りの値を持つ。いくつかの特徴量は値が集中していたり、離れたところのデータが多かったりしている。

import matplotlib.pyplot as plt
from sklearn.datasets import load_boston

boston_ds = load_boston()

X = boston_ds.data
y = boston_ds.target
feature_names = boston_ds.feature_names

n_features = X.shape[1]

fig, axs = plt.subplots(3, 5, figsize=(13, 7))
plt.subplots_adjust(hspace=0.4)
axs_1d = axs.reshape(1, -11)[0]

for ax, nf in zip(axs_1d, range(n_features)):
    ax.hist(X[:, nf], ec='k')
    ax.set_title(feature_names[nf])

axs_1d[n_features].hist(y, ec='k')
axs_1d[n_features].set_title("MEDV")

axs_1d[-1].axis('off')

plt.show()

import matplotlib.pyplot as plt

from sklearn.datasets import load_boston

boston_ds = load_boston()

X = boston_ds.data

y = boston_ds.target

feature_names = boston_ds.feature_names

n_features = X.shape[1]

fig, axs = plt.subplots(3, 5, figsize=(13, 7))

plt.subplots_adjust(hspace=0.4)

axs_1d = axs.reshape(1, -11)[0]

for ax, nf in zip(axs_1d, range(n_features)):

ax.hist(X[:, nf], ec='k')

ax.set_title(feature_names[nf])

axs_1d[n_features].hist(y, ec='k')

axs_1d[n_features].set_title("MEDV")

axs_1d[-1].axis('off')

plt.show()

各特徴量と価格の関係

13の特徴量1つ1つと価格の関係を散布図で見てみる。

比較的明らかな関係がみられるのはRM(1戸あたり部屋数)とLATAT(下位層の人口比率)で、この2つは特徴量自体の分布が比較的”整っている”。

NOX(NOx濃度)も特徴量の分布はそこそこなだらかだが、散布図では強い相関とは言い難い。

AGE(古い物件の比率)とDIS(職業紹介所への距離)はそれぞれ分布が単調減少／単調増加で、特徴量の大小と価格の高低の関係はある程度予想通りだがかなりばらついている。いずれの指標についてもMDEVがある値以下で密度が高くなっているように見えるのは興味深い。

import matplotlib.pyplot as plt
from sklearn.datasets import load_boston

boston_ds = load_boston()

X = boston_ds.data
target = boston_ds.target
feature_names = boston_ds.feature_names

n_features = X.shape[1]

fig, axs = plt.subplots(3, 5, figsize=(13, 7))
fig.subplots_adjust(hspace=0.4, wspace=0.4)

axs_1d = axs.reshape(1, -1)[0]

for ax, nf in zip(axs_1d, range(n_features)):
    ax.scatter(X[:, nf], target, s=5)
    ax.set_xlabel(feature_names[nf])
    ax.set_ylabel("MDEV")

for i in range(-2, 0): axs_1d[i].axis('off')

plt.show()

import matplotlib.pyplot as plt

from sklearn.datasets import load_boston

boston_ds = load_boston()

X = boston_ds.data

target = boston_ds.target

feature_names = boston_ds.feature_names

n_features = X.shape[1]

fig, axs = plt.subplots(3, 5, figsize=(13, 7))

fig.subplots_adjust(hspace=0.4, wspace=0.4)

axs_1d = axs.reshape(1, -1)[0]

for ax, nf in zip(axs_1d, range(n_features)):

ax.scatter(X[:, nf], target, s=5)

ax.set_xlabel(feature_names[nf])

ax.set_ylabel("MDEV")

for i in range(-2, 0): axs_1d[i].axis('off')

plt.show()

2つの特徴量と価格の関係

個々の特徴量ごとの、価格との相関がある程度が明確だったRMとLSTATについて価格との関係を3次元で見てみる。

それぞれの相関がある程度明確なので、3次元でも一つの帯のようになっている。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.datasets import load_boston

boston_ds = load_boston()

X_df = pd.DataFrame(boston_ds.data, columns=boston_ds.feature_names)
x = np.array(X_df['RM'])
y = np.array(X_df['LSTAT'])
z = boston_ds.target

fig = plt.figure(figsize=(12, 4.8))

ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(x, y, z)
ax1.set_xlabel("RM")
ax1.set_ylabel("LSTAT")
ax1.set_zlabel("MDEV")

ax2 = fig.add_subplot(122, projection='3d')
ax2.scatter(x, y, z)
ax2.set_xlabel("RM")
ax2.set_ylabel("LSTAT")
ax2.set_zlabel("MDEV")

plt.show()

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from sklearn.datasets import load_boston

boston_ds = load_boston()

X_df = pd.DataFrame(boston_ds.data, columns=boston_ds.feature_names)

x = np.array(X_df['RM'])

y = np.array(X_df['LSTAT'])

z = boston_ds.target

fig = plt.figure(figsize=(12, 4.8))

ax1 = fig.add_subplot(121, projection='3d')

ax1.scatter(x, y, z)

ax1.set_xlabel("RM")

ax1.set_ylabel("LSTAT")

ax1.set_zlabel("MDEV")

ax2 = fig.add_subplot(122, projection='3d')

ax2.scatter(x, y, z)

ax2.set_xlabel("RM")

ax2.set_ylabel("LSTAT")

ax2.set_zlabel("MDEV")

plt.show()

Boston house‐pricesデータセット

2020-03-25 / tau / コメントする

概要

Boston house-pricesデータセットは、カーネギーメロン大学のStatLibライブラリーから取得したもので、持家の価格とその持家が属する地域に関する指標からなる。

ボストンの各地域にある506の持家の価格の中央値に対して、その地域の犯罪発生率やNOx濃度など13の指標が得られる。

ここではPythonのscikit-learnにあるbostonデータの使い方をまとめる。

データの取得とデータ構造

Pythonで扱う場合、scikit-learnのdatasetsモジュールにあるload_breast_cancer()でデータを取得できる。データはBunchクラスのオブジェクト。

from sklearn.datasets import load_boston

boston_ds = load_boston()

for key, value in zip(boston_ds.keys(), boston_ds.values()):
    print("{}:\n{}\n".format(key, value))

from sklearn.datasets import load_boston

boston_ds = load_boston()

for key, value in zip(boston_ds.keys(), boston_ds.values()):

print("{}:\n{}\n".format(key, value))

データセットの構造は辞書型で、506の地域に関する13の特徴量と、当該地域における持家住宅の1000ドル単位の価格などのデータ。

data:
[[6.3200e-03 1.8000e+01 2.3100e+00 ... 1.5300e+01 3.9690e+02 4.9800e+00]
 [2.7310e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9690e+02 9.1400e+00]
 [2.7290e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9283e+02 4.0300e+00]
 ...
 [6.0760e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 5.6400e+00]
 [1.0959e-01 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9345e+02 6.4800e+00]
 [4.7410e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 7.8800e+00]]

target:
[24.  21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 15.  18.9 21.7 20.4
 18.2 19.9 23.1 17.5 20.2 18.2 13.6 19.6 15.2 14.5 15.6 13.9 16.6 14.8
 18.4 21.  12.7 14.5 13.2 13.1 13.5 18.9 20.  21.  24.7 30.8 34.9 26.6
 .....
 16.7 12.  14.6 21.4 23.  23.7 25.  21.8 20.6 21.2 19.1 20.6 15.2  7.
  8.1 13.6 20.1 21.8 24.5 23.1 19.7 18.3 21.2 17.5 16.8 22.4 20.6 23.9
 22.  11.9]

feature_names:
['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO'
 'B' 'LSTAT']

DESCR:
.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

.....

   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.
   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.


filename:
C:\Users\tomo\AppData\Local\Programs\Python\Python37-32\lib\site-packages\sklearn\datasets\data\boston_house_prices.csv

data:

[[6.3200e-03 1.8000e+01 2.3100e+00 ... 1.5300e+01 3.9690e+02 4.9800e+00]

[2.7310e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9690e+02 9.1400e+00]

[2.7290e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9283e+02 4.0300e+00]

...

[6.0760e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 5.6400e+00]

[1.0959e-01 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9345e+02 6.4800e+00]

[4.7410e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 7.8800e+00]]

target:

[24. 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 15. 18.9 21.7 20.4

18.2 19.9 23.1 17.5 20.2 18.2 13.6 19.6 15.2 14.5 15.6 13.9 16.6 14.8

18.4 21. 12.7 14.5 13.2 13.1 13.5 18.9 20. 21. 24.7 30.8 34.9 26.6

.....

16.7 12. 14.6 21.4 23. 23.7 25. 21.8 20.6 21.2 19.1 20.6 15.2 7.

8.1 13.6 20.1 21.8 24.5 23.1 19.7 18.3 21.2 17.5 16.8 22.4 20.6 23.9

22. 11.9]

feature_names:

['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO'

'B' 'LSTAT']

DESCR:

.. _boston_dataset:

Boston house prices dataset

---------------------------

**Data Set Characteristics:**

:Number of Instances: 506

:Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

.....

- Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.

- Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.

filename:

C:\Users\tomo\AppData\Local\Programs\Python\Python37-32\lib\site-packages\sklearn\datasets\data\boston_house_prices.csv

データのキーは以下のようになっている。

from sklearn.datasets import load_boston

boston_ds = load_boston()

print(boston_ds.keys())

# dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])

from sklearn.datasets import load_boston

boston_ds = load_boston()

print(boston_ds.keys())

# dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])

データの内容

`'data'`～特徴量データセット

506の地域における13の指標を特徴量として格納した2次元配列。列のインデックスが特徴量の番号に対応している。

data:
[[6.3200e-03 1.8000e+01 2.3100e+00 ... 1.5300e+01 3.9690e+02 4.9800e+00]
 [2.7310e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9690e+02 9.1400e+00]
 [2.7290e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9283e+02 4.0300e+00]
 ...
 [6.0760e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 5.6400e+00]
 [1.0959e-01 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9345e+02 6.4800e+00]
 [4.7410e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 7.8800e+00]]

data:

[[6.3200e-03 1.8000e+01 2.3100e+00 ... 1.5300e+01 3.9690e+02 4.9800e+00]

[2.7310e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9690e+02 9.1400e+00]

[2.7290e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9283e+02 4.0300e+00]

...

[6.0760e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 5.6400e+00]

[1.0959e-01 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9345e+02 6.4800e+00]

[4.7410e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 7.8800e+00]]

`'target'`～住宅価格

506の地域における持家住宅の1000ドル単位の価格中央値

target:
[24.  21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 15.  18.9 21.7 20.4
 18.2 19.9 23.1 17.5 20.2 18.2 13.6 19.6 15.2 14.5 15.6 13.9 16.6 14.8
 18.4 21.  12.7 14.5 13.2 13.1 13.5 18.9 20.  21.  24.7 30.8 34.9 26.6
 .....
 16.7 12.  14.6 21.4 23.  23.7 25.  21.8 20.6 21.2 19.1 20.6 15.2  7.
  8.1 13.6 20.1 21.8 24.5 23.1 19.7 18.3 21.2 17.5 16.8 22.4 20.6 23.9
 22.  11.9]

target:

[24. 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 15. 18.9 21.7 20.4

18.2 19.9 23.1 17.5 20.2 18.2 13.6 19.6 15.2 14.5 15.6 13.9 16.6 14.8

18.4 21. 12.7 14.5 13.2 13.1 13.5 18.9 20. 21. 24.7 30.8 34.9 26.6

.....

16.7 12. 14.6 21.4 23. 23.7 25. 21.8 20.6 21.2 19.1 20.6 15.2 7.

8.1 13.6 20.1 21.8 24.5 23.1 19.7 18.3 21.2 17.5 16.8 22.4 20.6 23.9

22. 11.9]

`'feature_names'`～特徴名

13種類の特徴量の名称。

CRIM：町ごとの人口当たり犯罪率
ZN：25,000平方フィート以上の区画の住居用途地区比率
INDUS：町ごとの小売り以外の産業用途地区比率
CHAS：チャールズ川に関するダミー変数（1：川沿い、0：それ以外）
NOX：NOx濃度（10ppm単位）
RM：1戸あたり部屋数
AGE：1940年より前に建てられた持家物件の比率
DIS：ボストンの5つの職業紹介所への重みづけ平均距離
RAD：放射道路へのアクセス性
TAX：10,000ドルあたりの固定資産税総額
PTRATIO：生徒対教師の比率
B：1000(Bk – 0.63)^2（Bkは待ちにおける黒人比率）
LSTAT：下位層の人口比率(%)

feature_names:
['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO'
 'B' 'LSTAT']

feature_names:

['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO'

'B' 'LSTAT']

`'filename'`～ファイル名

CSVファイルのフルパス名が示されている。1行目にはデータ数、特徴量数が並んでおり、2行目に13の特徴量とターゲットの住宅価格、その後に506行のレコードに対する13列の特徴量と1列のターゲットデータが格納されている。このファイルにはDESCRに当たるデータは格納されていない。

'C:...\lib\site-packages\sklearn\datasets\data\boston_house_prices.csv'

1	'C:...\lib\site-packages\sklearn\datasets\data\boston_house_prices.csv'

`'DESCR'`～データセットの説明

データセットの説明。print(breast_ds_dataset['DESCR'])のようにprint文で整形表示される。

レコード数506個
属性は、13の数値／カテゴリー属性と、通常はターゲットに用いられる中央値

.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pupil-teacher ratio by town
        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
        - LSTAT    % lower status of the population
        - MEDV     Median value of owner-occupied homes in $1000's

    :Missing Attribute Values: None

    :Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/


This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980.   N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.   
     
.. topic:: References

   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.
   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.

.. _boston_dataset:

Boston house prices dataset

---------------------------

**Data Set Characteristics:**

:Number of Instances: 506

:Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

:Attribute Information (in order):

- CRIM per capita crime rate by town

- ZN proportion of residential land zoned for lots over 25,000 sq.ft.

- INDUS proportion of non-retail business acres per town

- CHAS Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)

- NOX nitric oxides concentration (parts per 10 million)

- RM average number of rooms per dwelling

- AGE proportion of owner-occupied units built prior to 1940

- DIS weighted distances to five Boston employment centres

- RAD index of accessibility to radial highways

- TAX full-value property-tax rate per $10,000

- PTRATIO pupil-teacher ratio by town

- B 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town

- LSTAT % lower status of the population

- MEDV Median value of owner-occupied homes in $1000's

:Missing Attribute Values: None

:Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.

https://archive.ics.uci.edu/ml/machine-learning-databases/housing/

This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic

prices and the demand for clean air', J. Environ. Economics & Management,

vol.5, 81-102, 1978. Used in Belsley, Kuh & Welsch, 'Regression diagnostics

...', Wiley, 1980. N.B. Various transformations are used in the table on

pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression

problems.

.. topic:: References

- Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.

データの利用

データの取得方法

bostonデータセットから各データを取り出すのに、以下の2つの方法がある。

辞書のキーを使って呼び出す（例：boston['DESCR']）
キーの文字列をプロパティーに指定する（例：boston.DESCR）

全レコードの特徴量データの取得

'data'から、506のレコードに関する13の特徴量が506行13列の2次元配列で得られる。13の特徴量は’feature_names’の13の特徴名に対応している。

from sklearn.datasets import load_boston

boston_ds = load_boston()

print(boston_ds.data)

# [[6.3200e-03 1.8000e+01 2.3100e+00 ... 1.5300e+01 3.9690e+02 4.9800e+00]
#  [2.7310e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9690e+02 9.1400e+00]
#  [2.7290e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9283e+02 4.0300e+00]
#  ...
#  [6.0760e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 5.6400e+00]
#  [1.0959e-01 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9345e+02 6.4800e+00]
#  [4.7410e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 7.8800e+00]]

from sklearn.datasets import load_boston

boston_ds = load_boston()

print(boston_ds.data)

# [[6.3200e-03 1.8000e+01 2.3100e+00 ... 1.5300e+01 3.9690e+02 4.9800e+00]

# [2.7310e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9690e+02 9.1400e+00]

# [2.7290e-02 0.0000e+00 7.0700e+00 ... 1.7800e+01 3.9283e+02 4.0300e+00]

# ...

# [6.0760e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 5.6400e+00]

# [1.0959e-01 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9345e+02 6.4800e+00]

# [4.7410e-02 0.0000e+00 1.1930e+01 ... 2.1000e+01 3.9690e+02 7.8800e+00]]

特定の特徴量のデータのみ取得

特定の特徴量に関する全レコードのデータを取り出すときにはX[:, n]の形で指定する。

from sklearn.datasets import load_boston

boston_ds = load_boston()

features = boston_ds.feature_names
X = boston_ds.data
n_feature = 10

feature = X[:, n_feature]

print("feature name : {}".format(features[n_feature]))
print("feature data :\n{}".format(feature))

# feature name : PTRATIO
# feature data :
# [15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 15.2 15.2 15.2 21.
#  21.  21.  21.  21.  21.  21.  21.  21.  21.  21.  21.  21.  21.  21.
#  21.  21.  21.  21.  21.  21.  21.  19.2 19.2 19.2 19.2 18.3 18.3 17.9
#  ...
#  20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.1 20.1
#  20.1 20.1 20.1 19.2 19.2 19.2 19.2 19.2 19.2 19.2 19.2 21.  21.  21.
#  21.  21. ]

from sklearn.datasets import load_boston

boston_ds = load_boston()

features = boston_ds.feature_names

X = boston_ds.data

n_feature = 10

feature = X[:, n_feature]

print("feature name : {}".format(features[n_feature]))

print("feature data :\n{}".format(feature))

# feature name : PTRATIO

# feature data :

# [15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 15.2 15.2 15.2 21.

# 21. 21. 21. 21. 21. 21. 21. 21. 21. 21. 21. 21. 21. 21.

# 21. 21. 21. 21. 21. 21. 21. 19.2 19.2 19.2 19.2 18.3 18.3 17.9

# ...

# 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.2 20.1 20.1

# 20.1 20.1 20.1 19.2 19.2 19.2 19.2 19.2 19.2 19.2 19.2 21. 21. 21.

# 21. 21. ]

waveデータセット – knn

2020-03-22 / tau / コメントする

概要

k-最近傍回帰の例として、scikit-learnのwaveデータにKNeighborsRegressorを適用してみた結果。

近傍点数とクラス分類の挙動

訓練データとして10個のwaveデータを訓練データとして与え、2つのテストデータの予測するのに、近傍点数を1, 2, 3と変えた場合の様子を見てみる。

近傍点数=1の場合

2つのテストデータの特徴量の値に最も近い特徴量を持つ訓練データが選ばれ、その属性値がそのままテストデータの属性値となっている。

近傍点数=2の場合

テストデータの特徴量に最も近い方から1番目、2番目の特徴量を持つ訓練データが選ばれ、それらの属性値の平均がテストデータの属性値となっている。

近傍点数=3の場合

同様に、テストデータの特徴量に最も近い3つの訓練データの属性の平均がテストデータの属性値となっている。

実行コード

上記の計算のコードは以下の通り。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsRegressor
from mglearn.datasets import make_wave

X_train, y_train = make_wave(n_samples=10)

reg = KNeighborsRegressor(n_neighbors=3)
reg.fit(X_train, y_train)

X_test = np.array([[-1], [1]])
y_pred = reg.predict(X_test)

neigh_dist, neigh_ind = reg.kneighbors(X=X_test)
print(neigh_ind)

fig, ax = plt.subplots(figsize=(8.0, 4.8))

xmin, xmax = -3, 3
ymin, ymax = -3, 3

ax.scatter(X_train, y_train, marker='o', s=20)
ax.scatter(X_test, y_pred, marker='*', s=120)

for test, pred, ind in zip(X_test, y_pred, neigh_ind):
    for neigh in ind:
        ax.plot([test, test], [ymin, ymax], c='gray', linestyle='dashed')
        ax.plot(
            [test[0], X_train[neigh, 0]], [pred, y_train[neigh]],
            color='k', linestyle='dotted')

for x, y in zip(X_train, y_train):
    ax.annotate("{:6.3f}".format(y), xy=(x[0] - 0.1, y + 0.08))
for x, y in zip(X_test, y_pred):
    ax.annotate("{:6.3f}".format(y), xy=(x[0] - 0.2, y - 0.3))

ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)

ax.set_xlabel("feature")
ax.set_ylabel("prediction")

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.neighbors import KNeighborsRegressor

from mglearn.datasets import make_wave

X_train, y_train = make_wave(n_samples=10)

reg = KNeighborsRegressor(n_neighbors=3)

reg.fit(X_train, y_train)

X_test = np.array([[-1], [1]])

y_pred = reg.predict(X_test)

neigh_dist, neigh_ind = reg.kneighbors(X=X_test)

print(neigh_ind)

fig, ax = plt.subplots(figsize=(8.0, 4.8))

xmin, xmax = -3, 3

ymin, ymax = -3, 3

ax.scatter(X_train, y_train, marker='o', s=20)

ax.scatter(X_test, y_pred, marker='*', s=120)

for test, pred, ind in zip(X_test, y_pred, neigh_ind):

for neigh in ind:

ax.plot([test, test], [ymin, ymax], c='gray', linestyle='dashed')

ax.plot(

[test[0], X_train[neigh, 0]], [pred, y_train[neigh]],

color='k', linestyle='dotted')

for x, y in zip(X_train, y_train):

ax.annotate("{:6.3f}".format(y), xy=(x[0] - 0.1, y + 0.08))

for x, y in zip(X_test, y_pred):

ax.annotate("{:6.3f}".format(y), xy=(x[0] - 0.2, y - 0.3))

ax.set_xlim(xmin, xmax)

ax.set_ylim(ymin, ymax)

ax.set_xlabel("feature")

ax.set_ylabel("prediction")

plt.show()

knnの精度

O’Reillyの”Pythonではじめる機械学習”中、KNeighborsRegressorのwaveデータに対する精度が計算されている。40サンプルのwaveデータを発生させ訓練データとテストデータに分け、テストデータに対するR²スコアが0.83となることが示されている。実際に計算してみると、確かに同じ値となる。

import numpy as np
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import train_test_split
from mglearn.datasets import make_wave

X_source, y_source = make_wave(n_samples=40)

X_train, X_test, y_train, y_test =\
    train_test_split(X_source, y_source, random_state=0)

reg = KNeighborsRegressor(n_neighbors=3)
reg.fit(X_train, y_train)
y_pred = reg.predict(X_test)

print("R^2 score:{:6.3f}".format(reg.score(X_test, y_test)))

# R^2 score: 0.834

import numpy as np

from sklearn.neighbors import KNeighborsRegressor

from sklearn.model_selection import train_test_split

from mglearn.datasets import make_wave

X_source, y_source = make_wave(n_samples=40)

X_train, X_test, y_train, y_test =\

train_test_split(X_source, y_source, random_state=0)

reg = KNeighborsRegressor(n_neighbors=3)

reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)

print("R^2 score:{:6.3f}".format(reg.score(X_test, y_test)))

# R^2 score: 0.834

これを見ると比較的高い精度のように見えるが、train_test_split()の引数random_stateを変化させてみると以下のように精度はばらつく。乱数系列が異なると精度が0.3未満の場合もあるが、全体としてみると0.6～0.7あたりとなりそうである。

import numpy as np
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import train_test_split
from mglearn.datasets import make_wave

X_source, y_source = make_wave(n_samples=40)

reg = KNeighborsRegressor(n_neighbors=3)

print("random_state -> R^2")

for random_state in range(0, 10):
    X_train, X_test, y_train, y_test =\
        train_test_split(X_source, y_source, random_state=random_state)
    reg.fit(X_train, y_train)
    y_pred = reg.predict(X_test)

    print("{} -> {:6.3f}".format(random_state, reg.score(X_test, y_test)))

# random_state -> R^2
# 0 ->  0.834
# 1 ->  0.581
# 2 ->  0.798
# 3 ->  0.281
# 4 ->  0.773
# 5 ->  0.738
# 6 ->  0.554
# 7 ->  0.494
# 8 ->  0.678
# 9 ->  0.801

import numpy as np

from sklearn.neighbors import KNeighborsRegressor

from sklearn.model_selection import train_test_split

from mglearn.datasets import make_wave

X_source, y_source = make_wave(n_samples=40)

reg = KNeighborsRegressor(n_neighbors=3)

print("random_state -> R^2")

for random_state in range(0, 10):

X_train, X_test, y_train, y_test =\

train_test_split(X_source, y_source, random_state=random_state)

reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)

print("{} -> {:6.3f}".format(random_state, reg.score(X_test, y_test)))

# random_state -> R^2

# 0 -> 0.834

# 1 -> 0.581

# 2 -> 0.798

# 3 -> 0.281

# 4 -> 0.773

# 5 -> 0.738

# 6 -> 0.554

# 7 -> 0.494

# 8 -> 0.678

# 9 -> 0.801

ためしにmake_wave(n_samples=1000)としてみると、結果は以下の通りとなり、精度は0.67程度（平均は0.677）と一定してくる。

random_state -> R^2
0 ->  0.679
1 ->  0.662
2 ->  0.682
3 ->  0.672
4 ->  0.680
5 ->  0.697
6 ->  0.712
7 ->  0.682
8 ->  0.661
9 ->  0.641

random_state -> R^2

0 -> 0.679

1 -> 0.662

2 -> 0.682

3 -> 0.672

4 -> 0.680

5 -> 0.697

6 -> 0.712

7 -> 0.682

8 -> 0.661

9 -> 0.641

予測カーブ

訓練データが少ない場合

40個のwaveデータに対して、n_neighborsを変化させたときの予測カーブを見てみる。

n_neighbors=1の時は、全ての訓練データを通るような線となる
n_neighborsが多くなるほど滑らかになる
n_neighborsがかなり大きくなると水平に近くなる
n_neighborsが訓練データ数と同じになると、予測線は水平になる（任意の特徴量に対して、全ての点の平均を計算しているため）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsRegressor
from mglearn.datasets import make_wave

X_train, y_train = make_wave(n_samples=40)
xmin = np.min(X_train[:, 0])
xmax = np.max(X_train[:, 0])
X_test = np.linspace(xmin, xmax, 200).reshape(-1, 1)

fig, axs = plt.subplots(2, 3, figsize=(11, 6.4))
plt.subplots_adjust(hspace=0.4, wspace=0.4)

n_neighbors_list=[1, 2, 8, 16, 32, 40]
axs_1d = axs.reshape(1, -1)[0]

for ax, n_neighbors in zip(axs_1d, n_neighbors_list):
    reg = KNeighborsRegressor(n_neighbors=n_neighbors)
    reg.fit(X_train, y_train)
    y_pred = reg.predict(X_test)

    ax.scatter(X_train[:, 0], y_train, zorder=2, s=20, color='tab:blue')
    ax.plot(X_test, y_pred, zorder=1, color='tab:orange')

    ax.set_title("n_neighbors={}".format(n_neighbors))
    ax.set_xlabel("feature")
    ax.set_ylabel("target")

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from sklearn.neighbors import KNeighborsRegressor

from mglearn.datasets import make_wave

X_train, y_train = make_wave(n_samples=40)

xmin = np.min(X_train[:, 0])

xmax = np.max(X_train[:, 0])

X_test = np.linspace(xmin, xmax, 200).reshape(-1, 1)

fig, axs = plt.subplots(2, 3, figsize=(11, 6.4))

plt.subplots_adjust(hspace=0.4, wspace=0.4)

n_neighbors_list=[1, 2, 8, 16, 32, 40]

axs_1d = axs.reshape(1, -1)[0]

for ax, n_neighbors in zip(axs_1d, n_neighbors_list):

reg = KNeighborsRegressor(n_neighbors=n_neighbors)

reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)

ax.scatter(X_train[:, 0], y_train, zorder=2, s=20, color='tab:blue')

ax.plot(X_test, y_pred, zorder=1, color='tab:orange')

ax.set_title("n_neighbors={}".format(n_neighbors))

ax.set_xlabel("feature")

ax.set_ylabel("target")

plt.show()

訓練データが多い場合

今度はwaveデータでn_samples=200と数を多くしてみる。データ数を多くするとその名の通り、上下に波打ちながら増加している様子が見られる。これに対してn_neighborsを変化させたのが以下の図。

n_neighbors=10～20あたりで滑らかに、かつ波打つ状況が曲線で再現されている。

n_samples=300として訓練データに200を振り分け、n_neighborsを変化させたときのスコアは以下の通り。n_neighbors=20あたりで精度が最もよさそうである。

あるデータが得られたとき、その科学的なメカニズムは置いておいて、とりあえずデータから予測値を再現したいときにはそれなりに使えるかもしれない。

n_neighbors -> R^2
5 ->  0.754
10 ->  0.788
15 ->  0.789
20 ->  0.792
25 ->  0.777
50 ->  0.737
100 ->  0.613
200 -> -0.022

n_neighbors -> R^2

5 -> 0.754

10 -> 0.788

15 -> 0.789

20 -> 0.792

25 -> 0.777

50 -> 0.737

100 -> 0.613

200 -> -0.022

import numpy as np
from sklearn.neighbors import KNeighborsRegressor
from sklearn.model_selection import train_test_split
from mglearn.datasets import make_wave

X_source, y_source = make_wave(n_samples=300)

X_train, X_test, y_train, y_test =\
    train_test_split(X_source, y_source, train_size=200, random_state=0)

n_neighbors_list = [5, 10, 15, 20, 25, 50, 100, 200]

print("n_neighbors -> R^2")

for n_neighbors in n_neighbors_list:
    reg = KNeighborsRegressor(n_neighbors=n_neighbors)
    reg.fit(X_train, y_train)
    y_pred = reg.predict(X_test)

    print("{} -> {:6.3f}".format(n_neighbors, reg.score(X_test, y_test)))

import numpy as np

from sklearn.neighbors import KNeighborsRegressor

from sklearn.model_selection import train_test_split

from mglearn.datasets import make_wave

X_source, y_source = make_wave(n_samples=300)

X_train, X_test, y_train, y_test =\

train_test_split(X_source, y_source, train_size=200, random_state=0)

n_neighbors_list = [5, 10, 15, 20, 25, 50, 100, 200]

print("n_neighbors -> R^2")

for n_neighbors in n_neighbors_list:

reg = KNeighborsRegressor(n_neighbors=n_neighbors)

reg.fit(X_train, y_train)

y_pred = reg.predict(X_test)

print("{} -> {:6.3f}".format(n_neighbors, reg.score(X_test, y_test)))

forgeデータセット – knn

2020-03-22 / tau / コメントする

概要

ここでは、Pythonのscikit-learnパッケージのKNeighborsClassifierクラスにmglearnパッケージのforgeデータを適用してknnの挙動を確認する。

近傍点数を変化させたときのクラス分類の挙動や学習率曲線についてみていく。

近傍点数によるクラス分類の挙動

近傍点数=1の場合

データセットとしてmglearnで提供されているforgeデータを用いて、近傍点数=1とした場合の、3つのテストデータのクラス判定を以下に示す。各テストデータに対して最も距離(この場合はユークリッド距離)が近い点1つが定まり、その点のクラステストデータのクラスとして決定している。

なお、いろいろなところで見かけるforgeデータセットの散布図は当該データセットの特徴量0(横軸)と特徴量1(縦軸)の最小値と最大値に合わせて表示しており、軸目盛の比率が等しくない。ここでは、距離計算に視覚上の齟齬が生じないように、縦軸と横軸の比率を同じとしている。

後の計算のために、このグラフ描画のコードを以下に示す。

import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsClassifier
from mglearn.datasets import make_forge

X, y = make_forge()

clfr = KNeighborsClassifier(n_neighbors=1)
clfr.fit(X, y)

col = ['blue', 'red']

test_points = [[9., 4.], [10., 3.], [11., 2.]]
nb_dist, nb_idx = clfr.kneighbors(test_points)
test_pred = clfr.predict(test_points)

fig, ax = plt.subplots()

ax.scatter(X[:, 0][y==0], X[:, 1][y==0], marker='o', c=col[0], label="class-0")
ax.scatter(X[:, 0][y==1], X[:, 1][y==1], marker='^', c=col[1], label="class-1")

ax.legend(loc="lower left")

for pts, cls, ids, dists in zip(test_points, test_pred, nb_idx, nb_dist):
    print(pts)
    ax.scatter(pts[0], pts[1], marker='*', s=150, c=col[cls])
    for id, dst in zip(ids, dists):
        ax.plot([pts[0], X[id, 0]], [pts[1], X[id, 1]], c='gray')
        print(" [{:7.4f}, {:7.4f}] - {:7.4f}".format(X[id, 0], X[id, 1], dst))

plt.show()

import matplotlib.pyplot as plt

from sklearn.neighbors import KNeighborsClassifier

from mglearn.datasets import make_forge

X, y = make_forge()

clfr = KNeighborsClassifier(n_neighbors=1)

clfr.fit(X, y)

col = ['blue', 'red']

test_points = [[9., 4.], [10., 3.], [11., 2.]]

nb_dist, nb_idx = clfr.kneighbors(test_points)

test_pred = clfr.predict(test_points)

fig, ax = plt.subplots()

ax.scatter(X[:, 0][y==0], X[:, 1][y==0], marker='o', c=col[0], label="class-0")

ax.scatter(X[:, 0][y==1], X[:, 1][y==1], marker='^', c=col[1], label="class-1")

ax.legend(loc="lower left")

for pts, cls, ids, dists in zip(test_points, test_pred, nb_idx, nb_dist):

print(pts)

ax.scatter(pts[0], pts[1], marker='*', s=150, c=col[cls])

for id, dst in zip(ids, dists):

ax.plot([pts[0], X[id, 0]], [pts[1], X[id, 1]], c='gray')

print(" [{:7.4f}, {:7.4f}] - {:7.4f}".format(X[id, 0], X[id, 1], dst))

plt.show()

概要は以下の通り。

5行目でforgeデータセットを準備
7行目で近傍点数を1で指定してクラス分類器を構築
8行目で訓練データとしてforgeデータを与える
12行目で3つのテストデータを準備
13行目でテストデータに対する近傍点のインデックスとテストデータまでの距離を獲得
14行目でテストデータのクラスを決定
18-19行目で訓練データの散布図を描画
23行目で、テストデータとそのクラス決定結果、クラス決定に用いられた点群のインデックス、テストデータと各点の距離を並行してループ
- 24行目でテストデータの座標を出力
- 25行目でテストデータを描画
- 26行目のループで、テストデータごとの近傍点に関する処理を実行
  - 27行目でテストデータと近傍点の間に直線を描画
  - 28行目で近傍点とテストデータからの距離を出力

出力結果は以下の通りで、各予測点に対して近傍点が1つ決定されている。

[9.0, 4.0]
 [ 8.6749,  4.4757] -  0.5762
[10.0, 3.0]
 [10.2403,  2.4554] -  0.5952
[11.0, 2.0]
 [11.5640,  1.3389] -  0.8689

[9.0, 4.0]

[ 8.6749, 4.4757] - 0.5762

[10.0, 3.0]

[10.2403, 2.4554] - 0.5952

[11.0, 2.0]

[11.5640, 1.3389] - 0.8689

近傍点数=3の場合

先の例で、コードの7行目で近傍点=3で指定してクラス分類器を構築する。

clfr = KNeighborsClassifier(n_neighbors=3)

1	clfr = KNeighborsClassifier(n_neighbors=3)

一般にknnでは、テストデータに対して複数の近傍点を指定する場合、各近傍点のクラスのうち最も多いものをテストデータのクラスとする(多数決)。

[9.0, 4.0]
 [ 8.6749,  4.4757] -  0.5762
 [ 9.4912,  4.3322] -  0.5930
 [ 8.1062,  4.2870] -  0.9387
[10.0, 3.0]
 [10.2403,  2.4554] -  0.5952
 [ 9.5017,  1.9382] -  1.1729
 [ 8.7337,  2.4916] -  1.3645
[11.0, 2.0]
 [11.5640,  1.3389] -  0.8689
 [10.2403,  2.4554] -  0.8858
 [10.0639,  0.9908] -  1.3765

[9.0, 4.0]

[ 8.6749, 4.4757] - 0.5762

[ 9.4912, 4.3322] - 0.5930

[ 8.1062, 4.2870] - 0.9387

[10.0, 3.0]

[10.2403, 2.4554] - 0.5952

[ 9.5017, 1.9382] - 1.1729

[ 8.7337, 2.4916] - 1.3645

[11.0, 2.0]

[11.5640, 1.3389] - 0.8689

[10.2403, 2.4554] - 0.8858

[10.0639, 0.9908] - 1.3765

近傍点数=2の場合

テストデータのクラスを近傍点のクラスの多数決で求めるとすると、近傍点数が偶数の時の処理が問題になる。KNeighborsClassifierの場合、偶数でクラス分類が拮抗する場合は、クラス番号が最も小さいものに割り当てられるらしい。実際、n_neighbors=2としたときの3つのテストデータのうち中央の点(10.0, 3.0)については、赤い点(10.24, 2.45)～class-1～距離0.5952の方が青い点(9.5017, 1.9382)～class-0～距離1.1729よりも距離は近いがクラス番号が0である青い点のクラスで判定されている。

[9.0, 4.0]
 [ 8.6749,  4.4757] -  0.5762
 [ 9.4912,  4.3322] -  0.5930
[10.0, 3.0]
 [10.2403,  2.4554] -  0.5952
 [ 9.5017,  1.9382] -  1.1729
[11.0, 2.0]
 [11.5640,  1.3389] -  0.8689
 [10.2403,  2.4554] -  0.8858

[9.0, 4.0]

[ 8.6749, 4.4757] - 0.5762

[ 9.4912, 4.3322] - 0.5930

[10.0, 3.0]

[10.2403, 2.4554] - 0.5952

[ 9.5017, 1.9382] - 1.1729

[11.0, 2.0]

[11.5640, 1.3389] - 0.8689

[10.2403, 2.4554] - 0.8858

偶数の点で多数決で拮抗した場合には、最も近い点のクラスで決定する、平均距離が近い方のクラスで決定するといった方法が考えられるが、この場合は必ず番号が小さなクラスが選ばれるため、若干結果に偏りがでやすいのでは、と考える。

決定境界

近傍点の数を変えた時の決定境界の変化を確認する。k近傍法はscikit-learnのKNeighborsClassifierクラスを利用する。

近傍点の数を1, 2, 3, …と変化させたときの決定境界の変化は以下の通り。

近傍点数が少ないときは訓練データにフィットするよう決定境界が複雑になるが、近傍点数が多いと決定境界は滑らかになる。特に近傍点数が訓練データの点数に等しいとき、全訓練データの多数決でクラス決定され、全領域で判定結果が同じとなる(この場合は近傍点数26が偶数なので、クラス番号の小さいclass-0で決定されている)。

この図を描画したコードを以下に示す。

7行目、引数で与えたAxesに対して決定境界を描く関数を定義
- 18行目、決定境界をcontourf()を利用して描いている
21行目、引数で与えたAxesに対してクラスごとに色分けした散布図を描く関数を定義
54行目、2次元配列のAxesを1次元配列として扱っている

import numpy as np
import matplotlib.pyplot as plt
from mglearn.datasets import make_forge
from sklearn.neighbors import KNeighborsClassifier


def draw_decision_boundary(ax, n_neighbors, X, y, X0_field, X1_field):
    clsfr = KNeighborsClassifier(n_neighbors=n_neighbors)

    clsfr.fit(X, y)

    y_predicted = np.empty((len(X1_field), len(X0_field)))

    for row, x1 in enumerate(X1_field):
        for col, x0 in enumerate(X0_field):
            y_predicted[row, col] = clsfr.predict(np.array([[x0, x1]]))

    ax.contourf(X0_field, X1_field, y_predicted, levels=1, alpha=0.5)


def draw_scatter(ax, X0, X1, xlim, ylim):
    ax.scatter(X0[y==0], X1[y==0], marker='o', s=40, label="class-0")
    ax.scatter(X0[y==1], X1[y==1], marker='^', s=40, label="class-1")

    ax.set_xlim(xlim[0], xlim[1])
    ax.set_ylim(ylim[0], ylim[1])

    ax.set_xlabel("feature 0")
    ax.set_ylabel("feature 1")

    ax.tick_params(labelbottom=False, labelleft=False)
    ax.tick_params(bottom=False, left=False)

    ax.legend(loc='lower right')


X, y = make_forge()

X0_scatter = X[:, 0]
X1_scatter = X[:, 1]

n_X0_field, n_X1_field = 20, 20
y_predicted = np.empty((n_X1_field, n_X0_field))

xlim = (7.5, 12.5)
ylim = (-1.5, 6.5)
X0_field = np.linspace(xlim[0], xlim[1], n_X0_field)
X1_field = np.linspace(ylim[0], ylim[1], n_X1_field)

fig, axs = plt.subplots(2, 3, figsize=(9.6, 6.4))
fig.subplots_adjust(hspace= 0.4)

n_neighbors_list = [1, 2, 3, 24, 25, 26]
axs_1d = axs.reshape(1, -1)[0]

for n_neighbors, ax in zip(n_neighbors_list, axs_1d):
    ax.set_title("neighbors={}".format(n_neighbors))
    draw_decision_boundary(ax, n_neighbors, X, y, X0_field, X1_field)
    draw_scatter(ax, X0_scatter, X1_scatter, xlim, ylim)

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mglearn.datasets import make_forge

from sklearn.neighbors import KNeighborsClassifier

def draw_decision_boundary(ax, n_neighbors, X, y, X0_field, X1_field):

clsfr = KNeighborsClassifier(n_neighbors=n_neighbors)

clsfr.fit(X, y)

y_predicted = np.empty((len(X1_field), len(X0_field)))

for row, x1 in enumerate(X1_field):

for col, x0 in enumerate(X0_field):

y_predicted[row, col] = clsfr.predict(np.array([[x0, x1]]))

ax.contourf(X0_field, X1_field, y_predicted, levels=1, alpha=0.5)

def draw_scatter(ax, X0, X1, xlim, ylim):

ax.scatter(X0[y==0], X1[y==0], marker='o', s=40, label="class-0")

ax.scatter(X0[y==1], X1[y==1], marker='^', s=40, label="class-1")

ax.set_xlim(xlim[0], xlim[1])

ax.set_ylim(ylim[0], ylim[1])

ax.set_xlabel("feature 0")

ax.set_ylabel("feature 1")

ax.tick_params(labelbottom=False, labelleft=False)

ax.tick_params(bottom=False, left=False)

ax.legend(loc='lower right')

X, y = make_forge()

X0_scatter = X[:, 0]

X1_scatter = X[:, 1]

n_X0_field, n_X1_field = 20, 20

y_predicted = np.empty((n_X1_field, n_X0_field))

xlim = (7.5, 12.5)

ylim = (-1.5, 6.5)

X0_field = np.linspace(xlim[0], xlim[1], n_X0_field)

X1_field = np.linspace(ylim[0], ylim[1], n_X1_field)

fig, axs = plt.subplots(2, 3, figsize=(9.6, 6.4))

fig.subplots_adjust(hspace= 0.4)

n_neighbors_list = [1, 2, 3, 24, 25, 26]

axs_1d = axs.reshape(1, -1)[0]

for n_neighbors, ax in zip(n_neighbors_list, axs_1d):

ax.set_title("neighbors={}".format(n_neighbors))

draw_decision_boundary(ax, n_neighbors, X, y, X0_field, X1_field)

draw_scatter(ax, X0_scatter, X1_scatter, xlim, ylim)

plt.show()

疑 問

考え方

解 釈

ξやμはどのように定まるのか

参考サイト

定義

定式化

行列による表示

概要

学習率曲線

概要

適用対象となる問題

線形回帰による場合

確率分布の仮定

Logistic回帰の定式化

確率分布関数の仮定

最尤推定

1変数の場合

定式化

LogisticRegressionモデル

正則化

2変数の場合

特徴量の係数

例題のデータ

計算表の準備

収束計算

概要

各特徴量の分布

各特徴量と価格の関係

2つの特徴量と価格の関係

概要

データの取得とデータ構造

データの内容

'data'～特徴量データセット

'target'～住宅価格

'feature_names'～特徴名

'filename'～ファイル名

'DESCR'～データセットの説明

データの利用

データの取得方法

全レコードの特徴量データの取得

特定の特徴量のデータのみ取得

概要

近傍点数とクラス分類の挙動

近傍点数=1の場合

近傍点数=2の場合

近傍点数=3の場合

実行コード

knnの精度

予測カーブ

訓練データが少ない場合

訓練データが多い場合

概要

近傍点数によるクラス分類の挙動

近傍点数=1の場合

近傍点数=3の場合

近傍点数=2の場合

決定境界

疑　問

解　釈

`'data'`～特徴量データセット

`'target'`～住宅価格

`'feature_names'`～特徴名

`'filename'`～ファイル名

`'DESCR'`～データセットの説明