Оптимизация оценки нелинейных параметров Ленгмюра в R

Question

Меня интересует использование nls, чтобы помочь в установке уравнения Лангмюра Y =(Qmax*k*X)/(1+(k*X)), аналогично тому, что было сделано в этой статье Fitting Non-linear Langmuir Isotherm in R. Параметр интересующего меня уравнения равен , который соответствует горизонтальной асимптоте (зеленая линия) приведенных данных сорбции ниже. Есть ли более надежный подход, отличный от nls, или способ улучшить мое использование nls, которое я мог бы использовать, чтобы получить значение как можно ближе к визуальной асимптоте (зеленая линия) около 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. С учетом сказанного было бы весьма признательно, что любые технические рекомендации, которые могут помочь мне улучшить мою технику выше, будут весьма полезны.

Answer 1

Не уверен, почему вы ожидаете быть, что низкий

Я переписал вашу зависимость в простейшей форме, удаляя умножение и заменить его с добавлением (а => 1/к) путем деления как числитель и знаменатель на 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 

Graph

UPDATE

В основном, слишком много точек при низких 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") 

я получить и a

[1] 4121.114 
[1] 74.89386 

следующим графиком