バックナンバーはこちら。
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 & y_1 & 1\\
x_2 & y_2 & 1\\
\vdots & \vdots & \vdots\\
x_n & y_n & 1\\
\end{bmatrix},
\vec{x}=
\begin{bmatrix}
\alpha\\
\beta\\
\gamma
\end{bmatrix},
\vec{b}=
\begin{bmatrix}
z_1\\
z_2\\
\vdots\\
z_n
\end{bmatrix}
\)
推定対象の多項式
\(
z=3x-2y+5
\)
今回は、これを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(1, n)
y = rand(1, n)
z = 3*x.-2*y.+5+rand(1, n)*2 .-1
A=[x' 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=5)
yp=range(0, 1, length=5)
xpm,ypm=meshgrid(xp,yp)
ax.plot_wireframe( xpm, ypm, X[1].*xpm.+X[2].*ypm.+X[3])
ax.view_init(elev=20, azim=230)
plt.show()処理結果
そして処理結果。
[2.9832030959622173; -1.9194487071821298; 4.9303170774938145;;]考察
演算結果としてはOKそうだね。
コードとしては・・・。
そっか、meshgridが無いから自作関数を置いておく必要があるのか。
3Dグラフの方は、PyPlotがmatplotlibのラッパーであることから、MATLABよりPython的な仕様になるね。
でもまぁ、基本的な流れは一緒だね。
まとめ
まとめだよ。
- 正規方程式による重回帰分析をJuliaで実施。
- 誤差はあるものの目的の係数の算出はできている。
- 3DグラフはPython寄りの仕様。
- PyPlotがmatplotlibのラッパーであるため。
バックナンバーはこちら。
コメント