よちよちpython

独習 python/Qpython/Pydroid3/termux

Numpyだけで回帰分析その3。poly1d()について

今回は、これまでNumpyだけで回帰分析(単回帰)するときに使ってきたpolyfit()poly1d()のうちのpoly1d()にしぼって軽くみていく。

なお、回帰分析の数学的説明、poly1d()のプログラムのアルゴリズム的な説明等は一切ありません。全然理解していません。



参考
numpy.poly1d — NumPy v1.17 Manual



実行環境


Androidスマホ
termux
Python3.8
JupyterNotebook

  • Pythonライブラリ
    • Numpy
    • matplotlib

np.poly1d() 多項式生成マシーン


Numpyのpoly1d(係数)は、

f(x) = x3-5x2+4x-9


のような変数Xのn次多項式と呼ばれる式を作ってくれます。

名前に1dとあることから察せられるように、1つのディメンション、独立変数が1個の式を扱う。変数はXと表示される。
独立変数(説明変数)が1個なので、これで回帰分析をできるのは単回帰分析に限ります。



この関数の使い方は、上の方程式を得たいなら次数の高い順から係数のリストをpoly1d([1,-5,4,-9])と引数で渡す。ない時は0で飛ばす。
戻り値で上のような多項式が表示として返ってくると同時に、その多項式が関数としても機能する。値を渡せば計算してくれる。

Numpyだけでの単回帰分析にpoly1d()polyfit()をセットで使う場合、polyfit()の戻り値poly1d()に渡すようにするので、手作業で係数を指定したり触ったりすることは特にないかと思う。知らんけど。



参考
多項式 Wikipedia

多項式と方程式と関数 : 算数・数学を10倍楽しく学習させる方法



poly1d()の使い方1


polyfit()とセットで単回帰分析する際の使われ方を実際に動かして見ていきます。

import numpy as np

# 多項式の係数を引数にセット
function = np.poly1d([1,-5,4,-9])

print(function)
   3     2
1 x - 5 x + 4 x - 9

マシーンがXの多項式を出力した。
xの右上にある「3」や「2」は「xの3乗」「xの2乗」の意味です。
表示がアレですが上のものと同じと思って話を続ける。

係数のリストを増やしてみる。

import numpy as np

function = np.poly1d([1,0,0,1,-5,4,-9])

print(function)
   6     3     2
1 x + 1 x - 5 x + 4 x - 9

poly1d(係数)の引数に渡すリストに左3つ分(1,0,0)を付け加えた。
6次の項が追加されている。



できた方程式を利用して次のようなことができる。

import numpy as np

function = np.poly1d([1,-5,4,-9])

# xを用意する
x=1

# 上で得た関数にxを代入しyを得る
y = function(x)

# できた多項式を表示
print(function)

# ワイは何や
print(y)
   3     2
1 x - 5 x + 4 x - 9
-9



xを配列で与えると?

import numpy as np

# 多項式生成
function = np.poly1d([1,-5,4,-9])

# xを配列で
x = [0,1,2]

# 関数に代入してyを求める
y = function(x)

# 結果
print(function)
print(x)
print(y)
   3     2
1 x - 5 x + 4 x - 9
[0, 1, 2]
[ -9  -9 -13]

適当にxの配列を作り直して、できあがった関数に代入し、グラフを描いてみましょう。

import numpy as np
import matplotlib.pyplot as plt

# xの多項式作成と表示
func = np.poly1d([-5,11,-8,-2])
print(func)

# 適当に配列xを作成
x = np.arange(-50,50)

# 関数に配列xを代入し配列yを求める
y = func(x)

# 配列の表示(出力省略)
#print(x)
#print(y)

# グラフ描画
plt.scatter(x,y)
# グラフ保存と表示
plt.savefig("poly1d_1.png")
plt.show()
    3      2
-5 x + 11 x - 8 x - 2

20200121233744

poly1d()のその他の使い方


