MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/
はじめに
の、
MATLAB,Python,Scilab,Julia比較 第2章 その74【多変量多項式回帰分析(関数項)③】
を書き直したもの。
正規方程式を用いた、多変量多項式回帰分析(関数項あり)について。
今回は、Python(NumPy)で演算してみる。
正規方程式、各パラメータ、推定対象の多項式再掲
まずは恒例の正規方程式、多変量多項式回帰分析(関数項あり)で想定するパラメータの再掲。
正規方程式
\(
x=(A^TA)^{-1}A^Tb
\)
多変量多項式回帰分析(関数項あり)に於ける各パラメータ
\(
A=
\begin{bmatrix}
x_1^2 & \cos(6x_1) & y_1^2 & \exp(2y_1) &1\\
x_2^2 & \cos(6x_2) & y_2^2 & \exp(2y_2) &1\\
\vdots & \vdots & \vdots\\
x_n^2 & \cos(6x_n) & y_n^2 & \exp(2y_n) &1\\
\end{bmatrix},
\vec{x}=
\begin{bmatrix}
\alpha\\
\beta\\
\gamma\\
\delta\\
\epsilon\\
\end{bmatrix},
\vec{b}=
\begin{bmatrix}
z_1\\
z_2\\
\vdots\\
z_n
\end{bmatrix}
\)
推定対象の多項式
\(
z=4x^2-5\cos(6x)+3y^2+\exp(2y)+2
\)
これをPython(NumPy)で実現する。
Pythonコード
Pythonコードは以下になる。
import numpy as np
import matplotlib.pyplot as plt
n = 100
x = np.random.rand(n, 1)
y = np.random.rand(n, 1)
z = 4*x**2 - 5*np.cos(6*x) + 3*y**2 + np.exp(2*y) + 2 + np.random.rand(n, 1) -0.5
A=np.block([x**2, np.cos(6*x), y**2, np.exp(2*y), np.ones((x.size,1))])
b=z
X=np.linalg.inv(A.T@A)@A.T@b
print(X)
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection='3d')
ax.scatter3D(x, y ,z)
xp=np.linspace(0, 1, 10)
yp=np.linspace(0, 1, 10)
xpm,ypm=np.meshgrid(xp,yp)
ax.plot_wireframe( xpm, ypm, X[0]*xpm**2 + X[1]*np.cos(6*xpm) + X[2]*ypm**2 + X[3]*np.exp(2*ypm) +X[4])
ax.view_init(elev=20, azim=230)
plt.show()
処理結果
処理結果は以下。
[[ 4.06591289]
[-5.05828082]
[ 2.26445071]
[ 1.12934576]
[ 1.80696875]]考察
やや誤差は出ている。
これもサンプル点数が少ないためと思われる。
10000点にすると以下になる。
[[ 4.0057331 ]
[-4.99665852]
[ 3.08997676]
[ 0.98286319]
[ 2.02269517]]まとめ
- 正規方程式による多変量多項式回帰分析(関数項あり)をPython(NumPy)で実施。
- 誤差はあるものの目的の係数の算出はできている。
- サンプル点数を増やせば、理想値に近付く。
MATLAB、Python、Scilab、Julia比較ページはこちら
コメント