#mynavi(機械学習の為の数学の基礎)
#setlinebreak(on);
#html(<script async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2/MathJax.js?config=TeX-MML-AM_CHTML"></script>)


* 目次 [#d8db4f2c]
#contents
- 関連
-- [[機械学習の為の数学の基礎]]
-- [[Python覚え書き]]
-- [[IPythonの使い方メモ]]
-- [[numpy入門]]
-- [[pandas入門]]
-- [[matplotlibでグラフ作成]]
-- [[scikit-learnを使用した重回帰分析]]
-- [[Chainer入門]]
-- [[Chainerで重回帰分析をやってみる]]

* 概要 [#g644cbdc]
#html(<div style="padding-left: 10px;">)
ここでは、最急降下法を使用して重回帰分析を出来るだけライブラリ等を使用せずにスクラッチで実装してみる。
#html(</div>)

* モデル式 [#ff473705]
#html(<div style="padding-left: 10px;">)

#html(){{
単回帰分析のモデル式が 
<math>
  <mover><mo>y</mo><mo>^</mo></mover>
  <mo>=</mo>
  <msub><mo>w</mo><mi>1</mi></msub>
  <mo>x</mo>
  <mo>+</mo>
  <msub><mo>w</mo><mi>0</mi></msub>
</math>
 であるのに対して、
}}

#html(){{
<div>
特徴量が複数ある重回帰分析の場合のモデル式は 
<math>
  <mover><mo>y</mo><mo>^</mo></mover>
  <mo>=</mo>
  <msub><mo>w</mo><mi>1</mi></msub>
  <msub><mo>x</mo><mi>1</mi></msub>
  <mo>+</mo>
  <msub><mo>w</mo><mi>2</mi></msub>
  <msub><mo>x</mo><mi>2</mi></msub>
  <mo>+</mo>
  <msub><mo>w</mo><mi>0</mi></msub>
  <msub><mo>x</mo><mi>0</mi></msub>
</math>
 となる。( 行列計算用に <math><msub><mo>x</mo><mi>0</mi></msub><mo>=</mo><mi>1</mi></math> とする )
</div>
}}

以下のようなデータがある時、
| x1 | x2 | y |h
| 70 | 68 | 87 |
| 75 | 82 | 86 |
| 80 | 71 |  93 |

#html(){{
<div style="padding: 10px 0;">
1行目は <math>
  <mover><mo>y</mo><mo>^</mo></mover>
  <mo>=</mo>
  <msub><mo>w</mo><mi>1</mi></msub>
  <mo>*</mo>
  <mi>70</mi>
  <mo>+</mo>
  <msub><mo>w</mo><mi>2</mi></msub>
  <mo>*</mo>
  <mi>68</mi>
  <mo>+</mo>
  <msub><mo>w</mo><mi>0</mi></msub>
  <mo>*</mo>
  <mi>1</mi>
</math> となり、
</div>

以下の行列計算で <math><mover><mo>y</mo><mo>^</mo></mover></math> を求める事ができる。
<div style="padding: 10px 0;">
<math>
    <mfenced open="[" close="]">
        <mtable>
            <mtr><mtd><mn>70</mn></mtd><mtd><mn>68</mn></mtd><mn>1</mn></mtd></mtr>
            <mtr><mtd><mn>75</mn></mtd><mtd><mn>82</mn></mtd><mn>1</mn></mtd></mtr>
            <mtr><mtd><mn>80</mn></mtd><mtd><mn>71</mn></mtd><mn>1</mn></mtd></mtr>
        </mtable>
    </mfenced>
    <mfenced open="[" close="]">
        <mtable>
            <mtr><mtd><mn><msub><mo>w</mo><mi>1</mi></msub></mn></mtd></mtr>
            <mtr><mtd><mn><msub><mo>w</mo><mi>2</mi></msub></mn></mtd></mtr>
            <mtr><mtd><mn><msub><mo>w</mo><mi>0</mi></msub></mn></mtd></mtr>
        </mtable>
    </mfenced>
    <mo>=</mo>
    <mfenced open="[" close="]">
        <mtable>
            <mtr><mtd><mi>70</mi><mo>*</mo><msub><mo>w</mo><mi>1</mi></msub><mo>+</mo><mi>68</mi><mo>*</mo><msub><mo>w</mo><mi>2</mi></msub><mo>+</mo><mi>1</mi><mo>*</mo><msub><mo>w</mo><mi>0</mi></msub></mtd></mtr>
            <mtr><mtd><mi>75</mi><mo>*</mo><msub><mo>w</mo><mi>1</mi></msub><mo>+</mo><mi>82</mi><mo>*</mo><msub><mo>w</mo><mi>2</mi></msub><mo>+</mo><mi>1</mi><mo>*</mo><msub><mo>w</mo><mi>0</mi></msub></mtd></mtr>
            <mtr><mtd><mi>80</mi><mo>*</mo><msub><mo>w</mo><mi>1</mi></msub><mo>+</mo><mi>71</mi><mo>*</mo><msub><mo>w</mo><mi>2</mi></msub><mo>+</mo><mi>1</mi><mo>*</mo><msub><mo>w</mo><mi>0</mi></msub></mtd></mtr>
        </mtable>
    </mfenced>