上のやり方で得られた多項式の関数が単回帰分析で出来上がったモデルで、モデルに新たな数値を代入して得られる値予測値
という感じで、回帰分析的なpoly1d()の使い方は以上だろうと思いますが、scipyオフィシャルドキュメントページのpoly1d()の欄には他の使い方も載っているので以下で一応触れておきます。



引数パラメータ r


  • poly1d()の引数「r」 : 「r=False」がデフォルト

r=Trueにすると、今までやっていたような引数に係数を渡すやり方ではなく、引数を根とする方程式ができあがる。

生成される多項式は、

  1. r=False デフォルトのとき
    • poly1d([a, b, c]) = aX**2 + bX + c

  2. r=True のとき
    • poly1d([a, b, c]) = (X-a)(X-b)(X-c)

引数は「多項式の各項の係数」→「方程式の根」と、全然意味が変わる。

上でグラフを表示したコードのpoly1d()の引数に「r=True」だけ加えて実行してみます。

import numpy as np
import matplotlib.pyplot as plt

# xの多項式作成と表示
func = np.poly1d([-5,11,-8,-2], r=True)
print(func)

# 適当に配列xを作成
x = np.arange(-50,50)

# 関数に配列xを代入し配列yを求める
y = func(x)

# 配列の表示(出力省略)
#print(x)
#print(y)

# グラフ描画
plt.scatter(x,y)
# グラフ保存と表示
plt.savefig("poly1d_2.png")
plt.show()
   4     3      2
1 x + 4 x - 99 x - 646 x - 880

20200121233825

引数パラメータ variable


多項式をprintさせるとき変数名はデフォルトで「x」になっていましたが、それを変更できる。

# 変数名の確認
func.variable
'x'
import numpy as np

# xの多項式作成と表示
func2 = np.poly1d([1 ,-5 , 6], variable="β")
print(func2)
   2
1 β - 5 β + 6

多項式の根を求める


多項式=0を満たす方程式の解の事を「根」という(らしい)。
多項式」と「方程式」と「関数」の違いなどあまり気にしないし、「解」と「根」の違いも勿論。まぁそもそも会話に出てくる単語ではないしw

根は、上の式なら因数分解して(β-2)(β-3)=0を満たすβすなわち2と3のこと。

コンピュータはどうやって解くんだろう?解の公式?行列かなんか使ってやるのかな?

func2.r
array([3., 2.])

多項式に数値を代入


f = np.poly1d()で求まった多項式fをf(x)と書くと数式で関数化されるので、 そこにxを代入すると計算結果が得られる。

func2(7)
20

係数を表示する


自分で与えたものを自分で確認…

func2.c
array([ 1, -5,  6])

次数の表示


何次の方程式か

func2.order
2

その他省略。
参照

numpy.poly1d — NumPy v1.17 Manual



おわりに


poly1d()に変数Xのn次多項式を係数を与える事で出してもらい、それを利用してグラフを描画したりして遊んで参りましたが、

「なぜ回帰分析にこのn次多項式と呼ばれている式が関わって来るのか」

については、次回以降に回します。
1800年頃になってようやくあの天才数学者ガウスらが発見したようなぐあいのようですし、だいぶん難しいように思う。

線形回帰の説明で最小2乗法ってのが出てくる。データの各点を通過する直線y=ax+bがあると仮定し、その仮想的直線と各点の距離が小さくなるような直線の係数a,bを求める計算の過程で、a,bを2変数とする式ができる。で、その式を偏微分したものが0を満たすa,bを連立方程式で解けば…

みたいな解説の何処にもxの多項式の導出っぽい箇所がない。線形回帰は非線形回帰分析の特殊解なんだろう、polyfit()で与えるパラメータ毎にアルゴリズムを変えるような面倒臭いことにはなっていないはず。知らんけど。
このへんでそろそろおわるかの。


このつづきは、poly1d()の引数になっているpolyfit()をみてみることにする。

参考
多項式回帰 Wikipedia