最小二乗法 - Remrinのpython攻略日記

最小二乗法について。
　
(2, 3), (4, 7), (9, 11)の3点をデータとして、最小二乗法を1次式で行う場合。
f:id:rare_Remrin:20170518140226p:plain
とすると、
f:id:rare_Remrin:20170518140256p:plain
これの最小値を求めることになり、偏微分をして
f:id:rare_Remrin:20170518140345p:plain
となるベクトルが係数ベクトルとなる。
　

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt

data = [(2, 3), (4, 7), (9, 11)]
a = np.matrix(data)
b = a[:, 1].copy()
a[:, 1] = 1

# 係数ベクトルを求める
x = ((a.T * a)**-1) * a.T * b
x= np.array(x)

# matplotlibで描画
datax, datay = np.split(np.array(data), 2, axis=1)
maxx, maxy = max(datax), max(datay)
plt.xlim(0, maxx*1.1)
plt.ylim(0, maxy*1.1)
plt.scatter(datax, datay)
plt.plot([0, maxx*1.1], [x[1], maxx*x[0]*1.1 + x[1]], color="r")
plt.text(1, maxy, "y={:.2f}x+{:.2f}".format(*x[:, 0]), fontsize=14)
plt.show()

　
f:id:rare_Remrin:20170518140520p:plain
　
2次式以上でも最小二乗法で近似できるようにプログラムを作ってみました。

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# 東京の過去140年間の平均気温のデータを取得
url = 'http://python-remrin.hatenadiary.jp/entry/2017/05/18/142423'
## DataFrameのリストを得る。header=0のオプション指定で、最初の行をheader扱い。
fetched = pd.io.html.read_html(url, header=0)
df = fetched[0]

# 2次元リスト化
data = []
for i in range(len(df)):
    x = []
    for j in range(len(df.iloc[0])):
        x.append(df.iloc[i, j])
    data.append(x)

data = np.array(data)
xlabel = data[0:, 0]
datasize = data.shape[0]
datax = np.arange(datasize)
datay = np.array(data[:, -1], dtype=float)

# n次式で近似するための行列を作る
n = 2
a = np.empty((n + 1, datasize), dtype=int)
for i in range(n + 1):
    a[i] = datax**(n - i)
a = np.matrix(a.T)
b = np.matrix(datay).T

## 係数ベクトルを求める
x = ((a.T * a)**-1) * a.T * b
x= np.array(x)

def f(a):
    fa = 0
    for i in range(n + 1):
        fa += x[i] * a**(n-i)   
    return fa

def textf():
    txt = "y="
    for i in range(n):
        txt +="{:+.3f}x^{}".format(x[i][0], n-i)
    txt += "{:+.1f}".format(x[n][0])
    return txt

print(x)
print(textf())

# matplotlibで描画
maxx, maxy = max(datax), max(datay)
plt.xlim(0, maxx*1)
plt.ylim(10, maxy*1.1)
plt.scatter(datax, datay)
plt.plot(datax, f(datax), color="r")                # 近似曲線、直線
plt.text(1, maxy, textf(), fontsize=14)             # 近似式
plt.xticks(np.arange(datasize)[::20], xlabel[::20]) # 横軸目盛りは20ごと
plt.title("mean temperature in Tokyo")
plt.show()