Main Content

非線形データ適合

ライブ スクリプトを開く

この例では、Optimization Toolbox™ のいくつかのアルゴリズムを使用して非線形関数をデータに当てはめる方法を説明します。

問題の設定

以下のデータを考えます。

Data = ...
  [0.0000    5.8955
   0.1000    3.5639
   0.2000    2.5173
   0.3000    1.9790
   0.4000    1.8990
   0.5000    1.3938
   0.6000    1.1359
   0.7000    1.0096
   0.8000    1.0343
   0.9000    0.8435
   1.0000    0.6856
   1.1000    0.6100
   1.2000    0.5392
   1.3000    0.3946
   1.4000    0.3903
   1.5000    0.5474
   1.6000    0.3459
   1.7000    0.1370
   1.8000    0.2211
   1.9000    0.1704
   2.0000    0.2636];

これらのデータ点をプロットしてみましょう。

t = Data(:,1);
y = Data(:,2);
% axis([0 2 -0.5 6])
% hold on
plot(t,y,'ro')
title('Data points')

% hold off

関数

y = c(1)*exp(-lam(1)*t) + c(2)*exp(-lam(2)*t)

をデータに近似したいと思います。

lsqcurvefit を使用した解法

関数 lsqcurvefit では、このような問題を簡単に解くことができます。

まず、パラメーターを変数 x の項で定義します。

x(1) = c(1)

x(2) = lam(1)

x(3) = c(2)

x(4) = lam(2)

次に、パラメーター x およびデータ t の関数として曲線を定義します。

F = @(x,xdata)x(1)*exp(-x(2)*xdata) + x(3)*exp(-x(4)*xdata);

次のように、初期点 x0 を任意に設定します。c(1) = 1, lam(1) = 1, c(2) = 1, lam(2) = 0:

x0 = [1 1 1 0];

ソルバーを実行し、結果の近似をプロットします。

[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)
Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

x = 1×4

    3.0068   10.5869    2.8891    1.4003


resnorm = 0.1477

exitflag = 3

output = struct with fields:
      firstorderopt: 7.8830e-06
         iterations: 6
          funcCount: 35
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 0.0096
            message: 'Local minimum possible....'
       bestfeasible: []
    constrviolation: []


hold on
plot(t,F(x,t))
hold off

fminunc を使用した解法

fminunc を使用して問題を解くには、残差の二乗和として目的関数を設定します。

Fsumsquares = @(x)sum((F(x,t) - y).^2);
opts = optimoptions('fminunc','Algorithm','quasi-newton');
[xunc,ressquared,eflag,outputu] = ...
    fminunc(Fsumsquares,x0,opts)
Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

xunc = 1×4

    2.8890    1.4003    3.0069   10.5862


ressquared = 0.1477

eflag = 1

outputu = struct with fields:
       iterations: 30
        funcCount: 185
         stepsize: 0.0017
     lssteplength: 1
    firstorderopt: 2.9662e-05
        algorithm: 'quasi-newton'
          message: 'Local minimum found....'

fminunc を使用して lsqcurvefit と同じ解が得られましたが、より多くの関数評価を行ったことに注意してください。fminunc のパラメーターは、lsqcurvefit のパラメーターとは逆の順序になります。より大きな lam は lam(1) ではなく lam(2) です。変数の順序は任意なので、これは当然です。

fprintf(['There were %d iterations using fminunc,' ...
    ' and %d using lsqcurvefit.\n'], ...
    outputu.iterations,output.iterations)
There were 30 iterations using fminunc, and 6 using lsqcurvefit.

fprintf(['There were %d function evaluations using fminunc,' ...
    ' and %d using lsqcurvefit.'], ...
    outputu.funcCount,output.funcCount)
There were 185 function evaluations using fminunc, and 35 using lsqcurvefit.

線形問題と非線形問題の分割

