Одновременная подгонка к данным для оценки общего параметра с данными NaN в Matlab - PullRequest
Новая
       15

Одновременная подгонка к данным для оценки общего параметра с данными NaN в Matlab

1 голос
/ 26 апреля 2019

У меня есть данные (y), измеренные в более чем 4 временных точках (t), для трех концентраций (c).
Таким образом, данные будут такими:

При c1: y11, y12, y13, y14 (измерено в 4 временных точках)

Под c2: y21, y12, y23, y24

при c3: y31, y32, y33, y34

То, что я пытаюсь сделать, - это оценить два параметра c и d, подгоняя одновременно ко всем этим показателям данных при разных концентрациях.
Однако некоторые из этих значений являются NaN.Так, например,

Под c2: y21, y12, NaN, NaN

под c3: y31, y32, y33, NaN

здесьэто код Matlab, который я написал.

%C contains the different concentration values (c1,c2,c3)
[fittedVals,errorVals]=lsqcurvefit(@(xEstimate,thours)model(xEstimate,t,C),initial,t,y,lb,ub);

function output= model(xEstimate,t,C)

intVals=repmat(10^5,3,1);%initial value of the ODE system

[~,values] = ode45(@(t,y)Equations(t,y,C),t,intVals);

function s=Equations(~,y,para)

    a=0.25;
    b=-0.1;
    k=xEstimate(1);
    d=xEstimate(2);
    concentration=para;%different concentrations

    s=zeros(3,1);
    s(1)=a*y(1)-y(1)*(((a-b)*(concentration(1)/k).^d)/((concentration(1)/k).^d-(b/(a)))); 
    s(2)=a*y(2)-y(2)*(((a-b)*(concentration(2)/k).^d)/((concentration(2)/k).^d-(b/(a)))); 
    s(3)=a*y(3)-y(3)*(((a-b)*(concentration(3)/k).^d)/((concentration(3)/k).^d-(b/(a)))); 



   end

output=values;
end

Этот код работает, когда данные не являются NaN, но при отсутствии данных выдает ошибку как:

Функция цели возвращает неопределенные значения в начальной точке.lsqcurvefit не может продолжить.

Что я могу сделать здесь, чтобы решить эту проблему?Должен ли я вводить данные времени и y как массивы ячеек?
Если это так, я не совсем понимаю, как изменить код для работы с массивами ячеек.
Любая помощь приветствуется.

1 Ответ

1 голос
/ 26 апреля 2019

Что вы можете сделать, это использовать lsqnonlin вместо lsqcurvefit.Это позволяет вам иметь больше гибкости в векторе ошибок, который вы хотите минимизировать.

Вы можете определить погрешность для каждого измерения для каждой концентрации, а затем объединить их с большим вектором ошибок, который вы хотите минимизировать.Простым примером является, например, модель синуса с различной амплитудой, частотой и фазой.Скажем, мы знаем фазу и хотим найти амплитуду и частоту.

Модель: 

function y = model(t, A, w, phi)
    y = A.*sin(w.*t+phi); 
end

Функция ошибки для процесса подгонки принимает измеренные данные и известные параметры.Определите y_estimated для определенного набора параметров (заданного lsqnonlin) и определите погрешность с помощью измеренных значений y_meas.Сделайте это для всех различных концентраций и объедините их в один вектор ошибок.Поскольку некоторые значения в y_meas являются NaN, и вы хотите их опустить, удалите их из вектора ошибок.

function err = errorFun(params, y_meas, t, phi)

    % get w and phi from params
    A = params(1:3);
    w = params(4:6);

    % simulate model with params to estimate
    yest = model(t, A, w, phi);

    % determine error vector
    err = y_meas-yest; 
    err = err(:);           % make one vector
    err(isnan(err)) = [];   % remove NaNs
end

Пример: 

t = 0:0.5:4*pi;
A = rand(3,1);
w =  rand(3,1);
phi = rand(3,1);

y_true = A.*sin(w.*t+phi); % three sinusoids

% measured data
y_meas = y_true;
y_meas(randi([1 numel(y_meas)], 10,1)) = NaN; % set some values to NaN

% optimize
% p = [A;w];
p0 = [A;w;]+0.1; % small deviation from initial parameters, for the example
[p_estimated,a] = lsqnonlin(@(p) errorFun(p,y_meas,t,phi), p0); 

A_est = p_estimated(1:3);
w_est = p_estimated(4:6);

disp([A A_est w w_est])
...