Pythonで頂点を通る近似曲線の係数の求め方

概要

2次元平面上の$n$個の点 $(x_1,y_1),\cdots,(x_n,y_n)$ を通る$m$次曲線 $y=a_mx^m+a_{m-1}x^{m-1}+\cdots+a_1x+a_0$ の係数 $a_0,\dots,a_m$ を求めたい。
$n\le m+1 $ かつ $x_1,\cdots,x_n$ の値がすべて異なる場合すべての点を通る曲線が存在するが、そうでない場合は点と曲線の誤差が最小となるような係数を取る必要がある。
その際にscipyのfminメソッドが使えそうだ。
ただしfmin関数は目的関数の局所最適解を求めるため、最適解となっていない場合があるため注意。

3次曲線のプログラム例

import numpy as np
from scipy.optimize import fmin

def f(a, x):
    # 3次曲線
    return a[0] + a[1]*x + a[2]*x**2 + a[3]*x**3

def obj_func(a,x,y):
    # 頂点と曲線との二乗誤差の和を目的関数とする
    return sum((f(a,x) - y)**2)

# 5点を与える
x = np.array([ -8.,  -3.,  2.,   6.,   9.])
y = np.array([ 30., -22., 15., -17., -25.])

# 係数の初期値
a= np.array([ 0., 0., 0., 0.])

# 最適な係数を求める
opt = fmin(obj_func, a, args=(x,y))

print(opt)
# [-8.20775435  3.2028276   0.2150416  -0.09395258]