MATLABのlsqcurvefitを用いた最小二乗法で、ヤコビアンを指定する方法

フィッティングしたい関数myfuncがスカラー変数xとNp個のパラメータを受け取り、スカラーyを返すような関数とする。この時、数学の定義としてはヤコビアンは1行Np列である。しかし、matlabのプログラミングとしてはxに配列を投げてまとめて処理することが多々ある。このときヤコビアンはどうしたら良いのか、matlabのヘルプドキュメントを読んでもよく分からない。結論としては、xにNx個の要素をもつ配列を入力する場合は、ヤコビアンはNx行Np列にすれば良いらしいが、根拠となる正式ドキュメントは見つからなかった（myfuncも戻り値がNy個の要素をもつ配列の場合にどうすれば良いのかは、調べてないので不明）。

以下に例を示す。次のスクリプトをtest_lsqcurvefit.mとして保存して実行すると、下図の結果を出力する。

function test_lsqcurvefit
p0 = [1.5,  1.5]; %define these yourself

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ここでダミーのデータを用意
xdata = linspace(0,1);
[ydata,~] = myfunc(p0, xdata);
a=-0.1; b = 0.1;
err = a + (b-a)*rand(1,length(xdata));
ydata = ydata + err;

err = a + (b-a)*rand(1,length(p0));
p0 = p0 + err;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% ここでフィッティングの条件を設定
options = optimset('TolFun' , 1e-12,'TolX' ,1e-12,...
    'Jacobian', 'on', 'LargeScale','on');
lb = [0 0];
ub = [inf inf];

%%% ここでフィッティング
parEst = lsqcurvefit(@myfunc,p0,xdata,ydata,lb,ub,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 結果を出力
figure
hold all
plot(xdata, ydata, 'b+')

time = 0:0.01:xdata(end);
[FitCurve, ~] = myfunc(p0,time);
plot(time, FitCurve, 'g-')

[FitCurve, ~] = myfunc(parEst,time);
plot(time, FitCurve, 'r-')

legend('data', 'initial curve', 'fitting result');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% フィッティングしたい関数（モデル）
function [Sig Jacob] = myfunc(p,xdata)
% This function is called by lsqcurvefit.
% x is a vector which contains the coefficients of the
% equation. xdata is the data set
% passed to lscurvefit.

A = p(1);
B = p(2);

Sig = A.*exp(-xdata./B);
Jp1 = exp(-xdata./B);
Jp2 = A.*(xdata./(B.^2)).*exp(-xdata./B);
Jacob = [Jp1; Jp2]';
end