MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/
はじめに
の、
MATLAB,Python,Scilab,Julia比較 第2章 その76【多変量多項式回帰分析(関数項)⑤】
を書き直したもの。
正規方程式を用いた、多変量多項式回帰分析(関数項あり)について。
今回は、Juliaで演算してみる。
正規方程式、各パラメータ、推定対象の多項式再掲
まずは恒例の正規方程式、多変量多項式回帰分析(関数項あり)で想定するパラメータの再掲。
正規方程式
\(
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()
処理結果
処理結果は以下。
[3.9599058103789586; -5.010636186883172; 2.137338779263279; 1.1442906477838677; 1.8263511178012835;;]
考察
これもちょっと誤差が載ってる感じ。
サンプル点数を10000点にすると以下のようになる。
[3.995277396089315; -4.99837832634382; 3.0528113303488764; 0.9930704870656972; 2.012192078222188;;]だいぶ良い。
やはりサンプル点数次第なところはある。
まとめ
- 正規方程式による多変量多項式回帰分析をJuliaで実施。
- 誤差はあるものの目的の係数の算出はできている。
- 誤差の出方はサンプル点数次第。
MATLAB、Python、Scilab、Julia比較ページはこちら
コメント