pythonでの多項式回帰の計算方法、ライブラリの使い方を調べてみました。
調べた結果を記載します。
多項式回帰とは
線形回帰をするには無理があるデータの場合、2次式、3次式、もしくはそれ以上で回帰をするのが自然な場合があります。この場合、多項式回帰 (polynomial regression) を使います。
多項式分析は、1変数のみの場合、特殊な重回帰分析です。
多項式回帰で検索するのと、フィッティングという言葉で検索すると色々情報が出てきます。
多項式回帰が計算できるライブラリ
- numpy.polyfit — NumPy v1.15 Manual
- Orthogonal distance regression (scipy.odr) — SciPy v1.1.0 Reference Guide
scipy,odr には、これという機能はないですが、駆使できれば多項式回帰もできるようです。 - sklearn.preprocessing.PolynomialFeatures — scikit-learn 0.19.2 documentation
実際に計算する
統計局ホームページ/日本の統計 2018－第2章 人口・世帯 の2- 1 人口の推移と将来人口（エクセル：42KB） の総人口データを使用します。
* numpy.polyfitを使う
まず、総人口データの散布図を描きます。
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
x = np.array([1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2020,2025,2030,2035,2045,2055,2065,2075,2085,2095])
y = np.array([55963,59737,64450,69254,71933,72147,84115,90077,94302,99209,104665,111940,117060,121049,123611,125570,126926,127768,128033,128084,128032,128057,127834,127593,127414,127237,127095,126933,125325,122544,119125,115216,106421,97441,88077,78564,70381,63125])
plt.scatter(x , y)
# `np.polyfit` で、4次式で近似
plt.plot(x, np.poly1d(np.polyfit(x, y, 4))(x), label='d=4', color="green")
# `np.polyfit` で、2次式で近似
plt.plot(x, np.poly1d(np.polyfit(x, y, 2))(x), label='d=2', color="red")
array([ 9.61741508e-04, -7.71943642e+00, 2.32189763e+04,
-3.10181359e+07, 1.55281205e+10])
4次式にすると、大分変化に追従します。
- scipy.odr を使う
python_tips/poly_lsq.py at master · tiagopereira/python_tips を python3 向けに編集しました。scipy odrpackを使って、多項式最小二乗法をデータに適用します。
from scipy.odr import odrpack as odr
from scipy.odr import models
def poly_lsq(x,y,n,verbose=False,itmax=200):
''' Performs a polynomial least squares fit to the data,
with errors! Uses scipy odrpack, but for least squares.
IN:
x,y (arrays) - data to fit
n (int) - polinomial order
verbose - can be 0,1,2 for different levels of output
(False or True are the same as 0 or 1)
itmax (int) - optional maximum number of iterations
OUT:
coeff - polynomial coefficients, lowest order first
err - standard error (1-sigma) on the coefficients
--Tiago, 20071114
'''
# http://www.scipy.org/doc/api_docs/SciPy.odr.odrpack.html
# see models.py and use ready made models!!!!
func = models.polynomial(n)
mydata = odr.Data(x, y)
myodr = odr.ODR(mydata, func,maxit=itmax)
# Set type of fit to least-squares:
myodr.set_job(fit_type=2)
if verbose == 2: myodr.set_iprint(final=2)
fit = myodr.run()
# Display results:
if verbose: fit.pprint()
if fit.stopreason[0] == 'Iteration limit reached':
print('(WWW) poly_lsq: Iteration limit reached, result not reliable!')
# Results and errors
coeff = fit.beta[::-1]
err = fit.sd_beta[::-1]
return coeff,err
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
x = np.array([1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2020,2025,2030,2035,2045,2055,2065,2075,2085,2095])
y = np.array([55963,59737,64450,69254,71933,72147,84115,90077,94302,99209,104665,111940,117060,121049,123611,125570,126926,127768,128033,128084,128032,128057,127834,127593,127414,127237,127095,126933,125325,122544,119125,115216,106421,97441,88077,78564,70381,63125])
plt.scatter(x , y)
coeff, err = poly_lsq(x, y , 4)
print("polynomial coefficients:", coeff)
print("standard error:", err)
# グラフ描画
plt.plot(x, np.poly1d(coeff)(x), label='d=4', color="green")
polynomial coefficients: [ -2.42404329e-06 9.72801674e-03 -9.73287786e+00 7.75065922e+00
1.00000000e+00]
standard error: [ 8.35504229e-08 3.35038696e-04 3.35776615e-01 0.00000000e+00
0.00000000e+00]
[<matplotlib.lines.Line2D at 0x1063f9f28>]
numpy.polyfit も、最小二乗法を使っていますが、得られる結果が若干異なります。
poly_lsq が何か補正をかけている? のかもしれません。
- sklearn.preprocessing.PolynomialFeatures を使う
Tech Tips: scikit-learnで線形回帰 を参考に以下作成しました。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
# T でx、y を入れ替える
x = np.array([[1920,1925,1930,1935,1940,1945,1950,1955,1960,1965,1970,1975,1980,1985,1990,1995,2000,2005,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2020,2025,2030,2035,2045,2055,2065,2075,2085,2095]]).T
y = np.array([[55963,59737,64450,69254,71933,72147,84115,90077,94302,99209,104665,111940,117060,121049,123611,125570,126926,127768,128033,128084,128032,128057,127834,127593,127414,127237,127095,126933,125325,122544,119125,115216,106421,97441,88077,78564,70381,63125]]).T
# train a linear regression model
regr = Pipeline([
('poly', PolynomialFeatures(degree=4)),
('linear', LinearRegression())
])
regr.fit(x, y)
# make predictions
xt = np.linspace(1920, 2100, num=300, dtype = 'int').reshape(300, 1)
yt = regr.predict(xt)
# plot samples and regression result
plt.plot(x, y, 'o')
plt.plot(xt, yt)
plt.show()
scipy.odr を使った場合の結果とほぼ同じになりました。
参考
以下、参考になりました。
* numpy - polynomial regression using python - Stack Overflow
* 多項式回帰入門。線形回帰に飽きたらない人へ | 西陣に住むデータ分析屋のブログ
* 多項式補間 - Wikipedia
* 多項式回帰と重回帰のどちらを使えばいいか判断基準はありますか【メンターが回答】 | TechAcademyマガジン
* Scipy: curve fitting · tiagopereira/python_tips Wiki
* 最小二乗法でカーブフィッティング。関数３つを使い比べ－python | コード７区
* Tech Tips: scikit-learnで線形回帰
以上です。
コメント