［Octave］重み付き非線形最小二乗法を行う

optimパッケージのnlinfit関数を使う。

以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * x * exp(b3 * x)（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

>> % 計算の準備
>> pkg load optim
>> x = [1 2 3 5 6 7 8 9]';
>> y = [1 4 10 25 32 49 64 75]';
>> w = repmat(1, 1, length(x))';
>> w(3) = 100;
>> w(5) = 100;
>> f = @(b, x) (b(1) + b(2) * x .* exp(b(3) * x));
>> b0 = [1 1 1];
>>
>> % 重みを考慮しない計算(図の青実線)
>> [b, res, jj, covm, rse] = nlinfit(x, y, f, b0);
>> yest1 = f(b, x);
>> disp(b)     % 推定した回帰係数
  -3.8908
   3.3950
   0.1085
>> disp(rse)   % 回帰値の標準誤差
6.2601
>> disp(covm)  % 推定した回帰係数の分散共分散行列
   6.0731e+00  -1.5657e+00   4.4934e-02
  -1.5657e+00   5.5201e-01  -1.7014e-02
   4.4934e-02  -1.7014e-02   5.3738e-04
>>
>> % 重みを考慮した計算(図の赤実線)
>> [b, res, jj, covm, rse] = nlinfit(x, y, f, b0, [], 'weights', w);
>> yest2 = f(b, x);
>> disp(b)     % 推定した回帰係数
  -0.030452
   2.060791
   0.159292
>> disp(rse)   % 回帰値の標準誤差
13.991
>> disp(covm)  % 推定した回帰係数の分散共分散行列
   1.3955e+00  -3.7178e-01   2.1961e-02
  -3.7178e-01   1.1270e-01  -6.9215e-03
   2.1961e-02  -6.9215e-03   4.3259e-04
>>
>> % 図の描画
>> plot(x, yest1, 'b', x, yest2, 'r', x, y, 'ko', 'MarkerFaceColor', 'k')