用于指数衰减函数的Matlab图

Question

我有9组患者的经验数据,数据以这种格式显示

               input = [10      -1      1
                        20      17956   1
                        30      61096   1
                        40      31098   1
                        50      18446   1
                        60      12969   1
                        95      7932    1
                        120     6213    1
                        188     4414    1
                        240     3310    1
                        300     3329    1
                        610     2623    1
                        1200    1953    1
                        1800    1617    1
                        2490    1559    1
                        3000    1561    1
                        3635    1574    1
                        4205    1438    1
                        4788    1448    1
                      ];
                      calibrationfactor_wellcounter =1.841201569;

这里,第一列描述时间值,下一列是浓度.如您所见,浓度会增加一段时间,然后随着时间的增加呈指数下降.

如果我绘制以下特征,我获得以下曲线

我想创建一个代表上面引用的相同行为的脚本.以下是我制定的脚本,其中浓度线性增加直到某个时间段和后果它以指数方式衰减,但是当我绘制此函数我获得线性特征时,请告诉我我的逻辑是否合适

function c_o = Sample_function(td,t_max,a1,a2,a3,b1,b2,b3)
   t =(0: 100 :5000); % time of the sample post injection in mins
   c =(0 : 2275.3 :113765);
   A_max= max(c);%Max value of Concentration (Peak of the curve)

   c_o = zeros(size(t));
   c_o(t>td & t<=t_max) = A_max*(t(t>td & t<=t_max)-td);

   c_o(t>t_max)=(a1*exp(-b1*(t(t>t_max)-t_max)))+(a2*exp(-b2*(t(t>t_max)-t_max)))+(a3*exp(-b3*(t(t>t_max)-t_max)));        

   fprintf('plotting Data ...\n');
   hold on;
   %figure ;
   plot(c_o,'erasemode','background');
   xlabel('time of the sample in minutes ');
   ylabel('Activity of the sample Ba/ml');
   title (' Input function: Activity sample VS time ');

   pause;
end

我得到的数字是 在上面的图中,衰减是线性的而不是指数的,让我知道如何获得3阶衰减这是我为了获得3阶衰变而编写的代码行

c_o(t>t_max)=(a1*exp(-b1*(t(t>t_max)-t_max)))+(a2*exp(-b2*(t(t>t_max)-t_max)))+(a3*exp(-b3*(t(t>t_max)-t_max)));

Answer 1

我已经提出了一个使用Matlab 曲线拟合工具箱功能的解决方案.拟合结果看起来非常好.但是,我发现它在很大程度上取决于参数的起始值的正确选择,因此必须手动仔细选择.

从你的变量开始input,让我们为拟合,时间和浓度定义独立变量和因变量,

t = input(:, 1);
c = input(:, 2);

并绘制它们:

plot(t, c, 'x')
axis([-100 5000 -2000 80000])
xlabel time
ylabel concentration

这些数据将用具有三个部分的函数建模:1)恒定0直到时间td,2)在td和之间线性增加tmax,3)在时间之后减少为三个不同指数的总和tmax.此外,功能是连续的,因此三件必须无缝地配合在一起.该模型作为Matlab函数的实现:

function c = model(t, a1, a2, a3, b1, b2, b3, td, tmax)

c = zeros(size(t));

ind = (t > td) & (t < tmax);
c(ind) = (t(ind) - td) ./ (tmax - td) * (a1 + a2 + a3);

ind = (t >= tmax);
c(ind) = a1 * exp(-b1 * (t(ind) - tmax)) ...
    + a2 * exp(-b2 * (t(ind) - tmax)) + a3 * exp(-b3 * (t(ind) - tmax));

模型参数似乎由曲线拟合工具箱内部处理为由参数名称按字母顺序排序的向量,因此为了避免混淆,我还在此函数的定义中按字母顺序对参数进行了排序.a1to a3和b1to b3分别是三个指数的幅度和反时限常数.

让模型适合数据:

ft = fittype('model(t, a1, a2, a3, b1, b2, b3, td, tmax)', 'independent', 't');
fo = fit(t, c, ft, ...
    'StartPoint', [20000, 20000, 20000, 0.01, 0.01, 0.01, 10, 30], ...
    'Lower', [0, 0, 0, 0, 0, 0, 0, 0])

如前所述,只有当算法得到合适的起始值时,拟合才能正常工作.我在这里选择了的振幅a1到a3数20000,这是大约三分之一的最大的数据中的,用于b1向b30.01对应于约100的时间常数,该数据最大值的时刻,30,对于一个值tmax,并10作为起始恒定时间的粗略估计td.

输出fit:

fo = 

     General model:
     fo(t) = model(t, a1, a2, a3, b1, b2, b3, td, tmax)
     Coefficients (with 95% confidence bounds):
       a1 =        2510  (-2.48e+07, 2.481e+07)
       a2 =   1.044e+04  (-7.393e+09, 7.393e+09)
       a3 =   6.506e+04  (-4.01e+11, 4.01e+11)
       b1 =   0.0001465  (7.005e-05, 0.0002229)
       b2 =     0.01049  (0.006933, 0.01405)
       b3 =     0.09134  (0.08623, 0.09644)
       td =       17.97  (-3.396e+07, 3.396e+07)
       tmax =       26.78  (-6.748e+07, 6.748e+07)

我无法决定这些价值观是否在生理上有意义.由于许多置信区间很大并且实际上包含0,所以估计似乎也没有太明确定义.文档对此并不清楚,但我认为置信区间是非同时的,这意味着有可能巨大的区间仅表示不同参数的估计之间的强相关性.

将数据与拟合模型一起绘制

plot(t, c, 'x')
hold all
ts = 0 : 5000;
plot(ts, model(ts, fo.a1, fo.a2, fo.a3, fo.b1, fo.b2, fo.b3, fo.td, fo.tmax))
axis([-100 5000 -2000 80000])
xlabel time
ylabel concentration

表明合身非常好:

更有趣的初始部分的特写:

注意,真实最大浓度(27,78000)的估计时间和值仅取决于对下面的数据减少部分的拟合,因为线性增加仅由一个数据点表征,其不构成约束.

结果表明,数据不足以获得模型参数的精确估计.您应该考虑增加数据的采样率,特别是时间500,或降低模型的复杂性,例如,仅使用两个指数的总和; 或两者.