同时拟合数据以在matlab中估计具有NaN数据的公共参数

Question

对于三种浓度（c），我有超过4个时间点（t）测量的数据（y）。 所以数据会像：

在c1：y11，y12，y13，y14（在4个时间点测量）

在c2：y21，y12，y23，y24

在c3：y31，y32，y33，y34

我想要做的是通过同时拟合不同的连接处的所有这些数据测量来估计两个参数c和d。 但是，其中一些值是NaN。所以，例如，

你好c2：y21，y12，什么，什么

在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无法继续。

我可以在这做什么来解决这个问题？我应该将时间和数据作为单元格数组输入吗？ 如果是这样，我不太明白如何更改代码以使用单元格数组。 任何帮助表示赞赏。

Answer 1

你可以做的是使用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])