【问题标题】:NLS optimization with constraints : equal area under the curve带约束的 NLS 优化:曲线下面积相等
【发布时间】:2020-01-11 03:30:28
【问题描述】:

我有一条中间有扰动的曲线(在 t=9 处),这导致我的数据出现下降(t9)。我想拟合一个指数函数并添加一个约束,即两者(下降和上升)的面积相等。

见图:

我可以使用 optim 拟合曲线,但我无法弄清楚约束。这应该是这样的:

其中 f(x) 是指数函数。

我已经尝试过 constrOptim,但我也愿意使用其他求解器。


y <-c(170, 160, 145, 127, 117, 74, 76, 78, 101, 115, 120, 70, 64, 65)
t <- seq(1,14,1)

# starting values:
lm <-lm(log(y) ~ log(t))

# Exp. Least-Squares minimization:
func <-function(pars) {
  a <- pars["a"]
  b <- pars["b"]
  fitted <- a*exp(b*t)
  sum((y-fitted)^2)  
} 

a <-lm$coefficients[[1]]
b <-lm$coefficients[[2]]
c <- 

result <- optim(c(a=a, b=b), func)

# final parameters:
a <- result$par["a"]
b <- result$par["b"]


# predict values:
pred <- a*exp(b*t)
dat = data.frame(y=y, t=t, pred=pred)

library(ggplot2)
ggplot(dat, aes(x=t, y=y)) +
  geom_line() +
  geom_line(data=dat, aes(x=t, y=pred), color='blue')

编辑:

我知道我需要将约束添加到上面的优化中。像这样:

i = 6:12

result <- optim(c(a=a, b=b), func, sum(y[i]-a*exp(b*t[i])=0)

但这似乎不起作用。 optim 函数不允许这种约束。

【问题讨论】:

    标签: r optimization constraints


    【解决方案1】:

    我可以建议将干扰建模为正弦波吗?

    fit <- nls(y ~ a * exp(b * t) + ifelse(d*(t - m) < 0 | d*(t - m) > 2 * pi, 0 , c * sin(d*(t - m))), 
               start = list(a = 170, b = -0.1, c = -20, d = 1, m = 5))
    summary(fit)
    
    dat <- data.frame(y=y, t=t)
    preddat <- data.frame(t = seq(min(t), max(t), length.out = 100))
    preddat$y <- predict(fit, preddat)
    preddat$y1 <- coef(fit)[["a"]] * exp(coef(fit)[["b"]] * preddat$t)
    
    
    library(ggplot2)
    ggplot(dat, aes(x=t, y=y)) +
      geom_line() +
      geom_line(data=preddat, color='blue') +
      geom_line(data=preddat, aes(y = y1), color='red', linetype = 2)
    

    我确信有信号处理背景的人可以想出比丑陋的ifelse hack 更好的东西。

    【讨论】:

    • 嗨@Roland,我们的目标实际上是平衡低迷和好转,看看没有干扰的曲线会是什么样子。
    • 嗨@Roland,感谢您的回复。不幸的是,我不能保证干扰周围的形状总是像这个例子一样漂亮和对称。所以这行不通。
    • 这不是我的问题的完美解决方案,但我认为这是一个好主意,我实际上将它用于其他事情。所以谢谢。
    猜你喜欢
    • 2015-03-07
    • 1970-01-01
    • 1970-01-01
    • 2021-09-12
    • 1970-01-01
    • 1970-01-01
    • 2021-08-07
    • 1970-01-01
    相关资源
    最近更新 更多