目次

概要

ここでは、最急降下法を使用して重回帰分析を出来るだけライブラリ等を使用せずにスクラッチで実装してみる。

モデル式

単回帰分析のモデル式が  y^ = w1 x + w0  であるのに対して、
特徴量が複数ある重回帰分析の場合のモデル式は  y^ = w1 x1 + w2 x2 + w0 x0  となる。( 行列計算用に x0=1 とする )

以下のようなデータがある時、

x1x2y
706887
758286
807193
1行目は y^ = w1 * 70 + w2 * 68 + w0 * 1 となり、
以下の行列計算で y^ を求める事ができる。
70681 75821 80711 w1 w2 w0 = 70*w1+68*w2+1*w0 75*w1+82*w2+1*w0 80*w1+71*w2+1*w0
この為、特徴量の行列をX、パラメータの行列をWとした場合、以下のモデル式で表す事が出来る。
y^ = X W

コスト関数

コスト関数は以下の通り。

j(W) = 12m i=1m ( y^ - yi ) 2 = 12m (XW-y) T (XW-y)

最急降下法

最急降下法(ベクトル化)のモデル式は以下の通り。

W := W - α 1m XT (XW-y)

正規化

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

z-socre Normalization (標準化)

平均が0、標準偏差が1となるように調整する。

x1 = x1-xmean xstd

min-max normalization

最大値が1最小値が0となるように調整する。

x1 = x1-xmean xmax-xmin

サンプル実装

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

面積築年数家賃
40.2476.2
72.4849.3
72.4849.3
43.8507.5
88.33013.3

filesample_rent1.csv

sample_rent1.py

"""重回帰分析サンプル."""
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()

結果

広さ: 60, 築年数: 10年 ≒ 8.1万円
広さ: 50, 築年数: 10年 ≒ 7.2万円
広さ: 40, 築年数: 10年 ≒ 6.3万円

sample_rent1.png


添付ファイル: filesample_rent1.png 7件 [詳細] filesample_rent1.csv 9件 [詳細]

トップ   差分 バックアップ リロード   一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-11-22 (金) 07:49:12 (20d)