GNU octave で非線形フィッティング

Februray 3, 2010

Gnuplot の fit 関数と同じアルゴリズム (Levenberg-Marquardt 法、繰り返し計算で最小二乗フィットを求める) を GNU octave でやる方法。

用意するもの

• Mac OS X 10.6.2
• GNU octave v. 3.2.3
• leasqr.m と dfdp.m (ここから)

準備

• octave にパスを通す。
• フィットしたいデータを格納したファイル、leasqr.m、dfdp.m を同じディレクトリに置く。
• データファイル中では、一行に空白文字区切りで一組のデータが入っているとする。

フィッティング実行

三角関数をデータに当てはめる例。

• 以下のように octave で実行する。
  function obj = phi(x,p) % フィットしたい関数の定義
    A = p(1);
    B = p(2);
    C = p(3);
    D = p(4);
    obj = (A * cos(B * x + C) + D);
  endfunction
  
  load xdata.txt; % x座標値/独立変数値、上の関数での x の値、ベクトル
  X = xdata;
  load datafile.txt; % データファイル (y座標値/従属変数値、obj の値、各行が X と同じサイズのベクトル)
  data = datafile; % ファイル中のデータを行列に入れる
  
  for i = length(data(:,1)) : 1 : 1
    A = 1; % パラメータの初期値、何か適切な値にする
    B = 1; % ただ全部 0 だと leasqr() が受け付けてくれない
    C = 1;
    D = 1;
    [f, x, kvg] = leasqr(X, data, [A B C D], "phi", []); % フィッティング実行
    Par(i,1:4) = x(1:4); % 求められた最適パラメータ
    Fit(i,1:N) = f; % 最適パラメータでの関数値
    if (kvg == 1) Suc(i) = kvg; % 最適化がうまくいったかどうか
    else Suc(i) = 0;
    endif
  endfor;
  

注意

• leasqr が返す kvg は 1 か 0 の値を取るスカラー値ということになっているが (1 は LM 法が収束したこと、0 は収束しなかったことを示す)、Zero matrix を返すことがあり、kvg をスカラー変数に代入しようとすると、実行時に 「型が違うよ」エラーになることがある。だから上のスクリプトでは if 文にしている。
• 上のスクリプトでは、leasqr が返す値の一部しか使っていない。