我希望人们在尝试使用 R 进行非线性参数估计时，每浪费一小时我就有一美元。

这是您的问题的解决方案以及通过 delta 方法计算的估计标准偏差，解决方案图在上面。

y0 85.557909  3.0989e-01

a1 125.20943 1.3766e+01

a2 1394.155 4.4952e+03

b1 0.062640298 4.1774e-03

b2 0.43392314 2.9936e-01

我使用了与我在那里发布的相同的 AD 模型生成器代码，但针对您的问题稍作修改。正如我所说的，这是一种解决这类问题的通用技术。

如何为非线性最小二乘拟合选择初始值

这是针对您的问题的 AD 模型生成器代码。

DATA_SECTION
  init_int n
  int mid
 !! mid=75;
  init_matrix data(1,n,1,3)
  vector t(1,n)
  vector P(1,n)
 !! t=column(data,2);
 !! P=column(data,3);   //use column 3
  number tmax
  number Pmax
 !! tmax=max(t);
 !! Pmax=max(P);
 !! t/=tmax;
 !! P/=Pmax;

PARAMETER_SECTION
  init_number L1(3)     //(3) means estimate in phase 3
  init_number Lmid(3)
  init_number Ln(3)

  vector L(1,3)
  init_number log_b1       // estimate in phase 1
  init_number log_diff    // estimate in phase 2 
  sdreport_number b1
  sdreport_number b2
  matrix M(1,3,1,3);
  objective_function_value f
  sdreport_vector v(1,3)
  sdreport_number real_y0
  sdreport_number real_a1
  sdreport_number real_a2
  sdreport_number real_b1
  sdreport_number real_b2
  vector pred(1,n);
PROCEDURE_SECTION
  L(1)=L1;
  L(2)=Lmid;
  L(3)=Ln;
  b1=exp(log_b1);    // this parameterization ensures that b1<b2
  b2=b1+exp(log_diff);  
  M(1,1)=exp(-b1*t(1));
  M(1,2)=exp(-b2*t(1));
  M(1,3)=1;
  M(2,1)=exp(-b1*t(mid));
  M(2,2)=exp(-b2*t(mid));
  M(2,3)=1;
  M(3,1)=exp(-b1*t(n));
  M(3,2)=exp(-b2*t(n));
  M(3,3)=1;

  v=solve(M,L);  // solve for standard parameters 
                 // v is vector corresponding to the parameters b_1 b_2 P_0

  pred=v(3)-v(1)*exp(-b1*t)-v(2)*exp(-b2*t);
  if (current_phase()<4)
    f+=norm2(P-pred);
  else   // use concentrated likelihood so that Hessian is correct
    f+=0.5*n*log(norm2(P-pred)); //concentrated likelihood

  real_y0=v(3)*Pmax;
  real_a1=v(1)*Pmax;
  real_a2=v(2)*Pmax;
  real_b1=b1/tmax;
  real_b2=b2/tmax;

REPORT_SECTION
  dvar_matrix tmp(1,2,1,n);
  dvar_vector real_t=t*tmax;
  dvar_vector real_pred=real_y0-real_a1*exp(-real_b1*real_t)
    -real_a2*exp(-real_b2*real_t);

  tmp(1)=real_t;
  tmp(2)=real_pred;

  ofstream ofs1("pred");
  ofs1 << trans(tmp)<< endl;

  tmp(2)=P*Pmax;

  ofstream ofs2("obs");
  ofs2 << trans(tmp)<< endl;


  report << "y0 " << setprecision(8) << real_y0 << endl;
  report << "a1 " << setprecision(8) << real_a1 << endl;
  report << "a2 " << setprecision(8) << real_a2 << endl;
  report << "b1 " << setprecision(8) << real_b1 << endl;
  report << "b2 " << setprecision(8) << real_b2 << endl;

编辑删除原来的第二点，这是错误的。两条评论：

1. 查看数据，大约在 80 秒和 140 秒时发生了一些奇怪的事情。如果这些是实验性的人工制品，那么拟合复杂的模型就没有意义了。如果这些数据显示出真实的趋势（不是实验伪影），那么您将需要一个正弦模型来适应这些曲折和弯曲。拟合双指数模型假设所有变化都是高斯的，但 80 秒和 140 秒的趋势似乎太大而不能只是随机变化。

2. 以下是GraphPad Prism的结果，与 Dave 使用 ADBuilder（ Prism 文件）得到的结果非常相似：

一个。查看 A2 和 tau2 的置信区间。它们交叉成负值，这是不可能的，而且范围很广。您的数据确实没有定义这些参数。

湾。游程测试的 P 值很小。这意味着数据系统地偏离曲线。一个点是在曲线之上还是之下并不是完全随机的。这与上面的第 1 点一致。

