バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia2-backnumber/
はじめに
正規方程式を用いた、多変量多項式回帰分析について。
今回は、Scilabで演算してみる。
登場人物
博識フクロウのフクさん
イラスト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 & x_1y_1 & y_1^2 & y_1 &1\\
x_2^2 & x_2y_2 & y_2^2 & y_2 &1\\
\vdots & \vdots & \vdots\\
x_n^2 & x_ny_n & y_n^2 & y_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-5xy+3y^2+y+2
\)
これをScilabで実現する。
Scilabコード
Scilabコードは以下になる。
n = 100;
x = rand(1, n);
y = rand(1, n);
z = 4*x.^2 -5*x.*y + 3*y.^2 + y + 2 + rand(1, n)-0.5;
A=[x'.^2 x'.*y' y'.^2 y' ones(length(x),1)];
b=z';
X=(A'*A)^-1 *A'*b;
disp(X);
scatter3d(x, y ,z);
xp=linspace(0, 1, 5);
yp=linspace(0, 1, 5);
[xpm,ypm]=meshgrid(xp,yp);
mesh( xp, yp, X(1)*xpm.^2+X(2)*xpm.*ypm+X(3)*ypm.^2+X(4)*ypm+X(5));処理結果
処理結果は以下。
3.7901925
-5.1485552
2.9919487
0.8778802
2.1465127考察
これは・・・MATLABのコピペか？！
まぁ、コピペだな。
scatter3をscatter3dに変えたとかはあるけど。
ちなみに、重回帰分析の時も話したけど、Scilabはversion5ではなくversoin6を使用している。
version5だと散布図が作れないんで。
そういうえば、そういうこと言ってたね。
求まった係数も正しいそうだし、ScilabもOKだね。
まとめ
まとめだよ。
- 正規方程式による多変量多項式回帰分析をScilabで実施。
- 誤差はあるものの目的の係数の算出はできている。
- コード自体はMATLABコードのコピペ。
- scatter3をscatter3dに書き換えた程度。
バックナンバーはこちら。
コメント