MATLAB,Python,Scilab,Julia比較 第2章 その76【多変量多項式回帰分析(関数項)⑤】

MATLAB,Python,Scilab,Julia比較 第2章 その76【多変量多項式回帰分析(関数項)⑤】 数値計算
MATLAB,Python,Scilab,Julia比較 第2章 その76【多変量多項式回帰分析(関数項)⑤】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia2-backnumber/

はじめに

正規方程式を用いた、多変量多項式回帰分析(関数項あり)について。
今回は、Juliaで演算してみる。

登場人物

博識フクロウのフクさん

指差しフクロウ

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

技術者太郎

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

正規方程式、各パラメータ、推定対象の多項式再掲

太郎くん
太郎くん

まずは恒例の正規方程式、多変量多項式回帰分析(関数項あり)で想定するパラメータの再掲。

正規方程式

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

フクさん
フクさん

これをJuliaで実現する。

Juliaコード

フクさん
フクさん

Juliaコードは以下になる。

using PyPlot

function meshgrid(xin,yin)
	nx=length(xin)
	ny=length(yin)
	xout=zeros(ny,nx)
	yout=zeros(ny,nx)
	for jx=1:nx
	    for ix=1:ny
	        xout[ix,jx]=xin[jx]
	        yout[ix,jx]=yin[ix]
	    end
	end
	return (x=xout, y=yout)
end

n = 100

x = rand(n, 1)
y = rand(n, 1)
z = 4*x.^2 - 5*cos.(6*x) .+ 3*y.^2 + exp.(2*y) .+ 2 + rand(n, 1) .-0.5

A=[x.^2  cos.(6*x)  y.^2  exp.(2*y)  ones(length(x),1)];
b=z
X=(A'*A)^-1 *A'*b
print(X)

fig, (ax) = plt.subplots(1, 
	figsize=(8, 8), 
	subplot_kw=Dict("projection" => "3d"))

ax.scatter3D(x, y ,z)
xp=range(0, 1, length=10)
yp=range(0, 1, length=10)

xpm,ypm=meshgrid(xp,yp)
ax.plot_wireframe( xpm, ypm, X[1]*xpm.^2 + X[2]*cos.(6*xpm) + X[3]*ypm.^2 + X[4]*exp.(2*ypm) .+ X[5])
ax.view_init(elev=20, azim=230)

plt.show()

処理結果

フクさん
フクさん

処理結果は以下。

正規方程式で多変量多項式回帰分析(関数項あり)(Julia)、Figure 1
[3.9599058103789586; -5.010636186883172; 2.137338779263279; 1.1442906477838677; 1.8263511178012835;;]

考察

太郎くん
太郎くん

これもちょっと誤差が載ってる感じか。
サンプル点数を10000点にするとどんな感じになる?

フクさん
フクさん

こんな感じになるな。

[3.995277396089315; -4.99837832634382; 3.0528113303488764; 0.9930704870656972; 2.012192078222188;;]
太郎くん
太郎くん

だいぶ良いねー。
やっぱりサンプル点数次第なところはあるのか。

まとめ

フクさん
フクさん

まとめだよ。

  • 正規方程式による多変量多項式回帰分析をJuliaで実施。
  • 誤差はあるものの目的の係数の算出はできている。
    • 誤差の出方はサンプル点数次第。

バックナンバーはこちら。

コメント

タイトルとURLをコピーしました