在 R 中优化非线性 Langmuir 参数估计
Optimizing Non-Linear Langmuir Parameter Estimation in R
我有兴趣使用 nls
来帮助拟合 Langmuir 方程 Y =(Qmax*k*X)/(1+(k*X))
,类似于 post 中所做的。我感兴趣的方程式的参数是 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 并且更改初始值似乎不会改变结果。这是一个可以通过与我在 base R 中采用的不同的非线性回归方法来解决的挑战,还是首先是 statistical/mathematical 问题?
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
免责声明:这是我的第一个 post 所以请多多包涵,我是 R 的新手。话虽如此,如果有任何技术建议可以帮助我改进上述技术,我们将不胜感激.
不确定为什么您认为 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 的曲线弯曲。使曲线弯曲的设计方法是增加重量。例如,如果我在读取数据后添加权重列 table:
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
与下图
我有兴趣使用 nls
来帮助拟合 Langmuir 方程 Y =(Qmax*k*X)/(1+(k*X))
,类似于 post 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 并且更改初始值似乎不会改变结果。这是一个可以通过与我在 base R 中采用的不同的非线性回归方法来解决的挑战,还是首先是 statistical/mathematical 问题?
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
免责声明:这是我的第一个 post 所以请多多包涵,我是 R 的新手。话虽如此,如果有任何技术建议可以帮助我改进上述技术,我们将不胜感激.
不确定为什么您认为 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 的曲线弯曲。使曲线弯曲的设计方法是增加重量。例如,如果我在读取数据后添加权重列 table:
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
与下图