R中的简单曲线拟合
simple curve fitting in R
我正在尝试找到适合我的数据。但到目前为止还没有运气。
尝试了对数的,与 drc 包不同的那些..但我相信一定有更好的我只是不知道类型。
另一方面 - 如果您能提供有关如何进行一般曲线搜索的建议,我将不胜感激。
library(drc)
df<-structure(list(x = c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55), y = c(0.1066, -0.6204, -0.2028, 0.2621, 0.4083,
0.4497, 0.6343, 0.7762, 0.8809, 1.0029, 0.8089, 0.7845, 0.8009,
0.9319, 0.9414, 0.9505, 0.9323, 1.0321, 0.9381, 0.8975, 1.0929,
1.0236, 0.9589, 1.0644, 1.0411, 1.0763, 0.9679, 1.003, 1.142,
1.1049, 1.2868, 1.1569, 1.1952, 1.0802, 1.2125, 1.3765, 1.263,
1.2507, 1.2125, 1.2207, 1.2836, 1.3352, 1.1311, 1.2321, 1.4277,
1.1645), w = c(898, 20566, 3011, 1364, 1520, 2376, 1923, 1934,
1366, 1010, 380, 421, 283, 262, 227, 173, 118, 113, 95, 69, 123,
70, 80, 82, 68, 83, 76, 94, 101, 97, 115, 79, 98, 84, 92, 121,
97, 102, 93, 92, 101, 74, 124, 64, 52, 63)), row.names = c(NA,
-46L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(`47` = 47L), class = "omit"))
fit <- drm(data = df,y ~ x,fct=LL.4(), weights = w)
plot(fit)
基本思路是了解所选函数的执行情况。采用您知道的功能(例如逻辑)并修改它。或者(甚至更好)查阅文献,看看人们在您的特定领域中使用了哪些功能。然后创建一个用户定义的模型,使用它来理解参数,定义好的起始值然后拟合它。
她是一个用户定义函数的快速和肮脏的例子(包 growthrates)。它肯定可以与 drc 类似地制作。
library("growthrates")
grow_userdefined <- function (time, parms) {
with(as.list(parms), {
y <- (K * y0)/(y0 + (K - y0) * exp(-mumax * time)) + shift
return(as.matrix(data.frame(time = time, y = y)))
})
}
fit <- fit_growthmodel(FUN=grow_userdefined,
p = c(y0 = -1, K = 1, mumax = 0.1, shift = 1),
time = df$x, y = df$y)
plot(fit)
summary(fit)
当然可以做得更好。由于我们在开始时没有指数开始,例如可以从一个简单的饱和函数而不是逻辑函数开始,例如类似 Monod 的东西。如前所述,首选方法是使用与应用程序域相关的函数。
1) 如果我们忽略权重,那么 y = a + b * x + c/x^2 似乎适合并且在系数中是线性的所以很容易适合。这似乎是向上倾斜的,所以我们从一条线开始,但随后我们需要抑制它,所以我们添加了一个倒数项。基于残差平方和的二次倒数比简单的倒数稍微好一点,所以我们改用它。
fm <- lm(y ~ x + I(1 / x^2), df)
coef(summary(fm))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.053856e+00 0.116960752 9.010341 1.849238e-11
## x 4.863077e-03 0.002718613 1.788808 8.069195e-02
## I(1/x^2) -1.460443e+02 16.518887452 -8.841049 3.160306e-11
x 项的系数在 5% 的水平上不显着——上面 table 中的 p 值是 8%——所以我们可以删除它,它几乎和给定的一样适合只有两个参数的模型。在下面的图中,具有 3 个参数的 fm 拟合是实线,而具有 2 个参数的 fm2 拟合是虚线。
fm2 <- lm(y ~ I(1 / x^2), df)
plot(y ~ x, df)
lines(fitted(fm) ~ x, df)
lines(fitted(fm2) ~ x, df, lty = 2)
2) 另一种方法是使用两条直线。这仍然是连续的,但在过渡点有一个不可微分的点。该模型有 4 个参数,每条线的截距和斜率。下面我们确实使用了权重。它具有基于数据外观的明显动机的优点。两条线相交处的断点可能具有重要意义,因为它是较高倾斜度的初始增长与较低倾斜度的后续增长之间的过渡点。
# starting values use lines fitted to 1st ten and last 10 points
fm_1 <- lm(y ~ x, df, subset = 1:10)
fm_2 <- lm(y ~ x, df, subset = seq(to = nrow(df), length = 10))
st <- list(a = coef(fm_1)[[1]], b = coef(fm_1)[[2]],
c = coef(fm_2)[[1]], d = coef(fm_2)[[2]])
fm3 <- nls(y ~ pmin(a + b * x, c + d * x), df, start = st, weights = w)
# point of transition
X <- with(as.list(coef(fm3)), (a - c) / (d - b)); X
## [1] 16.38465
Y <- with(as.list(coef(fm3)), a + b * X); Y
## [1] 0.8262229
plot(y ~ x, df)
lines(fitted(fm3) ~ x, df)
我正在尝试找到适合我的数据。但到目前为止还没有运气。 尝试了对数的,与 drc 包不同的那些..但我相信一定有更好的我只是不知道类型。 另一方面 - 如果您能提供有关如何进行一般曲线搜索的建议,我将不胜感激。
library(drc)
df<-structure(list(x = c(10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,
52, 53, 54, 55), y = c(0.1066, -0.6204, -0.2028, 0.2621, 0.4083,
0.4497, 0.6343, 0.7762, 0.8809, 1.0029, 0.8089, 0.7845, 0.8009,
0.9319, 0.9414, 0.9505, 0.9323, 1.0321, 0.9381, 0.8975, 1.0929,
1.0236, 0.9589, 1.0644, 1.0411, 1.0763, 0.9679, 1.003, 1.142,
1.1049, 1.2868, 1.1569, 1.1952, 1.0802, 1.2125, 1.3765, 1.263,
1.2507, 1.2125, 1.2207, 1.2836, 1.3352, 1.1311, 1.2321, 1.4277,
1.1645), w = c(898, 20566, 3011, 1364, 1520, 2376, 1923, 1934,
1366, 1010, 380, 421, 283, 262, 227, 173, 118, 113, 95, 69, 123,
70, 80, 82, 68, 83, 76, 94, 101, 97, 115, 79, 98, 84, 92, 121,
97, 102, 93, 92, 101, 74, 124, 64, 52, 63)), row.names = c(NA,
-46L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(`47` = 47L), class = "omit"))
fit <- drm(data = df,y ~ x,fct=LL.4(), weights = w)
plot(fit)
基本思路是了解所选函数的执行情况。采用您知道的功能(例如逻辑)并修改它。或者(甚至更好)查阅文献,看看人们在您的特定领域中使用了哪些功能。然后创建一个用户定义的模型,使用它来理解参数,定义好的起始值然后拟合它。
她是一个用户定义函数的快速和肮脏的例子(包 growthrates)。它肯定可以与 drc 类似地制作。
library("growthrates")
grow_userdefined <- function (time, parms) {
with(as.list(parms), {
y <- (K * y0)/(y0 + (K - y0) * exp(-mumax * time)) + shift
return(as.matrix(data.frame(time = time, y = y)))
})
}
fit <- fit_growthmodel(FUN=grow_userdefined,
p = c(y0 = -1, K = 1, mumax = 0.1, shift = 1),
time = df$x, y = df$y)
plot(fit)
summary(fit)
当然可以做得更好。由于我们在开始时没有指数开始,例如可以从一个简单的饱和函数而不是逻辑函数开始,例如类似 Monod 的东西。如前所述,首选方法是使用与应用程序域相关的函数。
1) 如果我们忽略权重,那么 y = a + b * x + c/x^2 似乎适合并且在系数中是线性的所以很容易适合。这似乎是向上倾斜的,所以我们从一条线开始,但随后我们需要抑制它,所以我们添加了一个倒数项。基于残差平方和的二次倒数比简单的倒数稍微好一点,所以我们改用它。
fm <- lm(y ~ x + I(1 / x^2), df)
coef(summary(fm))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.053856e+00 0.116960752 9.010341 1.849238e-11
## x 4.863077e-03 0.002718613 1.788808 8.069195e-02
## I(1/x^2) -1.460443e+02 16.518887452 -8.841049 3.160306e-11
x 项的系数在 5% 的水平上不显着——上面 table 中的 p 值是 8%——所以我们可以删除它,它几乎和给定的一样适合只有两个参数的模型。在下面的图中,具有 3 个参数的 fm 拟合是实线,而具有 2 个参数的 fm2 拟合是虚线。
fm2 <- lm(y ~ I(1 / x^2), df)
plot(y ~ x, df)
lines(fitted(fm) ~ x, df)
lines(fitted(fm2) ~ x, df, lty = 2)
2) 另一种方法是使用两条直线。这仍然是连续的,但在过渡点有一个不可微分的点。该模型有 4 个参数,每条线的截距和斜率。下面我们确实使用了权重。它具有基于数据外观的明显动机的优点。两条线相交处的断点可能具有重要意义,因为它是较高倾斜度的初始增长与较低倾斜度的后续增长之间的过渡点。
# starting values use lines fitted to 1st ten and last 10 points
fm_1 <- lm(y ~ x, df, subset = 1:10)
fm_2 <- lm(y ~ x, df, subset = seq(to = nrow(df), length = 10))
st <- list(a = coef(fm_1)[[1]], b = coef(fm_1)[[2]],
c = coef(fm_2)[[1]], d = coef(fm_2)[[2]])
fm3 <- nls(y ~ pmin(a + b * x, c + d * x), df, start = st, weights = w)
# point of transition
X <- with(as.list(coef(fm3)), (a - c) / (d - b)); X
## [1] 16.38465
Y <- with(as.list(coef(fm3)), a + b * X); Y
## [1] 0.8262229
plot(y ~ x, df)
lines(fitted(fm3) ~ x, df)