</math>
</div>

この為、特徴量の行列をX、パラメータの行列をWとした場合、以下のモデル式で表す事が出来る。
<div style="padding: 10px 0;">
<math>
  <mover><mo>y</mo><mo>^</mo></mover>
  <mo>=</mo>
  <mo>X</mo>
  <mo>W</mo>
</math>
</div>

}}

#html(</div>)
// モデル式

* コスト関数 [#u868c2e4]
#html(<div style="padding-left: 10px;">)

コスト関数は以下の通り。
#html(){{
<div>
<math>
    <mrow>
    <mo>j</mo><mo>(</mo><mo>W</mo><mo>)</mo>
    <mo>=</mo>
    </mrow>
    <mfrac><mi>1</mi><mo>2m</mo></mfrac>
    <munderover><mi>&sum;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover>
    <msup><mrow><mo>(</mo>  <mover><mi>y</mi><mo>^</mo></mover>  <mo>-</mo>  <msub><mi>y</mi><mi>i</mi></msub>  <mo>)</mo> </mrow><mi>2</mi></msup>
    <mo>=</mo>
    <mfrac><mi>1</mi><mo>2m</mo></mfrac>
    <msup>
      <mrow><mo>(</mo><mo>XW</mo><mo>-</mo><mo>y</mo><mo>)</mo></mrow>
      <mo>T</mo>
    </msup>
    <mrow><mo>(</mo><mo>XW</mo><mo>-</mo><mo>y</mo><mo>)</mo></mrow>
</math>
</div>
}}

#html(</div>)
// コスト関数

* 最急降下法 [#g528f800]
#html(<div style="padding-left: 10px;">)

最急降下法(ベクトル化)のモデル式は以下の通り。
#html(){{
<div>
<math>
  <mo>W</mo>
  <mo>:=</mo>
  <mo>W</mo>
  <mo>-</mo>
  <mo>α</mo>
  <mfrac><mo>1</mo><mo>m</mo></mfrac>
  <msup><mo>X</mo><mo>T</mo></msup>
  <mo>(</mo><mo>XW</mo><mo>-</mo><mo>y</mo><mo>)</mo>
</math>
</div>
}}

#html(</div>)
// 最急降下法

* 正規化 [#n73377f0]
#html(<div style="padding-left: 10px;">)

特徴量の範囲が列ごとに異なる場合、データを正規化する事によって、正しく分析を行う事ができる。

** z-socre Normalization (標準化) [#j3f49df1]
#html(<div style="padding-left: 10px;">)
平均が0、標準偏差が1となるように調整する。
#html(){{
<div style="padding-bottom: 10px;">
<math>
  <msub><mo>x</mo><mi>1</mi></msub>
  <mo>=</mo>
  <mfrac>
    <mrow><msub><mo>x</mo><mi>1</mi></msub><mo>-</mo><msub><mo>x</mo><mi>mean</mi></msub></mrow>
    <msub><mo>x</mo><mi>std</mi></msub>
  </mfrac>
</math>
</div>
}}
#html(</div>)
// z-socre Normalization

