読者です 読者をやめる 読者になる 読者になる

Remrinのpython攻略日記

python3に入門しました。python3についてあれこれとサンプルコードとか。

最小二乗法

最小二乗法について。
 
(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()

 
東京の過去140年の平均気温を最小二乗法で近似。
1次式で最小二乗法
 
f:id:rare_Remrin:20170518163251p:plain

2次式で最小二乗法
 
f:id:rare_Remrin:20170518163301p:plain
 
ただし、3次以上だとうまくいかないので、原因調査中
 
参考:
正規方程式の導出と計算例 | 高校数学の美しい物語