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

【问题标题】：Simultaneously fitting to data to estimate a common parameter with NaN data in matlab在matlab中同时拟合数据以估计具有NaN数据的公共参数
【发布时间】：2023-03-15 02:44:02
【问题描述】：

我有超过 4 个时间点 (t) 测量的数据 (y)，用于三个浓度 (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 数据输入为元胞数组吗？
如果是这样，我不太明白如何更改代码以使用元胞数组。
任何帮助表示赞赏。

【问题讨论】：

这完全取决于您希望如何适应 NaN...您要省略曲线的那些部分吗？您想从非 NaN 值进行插值/外推吗？
@Wolfie 我想从曲线中省略这些。我不想在估算过程中使用这些。

标签： matlab optimization missing-data

【解决方案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])

【讨论】：