Numpy.polynomialで回帰分析

#はじめに
PythonのNumpyパッケージ np.polyfitやnp.poly1dを利用するとn次式で 2 変数の回帰分析を行うことができます。
Numpy Version 1.14より、NumPy polynomialパッケージが導入され、従来のnp.polyfit や np.poly1dについては互換性を維持するされるそうですが、オープンシステムの特性上、ある日突然サポートが終了することがあるため、新しいAPIが出てきた以上、新規作成または長期メンテ・システムについては新しいAPIに対応させないと厳しいかなと感じたこと。
また、NumPy polynomialパッケージの記事がQiitaにもあったとうろ覚え記憶はあったのですが、あたらめて調べると見つけられませんでした。
Numpy Version 1.14で導入されたNumpy.polynomialサブパッケージでの回帰分析について下記２ファンクションについて書こうと思います。自分ごときが書くのもおこがましくはありますが。

1. numpy.polynomial.polynomial.polyfit()
2. numpy.polynomial.polynomial.polyval()

関数の詳しい仕様は上記リンクをたどってみて下さい。また回帰分析についての説明については、他に沢山の説明がありますので、ぜひググってみて下さい。

#回帰分析
今回はテストデータを3次曲線とし、下記式に近似させます。

y = a+bx+cx^2+dx^3

##polynomial　API

mt = np.polynomial.polynomial.polyfit(x, y, 3)
prod_y = lambda x, m:m[0]+m[1]*x+m[2]*x**2+m[3]*x**3)(x, mt)

3は曲線の次数です。線形の場合は1を設定します。
prod_y は回帰分析で求められた推定値です。mtより求められます。
注意事項としてnp.polyfit()とpolynomial.polyfit()では返却リストの並びが逆になっています。
3次であれば下記のように求められます。

prod_y = mt[0]+mt[1]x+mt[2]x^2+mt[3]x^3

###polynomial.polyvalで求める

mt = np.polynomial.polynomial.polyfit(x, y, 3)
prod_y = np.polynomial.polynomial.polyval(x, mt)

polynomial.polyval()を使う事でも求められます。
pythonですべて処理する場合は、polyval()を使えば済みますが、例えば計算結果をPLCなどに持って行って制御につかうなどの場合には、polyfit()返却値の内容を知っていても良いかなと思います。

##Pythonサンプルコード

import numpy as np
import matplotlib.pyplot as plt

# データ
x = np.array([0.0, 1.06666667, 2.13333333, 3.2, 4.26666667, 5.33333333,
     6.4, 7.46666667, 8.53333333, 9.6,10.66666667,11.73333333,
    12.8,13.86666667,14.93333333,16.0 ])
y = np.array([12.19088724, 9.89814203, 7.7212389 , 6.45078335, 6.922034  , 5.63559657,
      7.16868661, 7.1853348 , 7.50113099, 8.84172386,10.03269881,11.03363723,
     14.14321256,16.35153431,17.53782209,23.08017944,])

# 生データをグラフプロット
plt.scatter(x, y,marker='.', label='データ')

# 旧API 従来方式による回帰分析
ef = np.polyfit(x, y, 3)
print('旧API numpy.polyfit:', ef)

# 新API polynomialによる回帰分析
mt = np.polynomial.polynomial.polyfit(x, y, 3)
print('新API polynomial.polyfit:', mt)

# グラフプロット
plt.plot(x, (lambda x, m:m[0]+m[1]*x+m[2]*x**2+m[3]*x**3)(x, mt), label='lambdaでプロット')
plt.plot(x, np.polynomial.polynomial.polyval(x, mt), linestyle="--", label='polynomial.polyval')
plt.legend(prop={'family':'MS Gothic'})
plt.grid()
plt.show()

データは3次曲線を想定してフィッティングを行います。

##結果グラフ

同じフィッティング結果で計算項を直接使って計算しプロットした場合と、polyval()でプロットした場合の２つの方法で、２本グラフを描いています。
グラフは重なっているので、同じ結果になっているかと思います。
いいかんじにフィッティングしていると思います。

#感想
実績重視でなかなか新しいものを使いたがらないのは技術を扱う人はその傾向があるかと思いますが、長年オープンソースに接して、ある日突然のAPIの仕様変更、動きの仕様変更に、「やられた」と思う事たびたび。バージョンアップに恐怖すら感じる場合もあります。
幸いPythonは10年の歳月で2から3へのアップを行っていたので割と安心はしていますが、パッケージ類については油断できません。
後、オープンソース系の情報は原本、公式Webページ、開発コミュニティに当たらないとダメですね。

#環境
Windows Anaconda (Miniforge3)
Python 3.9.4 packaged by conda-forge
Numpy 1.20.2
matplotlib 3.4.2

#参考
Emotion Explorer - Python numpy.polynomial で最小二乗法
Emotion Explorer - Python numpy.polynomial で最小二乗法(2) 線形の単回帰
Transition notice - Polynomials
numpy.polynomial.polynomial.polyfit()
numpy.polynomial.polynomial.polyval()