如何解决在 R 中使用 `fft()` 的带有卷积的非线性模型的错误曲线图
data <- read.table("https://gitlab.com/-/snippets/1999226/raw",sep = "\t")
colnames(data) <- c("irf","decay")
data$t <- 1:nrow(data)
data <- data[367:1279,]
data$t <- data$t - 366
模型使用fft()
进行卷积,如下所示。
ExpDec2RC <- function(t,irf,i0,shift,alpha1,tau1,alpha2,tau2) {
irfbins <- seq_along(t) - 1 # start from 0
irfcounts <- irf
n <- length(irf)
irfshifted <- ((1 - shift + floor(shift)) *
irfcounts[((((irfbins - floor(shift) - 1) %% n) + n) %% n) + 1]) +
((shift - floor(shift)) *
irfcounts[((((irfbins - ceiling(shift) - 1) %% n) + n) %% n) + 1])
irfnormalised <- irfshifted / sum(irfshifted)
ft_ymodel = fft((alpha1 * exp(-t / tau1)) + (alpha2 * exp(-t / tau2)))
ft_irfnormalised = fft(irfnormalised)
ftcp = ft_ymodel*ft_irfnormalised
i0 + Re(fft(z = ftcp,inverse = TRUE)/length(ftcp))
}
通过一组拟合参数,我可以得到如下拟合值。
data$fit <- ExpDec2RC(t = data$t,irf = data$irf,i0 = 11.7764066,shift = 2.00000973,alpha1 = 58414.9896,tau1 = 1.85542518,alpha2 = 1042.19136,tau2 = 75.4146475)
但是,当我使用 curve()
绘制相同的曲线时,我得到的曲线与拟合值完全不同。为什么会这样以及如何使用 curve()
获得正确的拟合值?
plot(data$t,data$decay,col="red",pch="o",log="y")
points(data$t,data$irf,col="blue",pch="o")
# plot fitted values
lines(data$t,data$fit,col="cyan",lwd = 5)
# plot with curve()
curve(ExpDec2RC(t = x,tau2 = 75.4146475),col = "yellow",add = TRUE,lwd = 5)
legend(600,9000,legend=c("Fitted values","With curve()"),col=c("cyan","yellow"),lty = 1,lwd = 5)
stat_function()
中的 ggplot2
也是如此
ggplot(data,aes(t,decay)) +
geom_point(alpha =0.25,colour = "red") +
geom_point(aes(t,irf),alpha =0.25,colour = "blue") +
geom_point(aes(t,fit),colour = "cyan") +
scale_y_log10() +
stat_function(fun = function (x) ExpDec2RC(t = x,colour = "yellow",inherit.aes = FALSE,alpha = 0.5,size= 2) +
theme_bw()
使用简单模型,拟合值和曲线完美重叠。
ExpDec2 <- function(t,tau2) {
i0 + (alpha1 * exp(-t / tau1)) + (alpha2 * exp(-t / tau2))
}
data$fit2 <- ExpDec2(t = data$t,i0 = 11.63431,alpha1 = 6618.36,tau1 = 8.475745201,alpha2 = 644.1506,tau2 = 103.1719055)
plot(data$t,log="y")
# plot fitted values
lines(data$t,data$fit2,lwd = 8)
# plot with curve()
curve(ExpDec2(t = x,tau2 = 103.1719055),lwd = 3)
legend(600,lwd = 5)
ggplot(data,fit2),colour = "cyan") +
scale_y_log10() +
stat_function(fun = function (x) ExpDec2(t = x,size= 2) +
theme_bw()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。