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

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

以下の例では、与えられた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 = @(x, b) b(1) + b(2) * x .* exp(b(3) * x);
>> b0 = [1 1 1]';
>>
>> % 重みを考慮しない計算(図の青実線)
>> [yest1, b, irlt, iter, corm, covm] = leasqr(x, y, b0, f, 1.0e-8, 100);
>> disp(b)     # 推定した回帰係数
  -3.8908
   3.3949
   0.1085
>> disp(irlt)  # 収束判定(収束したら1)
1
>> disp(covm)  # 推定した回帰係数の分散共分散行列
   6.0731e+00  -1.5657e+00   4.4934e-02
  -1.5657e+00   5.5200e-01  -1.7014e-02
   4.4934e-02  -1.7014e-02   5.3738e-04
>>
>> % 重みを考慮した計算(図の赤実線)
>> [yest2, b, irlt, iter, corm, covm] = leasqr(x, y, b0, f, 1.0e-8, 100, w);
>> disp(b)     # 推定した回帰係数
   0.2421
   1.9985
   0.1623
>> disp(irlt)  # 収束判定(収束したら1)
1
>> disp(covm)  # 推定した回帰係数の分散共分散行列
   6.8549e-01  -2.3642e-01   1.6110e-02
  -2.3642e-01   8.1785e-02  -5.5772e-03
   1.6110e-02  -5.5772e-03   3.8044e-04
>>
>> % 図の描画
>> plot(x, yest1, 'b', x, yest2, 'r', x, y, 'ko', 'MarkerFaceColor', 'k')