pythonで制限付きの多項式フィットを簡単にする方法

簡単な多項式フィット

多項式でデータをフィットしたいが、あるデータ点は必ず通過して欲しい場合に、簡単にやれる方法を紹介する。

コードの概要

numpyのpolyfitを使うときに、w というオプションで重みを指定できるので、固定したい点の重みを大きくする。ここでは、データは原点を通ってないが、フィットするときは原点を無理やり通したい場合の例を示します。元データ x, y に対して、numpy.appendでデータを追加して、その点の重みを1e3にして、他は1にしてあります。

サンプルコード

qiita_polyfit_limited.py

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt

x = np.array([-10, -8,-6,-4,-2,0,2,4,6,8,10])
y = np.array([-15,-10,-3,-2,-1,3,1,3,4,8,12])
weight = np.ones(len(x))
npoly=4

x_add = np.append(x,0)
y_add = np.append(y,0)
weight_add = np.append(weight,1e3)

coef=np.polyfit(x, y, npoly, w=weight)
yfit = np.poly1d(coef)(x)         

coef_add=np.polyfit(x_add, y_add, npoly, w=weight_add)
yfit_add = np.poly1d(coef_add)(x)         

F = plt.figure(figsize=(7,7))
ax = plt.subplot(1,1,1)

plt.errorbar(x,y,fmt="o",label="data",alpha=0.8,ms=3)
plt.errorbar(x,yfit,fmt="--",label="fit (no limit), order = "  +str(npoly),alpha=0.8,ms=3)
plt.errorbar(x,yfit_add,fmt="--",label="fit (fixed at 0,0), order = "  +str(npoly),alpha=0.8,ms=3)

plt.grid(linestyle='dotted',alpha=0.5)
plt.legend(numpoints=1, frameon=False, loc="best") 

outfile = "qiita_polyfitlimited.png"
plt.savefig(outfile)

フィット結果

何も制限をしないと原点を通過しないが、制限付きだと原点を通過しているのがわかる。

厳密な条件付きのフィットではないので、原点でゼロではなくて、小さい値が残る。もし、厳密に原点で0にしたい場合は、制限付きのフィットの関数を自分で書いて行う必要がある。その場合は、pythonでパラメータの制限付きでfittingを行う方法
を参照ください。