Last active
November 7, 2015 23:34
-
-
Save aidiary/7ca6c1cd25d9a1ad2345 to your computer and use it in GitHub Desktop.
多項式曲線フィッティング
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#coding:utf-8 | |
import numpy as np | |
import sys | |
from pylab import * | |
# M次多項式近似 | |
M = 3 | |
def y(x, wlist): | |
ret = wlist[0] | |
for i in range(1, M+1): | |
ret += wlist[i] * (x ** i) | |
return ret | |
# 訓練データからパラメータを推定 | |
def estimate(xlist, tlist): | |
# M次多項式のときはM+1個のパラメータがある | |
A = [] | |
for i in range(M+1): | |
for j in range(M+1): | |
temp = (xlist**(i+j)).sum() | |
A.append(temp) | |
A = array(A).reshape(M+1, M+1) | |
T = [] | |
for i in range(M+1): | |
T.append(((xlist**i) * tlist).sum()) | |
T = array(T) | |
# パラメータwは線形連立方程式の解 | |
wlist = np.linalg.solve(A, T) | |
return wlist | |
def example1(): | |
# 訓練データ | |
# sin(2*pi*x)の関数値にガウス分布に従う小さなランダムノイズを加える | |
xlist = np.linspace(0, 1, 10) | |
tlist = np.sin(2*np.pi*xlist) + np.random.normal(0, 0.2, xlist.size) | |
# 訓練データからパラメータを推定 | |
wlist = estimate(xlist, tlist) | |
print wlist | |
# 連続関数のプロット用X値 | |
xs = np.linspace(0, 1, 500) | |
ideal = np.sin(2*np.pi*xs) # 理想曲線 | |
model = [y(x, wlist) for x in xs] # 推定モデル | |
plot(xlist, tlist, 'bo') # 訓練データ | |
plot(xs, ideal, 'g-') # 理想曲線 | |
plot(xs, model, 'r-') # 推定モデル | |
xlim(0.0, 1.0) | |
ylim(-1.5, 1.5) | |
show() | |
if __name__ == "__main__": | |
example1() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 訓練データからパラメータを推定 | |
def estimate(xlist, tlist): | |
# M次多項式のときはM+1個のパラメータがある | |
A = [] | |
for i in range(M+1): | |
for j in range(M+1): | |
temp = (xlist**(i+j)).sum() | |
A.append(temp) | |
A = array(A).reshape(M+1, M+1) | |
T = [] | |
for i in range(M+1): | |
T.append(((xlist**i) * tlist).sum()) | |
T = array(T) | |
# パラメータwは線形連立方程式の解 | |
wlist = np.linalg.solve(A, T) | |
return wlist |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# M次多項式近似 | |
M = 3 | |
def y(x, wlist): | |
ret = wlist[0] | |
for i in range(1, M+1): | |
ret += wlist[i] * (x ** i) | |
return ret |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 訓練データからパラメータを推定 | |
def estimate(xlist, tlist, lam): | |
# M次多項式のときはM+1個のパラメータがある | |
A = [] | |
for i in range(M+1): | |
for j in range(M+1): | |
temp = (xlist**(i+j)).sum() | |
if i == j: temp += lam # 対角成分にλを足します | |
A.append(temp) | |
A = array(A).reshape(M+1, M+1) | |
... | |
def example1(): | |
... | |
wlist = estimate(xlist, tlist, np.exp(-18.0)) | |
... |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#coding:utf-8 | |
import numpy as np | |
from pylab import * | |
# 訓練データの数 | |
N = 10 | |
# N個の訓練データを生成 | |
xlist = np.linspace(0, 1, N) # 0から1の間から均等にN点を抽出 | |
tlist = np.sin(2 * np.pi * xlist) + np.random.normal(0, 0.2, xlist.size) # sinにN(0, 0.2)の乱数を加える | |
# オリジナルのsin用データ | |
xs = np.linspace(0, 1, 1000) # sinは連続関数なのでxを細かくとる | |
ideal = np.sin(2 * np.pi * xs) # xsに対応するsinを求める | |
# 訓練データとオリジナルのsinデータをプロット | |
plot(xlist, tlist, 'bo') # 訓練データを青い(b)の点(o)で描画 | |
plot(xs, ideal, 'g-') # sin曲線を緑(g)の線(-)で描画 | |
xlim(0.0, 1.0) # X軸の範囲 | |
ylim(-1.5, 1.5) # Y軸の範囲 | |
show() # グラフを表示 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment