Python numpyのpolyfitで二次元データを二次関数(y=ax^2+bx+c)でフィッティング、係数と誤差を得る

phst.hateblo.jp

の二次関数版

import numpy as np

# データ生成
x = np.array([10.09659323, 10.40523609, 10.72069954, 11.03063238, 11.33665493, 11.6391395,
              11.9381814,  12.23925842, 12.5349865,  12.82765475, 13.11407145, 13.40273189,
              13.69487096, 13.97883462, 14.28398633, 14.57440393, 14.84567687, 15.1130868,
              15.37613046, 15.63709515, 15.89563802, 16.15306505, 16.40719798, 16.66447298,
              16.91974601, 17.17561009, 17.4268742 , 17.67091169, 17.90822826, 18.14299743,
              18.36384415, 18.65457429, 18.87981717, 19.12584329, 19.35925588, 19.59113967,
              19.80996694, 20.0201713,  20.25120795, 20.50498234, 20.74189141, 20.9711738,
              21.18167175])
y = np.array([40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 
              50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 
              60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 
              70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 
              80, 81, 82])
yerr = np.array([0.07576136, 0.07513219, 0.07056666, 0.0683801 , 0.07147065, 0.0664667,
                 0.07014482, 0.06589157, 0.07052181, 0.06397122, 0.06861484, 0.0705902,
                 0.07197119, 0.07523888, 0.0782734 , 0.06497902, 0.07018011, 0.06930904,
                 0.07239904, 0.07280576, 0.07422179, 0.07614255, 0.07725881, 0.07820666,
                 0.07855506, 0.07713332, 0.08082234, 0.0814835 , 0.08801933, 0.08666913,
                 0.12859639, 0.10607327, 0.14610974, 0.09157884, 0.08456184, 0.08410188,
                 0.08778977, 0.10517566, 0.162944  , 0.07941121, 0.08081205, 0.09473661,
                 0.08497694])

# np.polyfitで係数と係数の誤差を計算, 2次関数なので3番目の引数は2, 誤差はw=1/sigmaで与える
p, cov =np.polyfit(x,y,2, w=1.0/yerr,cov=True)
a = p[0]
b = p[1]
c = p[2]
sigma_a = np.sqrt(cov[0,0])
sigma_b = np.sqrt(cov[1,1])
sigma_c = np.sqrt(cov[2,2])

# グラフ生成
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [10, 4]

xx=np.linspace(x[0],x[-1],100)
yy=np.polyval(p,xx) #結果からyを計算

plt.errorbar(x,y,yerr=yerr,fmt='o',label='data',capsize=5,ecolor='black',markeredgecolor="black",color='w')
plt.plot(xx,yy,label="y = ({:.3f} $\pm$ {:.3f}) x^2 ({:.3f} $\pm$ {:.3f}) x + ({:.3f} $\pm$ {:.3f})".format(a,sigma_a,b,sigma_b,c,sigma_c))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid()
plt.show()