pythonでの多項式回帰の計算方法、ライブラリの使い方を調べてみました。
調べた結果を記載します。


多項式回帰とは

線形回帰をするには無理があるデータの場合、2次式、3次式、もしくはそれ以上で回帰をするのが自然な場合があります。この場合、多項式回帰 (polynomial regression) を使います。
多項式分析は、1変数のみの場合、特殊な重回帰分析です。
多項式回帰で検索するのと、フィッティングいう言葉で検索すると色々情報が出てきます。


多項式回帰が計算できるライブラリ

実際に計算する

統計局ホームページ/日本の統計 2018-第2章 人口・世帯2- 1 <wbr>人口の<wbr>推移と<wbr>将来人口(エクセル: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])

png

4次式にすると、大分変化に追従します。

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>]

png

numpy.polyfit も、最小二乗法を使っていますが、得られる結果が若干異なります。
poly_lsq何か補正をかけている? のかもしれません。

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()

png

scipy.odr を使った場合の結果とほぼ同じになりました。


参考

以下、参考になりました。
* numpy - polynomial regression using python - Stack Overflow
* 多項式回帰入門。線形回帰に飽きたらない人へ | 西陣に住むデータ分析屋のブログ * 多項式補間 - Wikipedia
* 多項式回帰と重回帰のどちらを使えばいいか判断基準はありますか【メンターが回答】 | TechAcademyマガジン * Scipy: curve fitting · tiagopereira/python_tips Wiki
* 最小二乗法でカーブフィッティング。関数3つを使い比べ-python | コード7区
* Tech Tips: scikit-learnで線形回帰

以上です。

コメント