微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

在 R 中使用 `fft()` 的带有卷积的非线性模型的错误曲线图

如何解决在 R 中使用 `fft()` 的带有卷积的非线性模型的错误曲线图

我有一个适合以下数据的非线性模型 ExpDec2RC

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)

enter image description here

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()

enter image description here

使用简单模型,拟合值和曲线完美重叠。

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)

enter image description here

ggplot(data,fit2),colour = "cyan") +
  scale_y_log10() +
  stat_function(fun = function (x) ExpDec2(t = x,size= 2) +
theme_bw()

enter image description here

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。