LoginSignup
13

More than 5 years have passed since last update.

[Python]多項式によるカーブフィッティング

Posted at

サンプルデータに対し多項式でカーブフィッティングしたいときはnumpyのpolyfitを使うと簡単にできる。曲線は下記の多項式で表される。polyfitは下の式の$a_n$を計算してくれる。

y = \sum^{N}_{n=0} a_n x^n

polyfitで返される値は次数が高い方ほうから並んでいる。数式にすると下記のような感じになる。

y = \sum^N_{n=0} a_n x^{N-n}

polyfitで取得した係数を使った曲線のグラフを書きたいときは,xの値をnumpy.linspace等で用意して、numpy.polyvalに係数とxを渡すとyを計算してくれる。

下図はサンプルデータと1次式でカーブフィッティングした結果
numpy.polyfit(x,y,1)とするだけで係数が得られる。x,yはデータで1は次数を表す。
download.png

w = np.polyfit(x,y,1)
xs = np.linspace(np.min(x),np.max(x),100)
ys = np.polyval(w,xs)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x, y,label='input')
ax.plot(xs, ys, 'r-', lw=2, label='polyfit')
ax.set_title('polyfit w = [%.2f, %.2f]' % (w[0],w[1]))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(loc='best',frameon=False)
ax.grid(True)

サンプルデータは下記のように作成した。

def getData1(numPoints, w,b, nstdev, xmin=0.0, xmax=1.0 ):
    x = scipy.stats.uniform.rvs(loc=xmin, scale=xmax-xmin,size=numPoints)
    n = scipy.stats.norm.rvs(size=numPoints, scale=nstdev)

    y = w * x + b + n
    return x, y

x, y = getData1(100, 1.5, 10, 10, xmin=0, xmax=100)

3次多項式でカーブフィッティングした結果
download.png

w = np.polyfit(x,y,3)

xs = np.linspace(np.min(x),np.max(x),100)
ys = np.polyval(w,xs)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x, y,label='input')
ax.plot(xs, ys, 'r-', lw=2, label='polyfit')
ax.set_title('polyfit w = [%.2f, %.2f, %.2f, %.2f]' % (w[0],w[1],w[2],w[3]))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.legend(loc='best',frameon=False)
ax.grid(True)

サンプルデータは下記のように作成した。

def getData2(numPoints, Amp,  freq, phase, nstdev, xmin=0.0, xmax=np.pi*2.0 ):
    x = scipy.stats.uniform.rvs(loc=xmin, scale=xmax-xmin,size=numPoints)
    n = scipy.stats.norm.rvs(size=numPoints, scale=nstdev)

    y = Amp*np.sin(freq * x + phase) + n
    return x, y

x, y = getData2(100, 10, 1.0, 0.0, 5, xmin=0, xmax=np.pi*2.0)

参考URL
https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html
https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html#numpy.polyval

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
13