MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/
はじめに
の、
MATLAB,Python,Scilab,Julia比較 第2章 その75【多変量多項式回帰分析(関数項)④】
を書き直したもの。
正規方程式を用いた、多変量多項式回帰分析(関数項あり)について。
今回は、Scilabで演算してみる。
正規方程式、各パラメータ、推定対象の多項式再掲
まずは恒例の正規方程式、多変量多項式回帰分析(関数項あり)で想定するパラメータの再掲。
正規方程式
\(
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
\)
これをScilabで実現する。
Scilab
Scilabコードは以下になる。
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;
disp(X);
scatter3d(x, y ,z);
xp=linspace(0, 1, 10);
yp=linspace(0, 1, 10);
[xpm,ypm]=meshgrid(xp,yp);
mesh( xp, yp, X(1)*xpm.^2 + X(2)*cos(6*xpm) + X(3)*ypm.^2 + X(4)*exp(2*ypm)+X(5));処理結果
処理結果は以下。
3.9633006
-4.9449502
5.5469140
0.5920788
2.4355476考察
これもちょっと誤差大き目。
恐らく、サンプル点数増やすと改善されるパターン。
10000点で再度実施した結果が以下になる。
3.9963030
-5.0002834
2.8950522
1.0155636
1.9884052かなり期待通りの数値が出てきた。
あと、毎回言ってるかもしれないけど、散布図都合でScilabのVersion6を使用している。
あと、コードもMATLABコードをコピペレベルで一緒。
まとめ
- 正規方程式による多変量多項式回帰分析をScilabで実施。
- 誤差はあるものの目的の係数の算出はできている。
- サンプル点数を増やせば誤差は減る。
- コード自体はMATLABコードのコピペ。
- scatter3をscatter3dに書き換えた程度。
MATLAB、Python、Scilab、Julia比較ページはこちら
コメント