第3回 MATLAB seminar 2021年06月22日

MATLAB/Scilabで理解する数値計算 の50ページくらいまでの基礎文法をもう一度復習した. とくに, ブーリアンが便利であることを強調した. 例えば, 条件にあう要素をとってきたり, 数えたりできる.

簡単な構文で, for, if文を練習した. 標準正規分布からサンプルを $n$ 個取り, 0.5を超えた回数を数えることで, 標準正規分布に従う確率変数が0.5を超える確率を計算した. その値を分布表と比較して, $n$ が増えると精度が高くなることを確認した.

nTry = 1000; % サンプルの個数
count = 0;

for I = 1:nTry
    n = randn(1); % 乱数を発生させる. 
    if n > 0.5
        count = count + 1;
    end
end
count/nTry % 0.5をこえた割合を計算

前回, 最小自乗法の基礎となる正規方程式を紹介し, 1次, 2次多項式によるフィッティングを行った. $$ \large{ \boldsymbol{p} = (\boldsymbol{X}^{\mathrm{T}}\boldsymbol{X})^{-1}\boldsymbol{X}^{\mathrm{T}}\boldsymbol{y}, } $$ $$ \large{ \boldsymbol{X} = \begin{pmatrix} 1 & x_{1}^{\,} & x_{1}^{2} & \ldots & x_{1}^{d} \\ 1 & x_{2}^{\,} & x_{2}^{2} & \ldots & x_{2}^{d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n}^{\,} & x_{n}^{2} & \ldots & x_{n}^{d} \end{pmatrix}. } $$ ここで,

• $\boldsymbol{x}$：説明変数 ( データ )
• $\boldsymbol{y}$：目的変数 ( データ )
• $n$：データの次元 ( 個数, 今回は50個 )
• $d$：フィッティングしたい多項式の次数, $d$ 次多項式
• $p$：$d$ 次多項式の係数. 切片含めると, $p+1$ ベクトルとなる.

今回は, これを函数にした. つまり, データセット $\boldsymbol{x}$, $\boldsymbol{y}$, 多項式の次数を渡せば, 係数 $p$ を計算してくれる函数を作った. さらに, フィッティングの結果を図示させる機能を加えた ( doPlotには, trueかfalseを入れる ) .

うまく行列を使うことで, コードが短くなることを確認した. 特に, 10行目など.

function [p] = myMle(x, y, deg, doPlot)
    
    X = x .^ (0:deg);
    p = inv(X'*X)*X'*y; % invは使わないほうがいい
    
    if doPlot
        minX = min(x);
        maxX = max(x);
        x_fit = linspace(minX-abs(minX/10), maxX+abs(maxX/10), 100)';
        y_fit = (x_fit .^ (0:deg)) * p;
        
        scatter(x, y); hold on
        plot(x_fit, y_fit); hold off    
    end
    
end

データ