在 R 中优化非线性 Langmuir 参数估计答案

【问题标题】：Optimizing Non-Linear Langmuir Parameter Estimation in R在 R 中优化非线性 Langmuir 参数估计
【发布时间】：2016-05-11 19:28:40
【问题描述】：

我有兴趣使用nls 来帮助拟合朗缪尔方程Y =(Qmax*k*X)/(1+(k*X))，类似于本文Fitting Non-linear Langmuir Isotherm in R 中所做的。我感兴趣的方程的参数是Qmax，它对应于下面绘制的吸附数据的水平渐近线（绿线）。除了nls 之外，是否有更强大的方法或改进我对nls 的使用的方法，我可以采用该方法来使Qmax 值尽可能接近Qmax=3200 周围的视觉渐近线（绿线）？

Lang <- nls(formula = Y ~ (Qmax*k*X)/(1+(k*X)),  data = data, start = list(Qmax = 3600, k = 0.015), algorith = "port")

使用以下数据：

     X        Y
1    3.08   84.735
2    5.13  182.832
3    6.67  251.579
4    9.75  460.077
5   16.30  779.350
6   25.10  996.540
7   40.80 1314.739
8   68.90 1929.422
9  111.00 2407.668
10 171.00 3105.850
11 245.00 3129.240
12 300.00 3235.000

我收到了Qmax = 4253.63（红线）- 大约 1000 个单位。使用上限和下限只会导致我设置上限的 Qmax，而更改初始值似乎不会改变结果。这是一个可以通过与我在基础 R 中采用的非线性回归方法不同的方法来解决的挑战，还是首先是统计/数学问题？

Plot of Non-linear Langmuir Isotherm

 summary(Lang)

    Formula: Y ~ (Qmax * k * X)/(1 + (k * X))

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
Qmax 4.254e+03  1.554e+02   27.37 9.80e-11 ***
k    1.209e-02  1.148e-03   10.53 9.87e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 99.14 on 10 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)

我对模型线性化的尝试不太成功：

z <- 1/data
plot(Y~X,z)
abline(lm(Y~X,z))
M <- lm(Y~X,z)

Qmax <- 1/coef(M)[1]
#4319.22

k <- coef(M)[1]/coef(M)[2]
#0.00695

免责声明：这是我的第一篇文章，所以请多多包涵，我对 R 比较陌生。话虽如此，任何可能帮助我改进上述技术的技术建议都将不胜感激。

【问题讨论】：

标签： r data-fitting nls nonlinear-optimization non-linear-regression

【解决方案1】：

不知道为什么你认为Qmax 会那么低

我以最简单的形式重写了您的依赖关系，删除了乘法并将其替换为加法 (a => 1/k)，方法是将提名者和分母都除以 k。结果在我看来很完美。

library(ggplot2)
library(data.table)

dt <- fread("R/Langmuir.dat", sep = " ")

Lang <- nls(formula = Y ~ (Qmax*X)/(a+X),  data = dt, start = list(Qmax = 3600, a = 100.0), algorithm = "port")
q <- summary(Lang)

Qmax <- q$coefficients[1,1]
a    <- q$coefficients[2,1]

f <- function(x, Qmax, a) {
    (Qmax*x)/(a+x)
}

p <- ggplot(data = dt, aes(x = X, y = Y))
p <- p + geom_point()
p <- p + xlab("T") + ylab("Q") + ggtitle("Langmuir Fit")
p <- p + stat_function(fun = function(x) f(x, Qmax=Qmax, a=a))
print(p)

print(Qmax)
print(a)

输出

4253.631
82.68501

图表

更新

基本上，低 X 的点太多，难以获得低 Qmax 的曲线弯曲。设计曲线弯曲的方法是添加权重。例如，如果我在读取数据表后添加权重列：

dt[, W := (as.numeric(N)/12.0)^3]

并使用权重运行nls

Lang <- nls(formula = Y ~ (Qmax*X)/(a+X),  data = dt, start = list(Qmax = 3600, a = 100.0), weights = dt$W, algorithm = "port")

我会得到Qmax 和a

[1] 4121.114
[1] 74.89386

下图

【讨论】：