C。如果您尝试将曲线外推至 Time=0，则 Y 值会变得非常负。根据科学背景，这可能表明该模型不是很明智。也许模型应该包含一个 (X - X0) 项，其中 X0 是一个常数（可能是 11.33，您的第一个时间点）？

请注意，在我的讨论中，我将忽略模型明显缺乏拟合并专注于获得估计的问题。

像这样的模型有时可能会出现一些问题（因为参数空间中的脊），但这里似乎没有这样的困难。

和你一样，我先做了一个指数拟合。然后我尝试的第一个两个指数拟合起作用了。

这是我所做的，按顺序：

首先，我尝试为 y0、a1 和 b1 的一些值绘制曲线，直到我发现可以正常工作的东西：

步骤 0：通过反复试验找到单指数模型的一些值

#Leaving out the values I tried that didn't work as well, I got to here:
 plot(X1r_1_green~Time,data1)
 t=10:189
 f1=90-100*exp(-0.05*t)
 lines(t,f1,col=4)

第 1 步：拟合一个指数模型：

fit1exp=nls(X1r_1_green ~ y0-a1*exp(-b1*Time),
                   data=data1,start=c(y0=90,a1=100,b1=.05))
summary(fit1exp)

Formula: X1r_1_green ~ y0 - a1 * exp(-b1 * Time)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
y0 8.542e+01  3.007e-01  284.11   <2e-16 ***
a1 1.447e+02  7.140e+00   20.26   <2e-16 ***
b1 6.804e-02  2.707e-03   25.13   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.271 on 172 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 3.197e-06

第 2 步：将它们缩小到 0，对第二个参数集进行粗略的猜测（我将它们设置在第一个参数集的 15-20% 的范围内，但对此不精确是值得的）：

 fit2exp=nls(X1r_1_green~ y0-a1*exp(-b1*Time)-a2*exp(-b2*Time),
                  data=data1,start=c(y0=80,a1=120,b1=.06,a2=20,b2=.01))

 summary(fit2exp)

Formula: X1r_1_green ~ y0 - a1 * exp(-b1 * Time) - a2 * exp(-b2 * Time)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
y0 8.771e+01  1.205e+01   7.281 1.18e-11 ***
a1 1.509e+02  9.122e+00  16.539  < 2e-16 ***
b1 7.379e-02  7.828e-03   9.427  < 2e-16 ***
a2 5.132e+00  6.504e+00   0.789    0.431    
b2 6.429e-03  3.356e-02   0.192    0.848    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.261 on 170 degrees of freedom

Number of iterations to convergence: 22 
Achieved convergence tolerance: 4.992e-06

我没想到第二个模型会先尝试；通常你不得不对这样的模型进行一些调整（或者更好的是，使用学科领域知识来建立更好的猜测——但我没有这些）。

在与初始拟合相同的图上绘制拟合：

蓝色单指数的原始逐眼拟合，红色的最终双指数 nls 拟合。还有一个绿色的单指数 nls-fit，但它几乎完全被红色曲线重叠 - 即第二个指数几乎没有增加这里的拟合。）

我花了更长的时间来编辑你的帖子，而不是让那个合适的工作，其中大部分只是通过肉眼找到蓝色曲线。

[虽然这似乎工作得很好，但尝试其他起始值并查看您是否已收敛到局部最优值，或者是否卡在一个非常伸展的山脊中并最终与您可能所在的位置相距一段距离，这确实是值得的. 我通常也会尝试更改默认收敛标准。]

如果您意识到您的模型是部分线性的，那么整个任务就会变得非常简单。您只需要非线性参数的起始值，并且由于这些与半衰期有关，因此很容易从图中做出一个不错的猜测：

tau1_ini <- log(2) / 18            
tau2_ini <- log(2) / 180

然后你可以适应：

fit <- nls(X1r_1_green ~ cbind(1, exp(-tau1 * Time), exp(-tau2 * Time)), 
           data = DF, algorithm = "plinear", 
           start = list(tau1 = tau1_ini, tau2 = tau2_ini))

summary(fit)

#Formula: X1r_1_green ~ cbind(1, exp(-tau1 * Time), exp(-tau2 * Time))
#
#Parameters:
#        Estimate Std. Error t value Pr(>|t|)    
#tau1   7.379e-02  7.828e-03   9.427  < 2e-16 ***
#tau2   6.429e-03  3.356e-02   0.192    0.848    
#.lin1  8.771e+01  1.205e+01   7.280 1.18e-11 ***
#.lin2 -1.509e+02  9.122e+00 -16.539  < 2e-16 ***
#.lin3 -5.132e+00  6.505e+00  -0.789    0.431    
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 3.261 on 170 degrees of freedom
#
#Number of iterations to convergence: 19 
#Achieved convergence tolerance: 6.29e-06

这与@Glen_b 的答案相同。

请注意，参数估计是强相关的。您的数据并不真正支持如此复杂的模型。