%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def f(x, a0, a1, a2):
    return a0 + a1*x + a2*x**2

xdata = np.array([1,2,3,4])
ydata = np.array([1,3,4,10])
plt.plot(xdata,ydata, 'o', color='r')

params, cov = curve_fit(f, xdata, ydata)
print(params)

x =np.linspace(0,4,20)
y = f(x,params[0],params[1],params[2])
plt.plot(x,y, color='b')

plt.grid()
plt.show()

[ 2.5 -2.2  1. ]

2次元の最小二乗フィット¶

以下のデータを

$$ f(x, y) = a_0+a_1 x+a_2 y+a_3 xy $$

で近似せよ

x,  y,  z
-1, -1,  2.00000
-1,  0,  0.50000
-1,  1, -1.00000
 0, -1,  0.50000
 0,  0,  1.00000
 0,  1,  1.50000
 1, -1, -1.00000
 1,  0,  1.50000
 1,  1,  4.00000

In [8]:

import numpy as np
z = np.array([2.00000,
0.50000,
-1.00000,
0.50000,
1.00000,
1.50000,
-1.00000,
1.50000,
4.00000])
x = []
y = []
for i in range(-1,2):
    for j in range(-1,2):
        x.append(i)
        y.append(j)

In [9]:

%matplotlib notebook
# 保存のためのコマンド．kernel restartがいるかも．
# import matplotlib
# matplotlib.use('Agg')

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(np.array(x),np.array(y),z) 

plt.show()

In [10]:

%matplotlib notebook

from pprint import pprint
import scipy.linalg as linalg

n = z.size
n_j = 4
bb=np.zeros([n])
A=np.zeros([n,n_j])
for i in range(0,n):
    A[i,0]=1
    A[i,1]=x[i]
    A[i,2]=y[i]
    A[i,3]=x[i]*y[i]
    bb[i]=z[i]

c, resid, rank, sigma = linalg.lstsq(A, bb)
pprint(c)

Ai = linalg.inv(np.dot(np.transpose(A),A))
b = np.dot(np.transpose(A),bb)
np.dot(Ai,b)

array([1. , 0.5, 0.5, 2. ])

Out[10]:

array([1. , 0.5, 0.5, 2. ])

In [11]:

%matplotlib notebook

def z_surf(xx,yy):
    val = c[0] + c[1]*xx + c[2]*yy
    val += c[3]*xx*yy
    return val

x1 = np.arange(-1, 1.25, 0.25)
y1 = np.arange(-1, 1.25, 0.25)
X, Y = np.meshgrid(x1, y1)
Z1 = z_surf(X,Y)

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(np.array(x),np.array(y),z, color='r') 
ax.plot_wireframe(X,Y,Z1) 

plt.show()

In [6]:

# plt.savefig("result.png")

In [ ]:

Table of Contents

1次元の線形最小二乗法¶

2次元の最小二乗フィット¶