** min-max normalization [#o4c1e2e1]
#html(<div style="padding-left: 10px;">)
最大値が1最小値が0となるように調整する。

#html(){{
<div>
<math>
  <msub><mo>x</mo><mi>1</mi></msub>
  <mo>=</mo>
  <mfrac>
    <mrow><msub><mo>x</mo><mi>1</mi></msub><mo>-</mo><msub><mo>x</mo><mi>mean</mi></msub></mrow>
    <mrow><msub><mo>x</mo><mi>max</mi></msub><mo>-</mo><msub><mo>x</mo><mi>min</mi></msub></mrow>
  </mfrac>
</math>
</div>
}}

#html(</div>)
// min-max normalization

#html(</div>)
// 正規化

* サンプル実装 [#ddf2ff35]
#html(<div style="padding-left: 10px;">)

以下のデータから、任意の部屋面積と築年数の時に家賃がいくらになるかを、重回帰分析を使用して推論する処理を実装する。

|面積|築年数|家賃|h
|40.24|7|6.2|
|72.48|4|9.3|
|72.48|4|9.3|
|43.85|0|7.5|
|88.33|0|13.3|
※ &ref(sample_rent1.csv);

sample_rent1.py
#mycode2(python){{
"""重回帰分析サンプル."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def normalize(X, data=None):
    """正規化(Z-score normalization)."""
    data = X if data is None else data
    m = data.shape[0]
    X_norm = np.zeros((data.shape[0], data.shape[1]))
    for i in range(data.shape[1]):
        X_norm[:, i] = (data[:, i] - float(np.mean(X[:, i]))) / float(np.std(X[:, i]))

    # x0を追加
    X_norm = np.column_stack([np.ones([m,1]), X_norm])

    return X_norm

def cost(x, y, w):
    """コスト関数."""
    xw = x.dot(w)
    return np.dot((xw - y).T, (xw - y)) / (2*m)
    # もしくは
    #diff = np.power((x.dot(w) - y), 2)
    #return diff.sum(axis=0) / (2 * len(y))

def gradient_descent(x, y, w, alpha, iter_num):
    """最急降下法."""
    m = len(y)
    costs = np.zeros((iter_num, 1))
    for i in range(iter_num):
        w = w - alpha * (1.0/m) * np.transpose(x).dot(x.dot(w) - y)
        costs[i] = cost(x, y, w)
    return w, costs

if __name__ == "__main__":

    # --------------------------
    # 最急降下法による重回帰分析
    # --------------------------

    # データ読み込み
    data = np.loadtxt("data/sample_rent1.csv", delimiter=",", skiprows=1)
    x = data[:, 1:3]
    y = data[:, 3:4]

    # データの個数
    m = len(y)

    #  正規化し、x0を追加
    X_norm = normalize(x)

    # 初期値
    w_int = np.zeros((3, 1))

    # 学習率
    alpha = 0.01

    # 学習回数
    iter_num = 1000

    # 最急降下法による分析の実行
    w, costs = gradient_descent(X_norm, y, w_int, alpha, iter_num)

    # --------------------------
    # 作成したモデルを使用して別のデータを予測してみる
    # --------------------------
    z = np.array([[60, 10], [50, 10], [40, 10]])
    result = normalize(x, z).dot(w)
    for i in range(z.shape[0]):
        print("広さ: {}, 築年数: {}年 ≒ {:0.1f}万円".format(z[i,0], z[i,1], result[i,0]))

    # --------------------------
    # グラフ表示
    # --------------------------
    fig = plt.figure(figsize=(10, 5))

    # 3Dグラフ
    ax = fig.add_subplot(1, 2, 1, projection='3d')
    ax.scatter(data[:, 1], data[:, 2], data[:, 3], color="#ef1234")
    ax.set_xlabel("Size")
    ax.set_ylabel("Age")
    ax.set_zlabel("Rent")

    # 学習毎のコスト
    ax2 = fig.add_subplot(1, 2, 2)
    ax2.plot(range(costs.size), costs[:, 0], "r")
    ax2.set_xlabel("iterations")
    ax2.set_ylabel("cost")
    ax2.grid(True)
    plt.show()
}}

結果
#myterm2(){{
広さ: 60, 築年数: 10年 ≒ 8.1万円
広さ: 50, 築年数: 10年 ≒ 7.2万円
広さ: 40, 築年数: 10年 ≒ 6.3万円
}}

&ref(sample_rent1.png);


#html(</div>)

トップ   差分 バックアップ リロード   一覧 単語検索 最終更新   ヘルプ   最終更新のRSS