ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

非線形データ近似

この例では、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)

hold on
plot(t,F(x,t))
hold off
Local minimum possible.

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




x =

    3.0068   10.5869    2.8891    1.4003


resnorm =

    0.1477


exitflag =

     3


output = 

    firstorderopt: 7.8802e-06
       iterations: 6
        funcCount: 35
     cgiterations: 0
        algorithm: 'trust-region-reflective'
          message: 'Local minimum possible.

lsqcurvefit stopped because t...'

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 default value of the function tolerance.




xunc =

    2.8890    1.4003    3.0069   10.5862


ressquared =

    0.1477


eflag =

     1


outputu = 

       iterations: 30
        funcCount: 185
         stepsize: 1
    firstorderopt: 2.9476e-05
        algorithm: 'medium-scale: Quasi-Newton line search'
          message: 'Local minimum found.

Optimization completed because t...'

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)
fprintf(['There were %d function evaluations using fminunc,' ...
    ' and %d using lsqcurvefit.'], ...
    outputu.funcCount,output.funcCount)
There were 30 iterations using fminunc, and 6 using lsqcurvefit.
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 default value of the function tolerance.




x2 =

   10.5861    1.4003


resnorm2 =

    0.1477


exitflag2 =

     3


output2 = 

    firstorderopt: 4.4085e-06
       iterations: 10
        funcCount: 33
     cgiterations: 0
        algorithm: 'trust-region-reflective'
          message: 'Local minimum possible.

lsqcurvefit stopped because t...'

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)

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)
Local minimum possible.

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




xbad =

  -22.9036    2.4792   28.0273    2.4791


resnormbad =

    2.2173


exitflagbad =

     3


outputbad = 

    firstorderopt: 0.0058
       iterations: 32
        funcCount: 165
     cgiterations: 0
        algorithm: 'trust-region-reflective'
          message: 'Local minimum possible.

lsqcurvefit stopped because t...'

The residual norm at the good ending point is 0.147723, and the residual norm at the bad ending point is 2.217300.

この情報は役に立ちましたか?