近似問題は、パラメーター c(1) および c(2) の線形であることに注意してください。これは、lam(1) および lam(2) の任意の値に対し、バックスラッシュ演算子を使用して、最小二乗問題を解く c(1) および c(2) の値を見つけることができるということです。

2 次元問題として問題を解き直し、lam(1) および lam(2) の最適な値を探します。c(1) および c(2) の値は、上記で説明したように、バックラッシュ演算子を使用して各ステップで計算されます。

type fitvector
function yEst = fitvector(lam,xdata,ydata)
%FITVECTOR  Used by DATDEMO to return value of fitting function.
%   yEst = FITVECTOR(lam,xdata) returns the value of the fitting function, y
%   (defined below), at the data points xdata with parameters set to lam.
%   yEst is returned as a N-by-1 column vector, where N is the number of
%   data points.
%
%   FITVECTOR assumes the fitting function, y, takes the form
%
%     y =  c(1)*exp(-lam(1)*t) + ... + c(n)*exp(-lam(n)*t)
%
%   with n linear parameters c, and n nonlinear parameters lam.
%
%   To solve for the linear parameters c, we build a matrix A
%   where the j-th column of A is exp(-lam(j)*xdata) (xdata is a vector).
%   Then we solve A*c = ydata for the linear least-squares solution c,
%   where ydata is the observed values of y.

A = zeros(length(xdata),length(lam));  % build A matrix
for j = 1:length(lam)
   A(:,j) = exp(-lam(j)*xdata);
end
c = A\ydata; % solve A*c = y for linear parameters c
yEst = A*c; % return the estimated response based on c

lsqcurvefit を使用して問題を解きます。2 次元の初期点 lam(1)、lam(2) から始めます。

x02 = [1 0];
F2 = @(x,t) fitvector(x,t,y);

[x2,resnorm2,~,exitflag2,output2] = lsqcurvefit(F2,x02,t,y)
Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

x2 = 1×2

   10.5861    1.4003


resnorm2 = 0.1477

exitflag2 = 3

output2 = struct with fields:
      firstorderopt: 4.4087e-06
         iterations: 10
          funcCount: 33
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 0.0080
            message: 'Local minimum possible....'
       bestfeasible: []
    constrviolation: []

2 次元解法の効率は、4 次元解法の効率と同等です。

fprintf(['There were %d function evaluations using the 2-d ' ...
    'formulation, and %d using the 4-d formulation.'], ...
    output2.funcCount,output.funcCount)
There were 33 function evaluations using the 2-d formulation, and 35 using the 4-d formulation.

初期推定に対してよりロバストな分割問題

元の 4 パラメーター問題に対して不適切な始点を選択すると、大域的でない局所的な解がもたらされます。分割 2 パラメーター問題に対して、同じ不適切な lam(1) および lam(2) 値をもつ始点を選択すると、大域的な解がもたらされます。これを示すために、相対的に不適切である局所的な解をもたらす始点をもつ元の問題を再実行し、大域的な解によって得られた近似と比較します。

x0bad = [5 1 1 0];
[xbad,resnormbad,~,exitflagbad,outputbad] = ...
    lsqcurvefit(F,x0bad,t,y)
Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

xbad = 1×4

  -22.9036    2.4792   28.0273    2.4791


resnormbad = 2.2173

exitflagbad = 3

outputbad = struct with fields:
      firstorderopt: 0.0056
         iterations: 32
          funcCount: 165
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 0.0021
            message: 'Local minimum possible....'
       bestfeasible: []
    constrviolation: []


hold on
plot(t,F(xbad,t),'g')
legend('Data','Global fit','Bad local fit','Location','NE')
hold off

fprintf(['The residual norm at the good ending point is %f,' ...
   ' and the residual norm at the bad ending point is %f.'], ...
   resnorm,resnormbad)
The residual norm at the good ending point is 0.147723, and the residual norm at the bad ending point is 2.217300.

関連するトピック