Matlab — подбор данных методом наименьших квадратов — функция стоимости с дополнительным ограничением

В настоящее время я работаю над некоторым кодом MatLab, чтобы подогнать экспериментальные данные к сумме экспонент, следуя методу, описанному в данный документ.

Согласно документу, данные должны следовать следующему уравнению (написанному псевдокодом):

y = sum(v(i)*exp(-x/tau(i)),i=1..n)

Здесь tau(i) — это набор n предопределенных констант. Количество констант также определяет размер суммирования и, следовательно, размер v. Например, мы можем попытаться подобрать сумму из 100 экспонент, каждая из которых имеет разные tau(i) для наших данных. Однако из-за характера подгонки и экспоненциальной суммы нам необходимо добавить еще одно ограничение к задаче и, следовательно, к функции стоимости используемого метода наименьших квадратов.

Обычно функция стоимости метода наименьших квадратов определяется следующим образом:

(y_data - sum(v(i)*exp(-x/tau(i)),i=1..n)^2

И это должно быть сведено к минимуму. Однако, чтобы предотвратить переобучение, которое сделало бы постоянный во времени спектр чрезвычайно зашумленным, в документе добавлено следующее ограничение на функцию стоимости:

|v(i) - v(i+1)|^2

Из-за этого дополнительного ограничения, насколько я знаю, обычные алгоритмы, такие как lsqcurvefit, больше нельзя использовать, и мне приходится использовать fminsearch для поиска минимума моей функции стоимости наименьших квадратов с ограничением. Функция, которая должна быть минимизирована, по моему мнению, следующая:

(y_data - sum(v(i)*exp(-x/tau(i)),i=1..n)^2 + sum(|v(j) - v(j+1)|^2,j=1..n-1)

Моя попытка закодировать это в MatLab заключается в следующем. Сначала мы определяем функцию в функциональном скрипте, затем используем fminsearch, чтобы минимизировать функцию и получить значения для v.

function res = funcost( v )
%FUNCOST Definition of the function that has to be minimised

%We define a function yvalues with 2 exponentials with known time-constants
% so we know the result that should be given by minimising.
xvalues = linspace(0,50,10000);
yvalues = 3-2*exp(-xvalues/1)-exp(-xvalues/10);

%Definition of 30 equidistant point in the logarithmic scale  
terms = 30;
termsvector = [1:terms];
tau = termsvector;
for i = 1:terms
    tau(i) = 10^(-1+3/terms*i);
end

%Definition of the regular function
res_1 = 3;

for i=1:terms
    res_1 =res_1+ v(i).*exp(-xvalues./tau(i));
end

res_1 = res_1-yvalues;

%Added constraint
k=1;
res_2=0;

for i=1:terms-1
    res_2 = res_2 + (v(i)-v(i+1))^2;
end

res=sum(res_1.*res_1) + k*res_2;

end

fminsearch(@funcost,zeros(30,1),optimset('MaxFunEvals',1000000,'MaxIter',1000000))

Однако этот код дает мне неточные результаты (никакой ошибки, просто неточные результаты), что наводит меня на мысль, что я либо допустил ошибку в кодировании, либо в интерпретации добавленного ограничения для метода наименьших квадратов.


matlab algorithm least-squares data-fitting
person Spears    schedule 09.08.2015    source источник
comment
Откуда вы знаете, что полученные результаты неточны?   -  person A Hernandez    schedule 09.08.2015
comment
Функция, которую я сейчас использую для генерации значений x и y, 3-2*exp(-t)-exp(-t/10), должна давать увеличенные v-амплитуды для tau = 1 и tau = 10. Однако это не тот случай, когда я минимизирую эту функцию и строю график (tau,v).   -  person Spears    schedule 09.08.2015
comment
И это было бы так, если бы вы не добавили дополнительный член к функции стоимости. Также fminsearch является локальным оптимизатором, и решение, вероятно, имеет локальный минимум, отличный от глобального минимума. Закомментируйте лишний термин и дайте ввод, близкий к v(i=10)=-2 и v(i=20)=-1   -  person A Hernandez    schedule 09.08.2015
comment
Спасибо за предложение! Возможно, это так (локальный минимум отличается от глобального минимума)! Если я закомментирую дополнительный термин, я получу решение, которое соответствует данным. Однако тогда амплитуды во временной области чрезвычайно зашумлены (около 2000-4000 в некоторых случаях) и довольно часто меняют знаки. Это также упоминается в документе, что ограничение должно быть включено, чтобы предотвратить это. Поскольку экспоненты не образуют ортогональной основы, я считаю, что не будет уникального решения, поэтому я думаю, что ограничение необходимо.   -  person Spears    schedule 09.08.2015

Ответы (1)


arrow_upward
0
arrow_downward

Я бы попытался ввести дополнительное ограничение следующим образом:

res_2 = max((v(1:(end-1))-v(2:end)).^2);

например вместо минимизации интегральной (суммированной) ошибки она делает minmax.

Вы также можете сделать это ограничение жестким,

if res_2 > some_number
    k = a_very_big_number;
else
    k=0; % or k = a_small_number
end;
person freude    schedule 09.08.2015
comment
Спасибо за предложение! Я попробую это, как только смогу (программа в настоящее время работает для некоторых других значений, поэтому не могу попробовать ее сразу). Я дам вам знать, если это поможет! - person Spears; 09.08.2015
comment
Попробовал ваше ограничение! Это помогает (код быстрее, и я получаю видимый пик при тау = 10, как и ожидалось). Однако пик на 1 едва виден, что также может быть связано с меньшим количеством точек в диапазоне 10^-2..10^0, где эта точка была бы видна. Хотя подгонка не так уж хороша и для более низких постоянных времени (например, не доходит до 0 при t = 0), что может быть связано с характером подгонки. - person Spears; 10.08.2015