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