非線形データ適合

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

問題の設定

問題は、関数を次のデータに当てはめることです。

xdata = ...
  [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 = xdata(:,1);
y = xdata(:,2);
plot(t,y,"o")
title("Data points")

次の関数

$$y(t)=c_1\exp (-{\lambda }_{1}t)+c_2\exp (-{\lambda }_{2}t)$$

をデータに当てはめてみます。

lsqcurvefit を使用した解法

lsqcurvefit を使用して関数をデータに当てはめるには、次のように、$$y(t)$$ 関数のパラメーターを変数 x にマッピングします。

x(1) = $$c_1$$

x(2) = $${\lambda }_1$$

x(3) = $$c_2$$

x(4) = $${\lambda }_2$$

関数ハンドルで、当てはめるパラメーター x とともにデータ xdata を定義します。

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

$$c_1$$ = 1、$${\lambda }_1$$ = 1、$$c_2$$ = 1、$${\lambda }_2$$ = 0 のように、初期点 x0 を任意に設定します。

x0 = [1 1 1 0];

lsqcurvefit ソルバーを実行し、結果の当てはめをプロットします。

[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.

<stopping criteria details>

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.↵↵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.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
       bestfeasible: []
    constrviolation: []

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

fminunc を使用した解法

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

Fsumsquares = @(x)sum((F(x,t) - y).^2);
[xunc,ressquared,eflag,outputu] = ...
    fminunc(Fsumsquares,x0)

Local minimum found.

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

<stopping criteria details>

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.↵↵Optimization completed because the size of the gradient is less than↵the value of the optimality tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The first-order optimality measure, 9.409720e-07, is less ↵than options.OptimalityTolerance = 1.000000e-06.'

fminunc を使用して lsqcurvefit と同じ解が得られますが、より多くの関数評価を行っていることに注意してください。fminunc のパラメーターは、lsqcurvefit のパラメーターとは逆の順序になります。より大きな $$\lambda$$ は $${\lambda }_1$$ ではなく $${\lambda }_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$$ について線形であることに注意してください。これは、$${\lambda }_1$$ と $${\lambda }_2$$ の任意の値に対して、バックスラッシュ演算子を使用して、最小二乗問題を解く $$c_1$$ と $$c_2$$ の値を見つけることができるということを意味します。

2 次元問題として問題を解き直し、$${\lambda }_1$$ および $${\lambda }_2$$ の最適な値を探します。$$c_1$$ および $$c_2$$ の値は、上記で説明したように、バックラッシュ演算子を使用して各ステップで計算されます。

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 the 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.
end

lsqcurvefit を使用して問題を解きます。2 次元の初期点 [$${\lambda }_1{\lambda }_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.

<stopping criteria details>

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.↵↵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.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
       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 パラメーター問題に対して、同じ不適切な $${\lambda }_1$$ および $${\lambda }_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.

<stopping criteria details>

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.↵↵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.↵↵<stopping criteria details>↵↵Optimization stopped because the relative sum of squares (r) is changing↵by less than options.FunctionTolerance = 1.000000e-06.'
       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,\r' ...
   '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.

参考

トピック

• さまざまな問題ベースのアプローチを使用した非線形データ適合
• 非線形最小二乗法とヤコビアンの有無
• lsqcurvefit を使用した非線形曲線近似