Нелинейная подгонка данных

Этот пример показов аппроксимацию нелинейной функции к данным с помощью несколького Optimization Toolbox™ алгоритмов.

Setup задачи

Примите во внимание следующие данные:

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')

Figure contains an axes. The axes with title Data points contains an object of type line.

% hold off

Мы хотели бы подогнать функцию

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

к данным.

Подход к решению с использованием lsqcurvefit

The 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.8852e-06
       iterations: 6
        funcCount: 35
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0096
          message: '...'

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

Figure contains an axes. The axes with title Data points contains 2 objects of type line.

Подход к решению с использованием 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: '...'

Заметьте, что fminunc найдено то же решение, что и lsqcurvefitно для этого потребовалось еще много вычисления функции. Параметры для fminunc находятся в противоположном порядке, как для lsqcurvefit; большим ламом является лама (2), а не лама (1). Это неудивительно, порядок переменных произволен.

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), которые решают задачу наименьших квадратов.

Теперь мы переделываем задачу как двумерную задачу, ища лучшие значения 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, начиная с двумерной начальной точки 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.4018e-06
       iterations: 10
        funcCount: 33
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0080
          message: '...'

Эффективность двумерного решения аналогична эффективности четырехмерного решения:

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.

Задача разделения более устойчива к первоначальным предположениям

Выбор плохой начальной точки для исходной четырехпараметрической задачи приводит к локальному решению, которое не является глобальным. Выбор начальной точки с одинаковыми значениями bad 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.0057
       iterations: 32
        funcCount: 165
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0021
          message: '...'

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

Figure contains an axes. The axes with title Data points contains 3 objects of type line. These objects represent Data, Global fit, Bad local fit.

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.

Похожие темы