Нелинейный подбор кривой данных

В этом примере показано, как соответствовать нелинейной функции к данным с помощью нескольких алгоритмов 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')

% 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, убегите (1) = 1, c (2) = 1, убегите (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.8876e-06
       iterations: 6
        funcCount: 35
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 0.0096
          message: '...'

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: '...'

Заметьте тот 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). Это означает для любых значений бегства (1) и бегства (2), мы можем использовать оператор обратной косой черты, чтобы найти значения c (1) и c (2), которые решают задачу наименьших квадратов.

Мы теперь переделываем проблему как двумерную проблему, ища оптимальные значения бегства (1) и убегаем (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, запуск с двумерного начального бегства точки (1), убегите (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.4071e-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.

Разделите проблема Более устойчива к Исходному предположению

Выбор плохой начальной точки для исходной проблемы с четырьмя параметрами приводит к локальному решению, которое не является глобальной переменной. Выбор начальной точки с тем же плохим бегством (1) и бегством (2) значения для проблемы 2D параметра разделения приводит к глобальному решению. Чтобы показать это, мы повторно выполняем исходную проблему со стартовой точкой, которая приводит к относительно плохому локальному решению, и сравните получившуюся подгонку с глобальным решением.

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

  -24.6714    2.4788   29.7951    2.4787

resnormbad = 2.2173

exitflagbad = 3

outputbad = struct with fields:
    firstorderopt: 5.1531e-05
       iterations: 38
        funcCount: 195
     cgiterations: 0
        algorithm: 'trust-region-reflective'
         stepsize: 2.6834e-04
          message: '...'

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.

Документация