重みによるモデル近似の改善

この例では、比較のために線形最小二乗法と重み付き最小二乗法の両方を使用して多項式モデルをデータに当てはめる方法を説明します。

関数randnを使用して、さまざまな正規分布からサンプル データを生成します。

rng("default")  % For reproducibility

rnorm = [];
idx = [];

for k=1:20
    r = k*randn([20,1]) + (1/20)*(k^3);
    rnorm = [rnorm;r];
    idx = [idx;ones(20,1).*k];
end

従属変数 rnorm には、20 個の正規分布からのサンプル データが含まれます。独立変数 idx には、rnorm 内の 2 つの要素が同じ正規分布からサンプリングされているかどうかを示す整数が含まれます。

3 次多項式モデルを idx と rnorm に当てはめます。モデルを当てはめるときに使用された係数推定値とアルゴリズムに関する情報を返します。

[p3fit,~,fitinfo] = fit(idx,rnorm,"poly3")

p3fit = 
     Linear model Poly3:
     p3fit(x) = p1*x^3 + p2*x^2 + p3*x + p4
     Coefficients (with 95% confidence bounds):
       p1 =     0.05156  (0.0438, 0.05932)
       p2 =    -0.03993  (-0.2875, 0.2076)
       p3 =      0.1418  (-2.124, 2.408)
       p4 =      0.0462  (-5.585, 5.678)

fitinfo = struct with fields:
        numobs: 400
      numparam: 4
     residuals: [400×1 double]
      Jacobian: [400×4 double]
      exitflag: 1
     algorithm: 'QR factorization and solve'
    iterations: 1

p3fit には、モデル係数 (95% 信頼区間) の推定値が含まれます。多項式モデルを当てはめる場合の既定の近似法は線形最小二乗法です。fitinfo には、係数をデータに当てはめるときに使用された近似アルゴリズムについての情報が含まれます。データ内の誤差は、fitinfo に保存される残差によって推定できます。

ステム プロットを使用して残差をプロットします。

stem(idx,fitinfo.residuals)
xlabel("idx")
ylabel("residuals")

残差のプロットは、idx 内の値が大きくなるにつれて、分散が大きくなっていることを示しています。idx のさまざまな値に対して分散が一定でないことは、重み付き最小二乗近似法の方がモデル係数の計算に適していることを示しています。

後で重みを保存するために、ゼロのベクトルを作成します。

W = zeros(length(rnorm),1);

ここで提供する重みは、残差分散を変換し、それが idx のさまざまな値に対して一定になるようにします。rnorm 内の要素ごとに、idx 内の対応する値の残差分散の逆数として重みを定義します。次に、重み付きのモデルを当てはめます。

for k=1:20
    rnorm_idx = rnorm(idx==k);
    recipvar = 1/var(rnorm_idx);
    w = (idx==k).*recipvar;
    W = W+w;
end
[wp3fit,~,wfitinfo] = fit(idx,rnorm,"poly3","Weights",W)

wp3fit = 
     Linear model Poly3:
     wp3fit(x) = p1*x^3 + p2*x^2 + p3*x + p4
     Coefficients (with 95% confidence bounds):
       p1 =     0.04894  (0.04419, 0.0537)
       p2 =     0.03601  (-0.08744, 0.1595)
       p3 =     -0.4262  (-1.253, 0.4009)
       p4 =      0.9836  (-0.1959, 2.163)

wfitinfo = struct with fields:
        numobs: 400
      numparam: 4
     residuals: [400×1 double]
      Jacobian: [400×4 double]
      exitflag: 1
     algorithm: 'QR factorization and solve'
    iterations: 1

wp3fit と wfitinfo には、重み付き最小二乗近似の結果が含まれます。

p3fit、wp3fit、rnorm を同じプロットに表示します。

plot(p3fit,idx,rnorm)
hold on
plot(wp3fit,"g")
xlabel("idx")
ylabel("rnorm")
legend(["rnorm","linear least-squares fit","weighted least-squares fit"])
hold off

プロットは、wp3fit が p3fit を密接に追従していることを示しています。

wp3fit の方が p3fit より適切な近似かどうかは、残差をプロットすることによって判断できます。

stem(idx,wfitinfo.residuals)
xlabel("idx")
ylabel("residuals")

この出力は、wp3fit 残差の方が、p3fit 残差より小さいことを示しています。また、wp3fit 残差の分散の方が、p3fit 残差の分散より、idx の値が変わっても類似しています。

参考

関数

• fit | stem

トピック

• 最小二乗近似の紹介