最小二乗法（実装編） - 波動と機械学習

今回は最小二乗法（理論編）で解説したものを実際にプログラミングによって実装していきたいと思います。

そのための開発環境は以下の通りです。

Python

jupyter notebook

Pythonはインタプリタ型のプログラミング言語で、可読性の高さや機械学習・解析学などに必要となる演算のモジュールが多数用意されていることから、機械学習では多くの人がこのプログラミング言語を使って開発を行なっている。

jupyter notebookはPythonを用いた開発をサポートするオープンソースの開発環境である。
コードを「セル」という小さい単位から実行することが可能で、可読性の高いことと無駄のない実行が可能であることが特徴である。

f:id:araitbs007:20190217062603j:plain — 実際に最小二乗法をjupyter notebookで記述したもの。

※以下にコピペできるようにサンプルコードを記載していきますが、それぞれセルごとにまとめて処理を記載します。
（要は、私と同じ開発環境jupyter notebookで真似しながら勉強する方はそのままコピペして実行すればおそらく動くということ。）

実際のコード

取り扱う問題についても前回：最小二乗法（理論編）と同じく身長と体重の問題です。

データと関数の準備

1.)まずは今回使用するモジュールをある程度まとめてimport

import matplotlib.pyplot as plt
import numpy as np
import math
import random
%matplotlib inline

※もしここでエラーが起きる場合は、パソコンの中にモジュールが導入されていないかもしれないので以下のコードを試してみてください。

pip install (ここに必要なモジュール名)

2.)身長から理想体重を計算する関数の定義↓
（私たちは体重が身長の２乗に比例することを知っている状態でスタート）

#身長から算出する理想体重
def BMI(x):
    y=x*x*22
    return y

3.)次にデータを勝手に作る。
世の中の人々が全員理想体重ということはないでしょうからデータにある程度乱数を加えます。↓
ちなみに乱数は正規乱数（正規分布に従う乱数）を利用しています。
正規分布についての解説は別の記事にて解説する予定です。

#データの作成（x：身長、y：体重）
x=np.array([])
y=np.array([])
for i in range(20):#データ点の数が20個用意される。
    tall=np.random.normal(1.75 , 0.1)
    x=np.insert(x,len(x),tall,axis=0)
    norm=np.random.normal(0 , 2)
    weight=BMI(x[i])+norm
    y=np.insert(y,len(y),weight,axis=0)

print(x)
print(y)

4.)データを確認

#データの可視化
plt.scatter(x[:], y[:],color="red")
plt.show()

f:id:araitbs007:20190217061632j:plain

5.)推定する直線の関数↓
ここからが本題。データに対して「それっぽい直線」を考えていきます。
（a,b）が分かれば今回は目標達成になります。

#推定するもっともらしい直線（a,b:推定パラメータ）
def f(a,b,x):
    y=a*x+b
    return y

6.)平均値を求める関数↓
今後分散や共分散を求めるのに頻繁に使うことになる関数です。

#平均値の導出
def means(x):
    total=0
    for i in range(len(x)):
        total+=x[i]
    means=total/len(x)
    return means

パラメータ（a,b）の導出

パラメータの導出法は２通り試していきます。

１つ目はXの分散とX,Yの共分散から求める方法です。

8.1)定義通りに数式を組み立てる↓

#パラメータa,bはxの分散とx、yの共分散によって導出可

#平均値
x_means =means(x)
y_means =means(y)

#偏差
x_dev=x[:]-x_means
y_dev=y[:]-y_means

#分散（コードの内容はXとXの内積計算と等価）
x_s2=x_dev.dot(x_dev)/len(x)
xy_s=x_dev.dot(y_dev)/len(x)

print(x_means)
print(y_means)
print(x_s2)
print(xy_s)

8.2)こちらでも分散と共分散は出せる↓

x_s2=x.dot(x)/len(x)-x_means**2
xy_s=x.dot(y)/len(x)-x_means*y_means
print(x_s2)
print(xy_s)

9.)パラメータの導出が終了。そして表示。

#直線の傾きと切片
a=xy_s/x_s2
b=y_means-a*x_means
print(a)
print(b)

10.)推定した直線をデータ点と共に表示

#新しい傾きと切片で図示
x_axis=np.arange(1.5,2.0,0.01)
plt.scatter(x[:], y[:],color="red")
plt.plot(x_axis[:],a*x_axis[:]+b)
plt.show()

f:id:araitbs007:20190217061637j:plain

２つ目は偏微分を利用した方法です。
解析的には先のような方法は邪道です。パラメータがどんな値をとるかは未知の場合が多いで、どんなデータがきてもいいような仕組みにしていきます。

11.)sympyモジュールの用意
sympyは微分・積分などの演算を行うためのモジュールです。
他のプログラミング言語では、極限値の導出についてモジュールがなく自分で定義・記述しなければならない場合がありますが、このモジュールは人間が計算するのと同様にシンボルを使って文字式の微分・積分や極限の計算ができます。

#モジュールの導入とシンボル（a,b）の定義
from sympy import *
symbol_a = Symbol('a')
symbol_b = Symbol('b')

12.)最小二乗法の定義通り、誤差の二乗和をとる。

def min2error(a,b):
    total=0
    for i in range(len(x)):
        error=y[i]-f(a,b,x[i])
        total+=error**2
    return total

13.)a,bそれぞれの偏微分・連立

#偏微分
print(diff(min2error(symbol_a,symbol_b),symbol_a))
print(diff(min2error(symbol_a,symbol_b),symbol_b))
#連立方程式を解く
p=solve([diff(min2error(symbol_a,symbol_b),symbol_a), diff(min2error(symbol_a,symbol_b),symbol_b)], [symbol_a, symbol_b])

14.)結果の表示
ここで１つ目の方法で求めたパラメータの値と同じになっていれば成功です。
そして、ここで作成したアルゴリズムはデータを差し替えればいろんなデータに対して最小二乗法が使えるようになります。

a=p[symbol_a]
b=p[symbol_b]
print(a)
print(b)

#新しい傾きと切片で図示
x_axis=np.arange(1.5,2.0,0.01)
plt.scatter(x[:], y[:],color="red")
plt.plot(x_axis[:],a*x_axis[:]+b)
plt.show()

f:id:araitbs007:20190217061637j:plain — 大丈夫そう…

いかがだったでしょうか？
最小二乗法というパラメータ導出のための「アイデア」を形にしてみました。
自分の立てた理論を元にモノ作りをする楽しさが少しは味わえたのではないでしょうか？
では、次回：勾配降下法でまたお会いしましょう。