具有顺序约束的 LS 回归答案

【问题标题】：LS regression with order constraint具有顺序约束的 LS 回归
【发布时间】：2017-02-26 14:08:03
【问题描述】：

我是 R 新手，非常感谢您的帮助。我有一个约束的优化问题。尽管有多种方法可以解决 R 中的优化问题，但我无法使用需要应用的约束正确表达我的问题。

假设我有以下三类数据：

A<-c(99.1,  96.5,   94.4,   92.7,   91.5,   91.3,   91.4,   90.1,   87.1,   82.6,   76.4)
B<-c(146.4, 140.2,  133.6,  126.5,  118.7,  109.4,  101.2,  101.8,  103.7,  102.5,  98.3)
C<-c(237.5, 213.9,  191,    168.9,  147.4,  124.9,  108.3,  95.7,   84.4,   73.5,   63)
t<-seq(1:11)
DT<-cbind.data.frame(t,A,B,C)

我想将指数函数 y(t) 拟合到每个类别中的数据点（最小化平方误差），以便 y(t)_c > y( t)_b > y(t)_a > 0 代表选定的 t [1;15]

【问题讨论】：

标签： r optimization regression least-squares simplex

【解决方案1】：

我不会尝试将约束合并到回归中。只需构建三个独立的回归：

fit_loga <- lm(y ~ log(A) + t, data=DT) 
fit_logb <- lm(y ~ log(B) + t, data=DT) 
fit_logc <- lm(y ~ log(C) + t, data=DT)

fit_a <- exp(A)
fit_b <- exp(B)
fit_c <- exp(C)

然后验证它们是否满足该范围（或至少在每个整数数据点）上的约束：(fit_c > fit_b) & (fit_b > fit_a)。只有没有，我们才会担心。比如把 t 中的一些其他术语加入到模型中，可能是exp(t)、I(t^2)、poly(t, <order>)...

编辑：我不知道solnp 包。

【讨论】：

感谢您的回复 smci！我想像你建议的那样分两步完成它，但是在这种情况下，当 C 明显与 B 和 A 相交时，你会怎么做呢？恐怕我需要对所有曲线保持相同的曲线形式假设（这里：指数）。因此，一个参数会影响另一个参数，这需要同时对它们进行优化。我完全同意您应该在对数空间中寻找线性曲线的参数，而不是尝试拟合指数曲线。
哦，废话。好吧，如果它们相交，那么我们如何应用该约束？不用做一些笨拙的事情，比如使用 min/max 进行夹紧？
我认为这个问题可以通过像 {solnp()} 这样可以处理不等式约束的优化函数来解决。但是，当直接应用于最小化给定约束的平方误差时，它不会收敛或导致不可行的解决方案
@Walle：我不知道solnp。建议您将代码和不可行的错误添加到问题详细信息和标题中。此外，这将在CrossValidated.SE
原来是solnp() 的格式问题，因为该函数需要估计所有参数，作为单个向量的输入。因此，我必须将必须的矩阵转换为向量并且它起作用了。从运行该函数时收到的错误消息中并不清楚，即non-numeric argument to binary operator 用于不可行的约束（太紧）。不知道我是应该发布代码还是关闭这个主题。

【解决方案2】：

我认为使用solnp 得到的错误消息主要是指约束不足。此外，正如文档中所述，需要将所有参数放在一个向量中。在对代码进行适当调整后，我能够直接实现y(t)_c > y(t)_b > y(t)_a > 0 的约束，而无需更改问题。最方便的方法是对求解器使用矩阵形式。使用上面的数据，我得到以下信息： Results shown here

library(data.table)
library(Rsolnp)

t<-as.vector(10:20)
DT<-cbind.data.frame(A,B,C)
tlogDT<-transpose(log(DT))

# min[log(y)'- Ax-b]
# Arr = [A1 A2 .. An b1 b2 .. bn]

gofn = function(arrin)
{
  arr = cbind(arrin[1:3],arrin[4:6])
  sum(
    (as.matrix(arr[,1])%*%t + arr[,2] - tlogDT)^2
  )
}

nocross=t #defines the times where the curves are not allowed to intersect
ineqfn2=function(arrin)
{
  #constrains:
  # 0<f_a(t)<f_b(t)<f_c(t), for some t,
  arr = cbind(arrin[1:3],arrin[4:6])
  nextarr=as.matrix(rbind(rep(0,2),arr[1:(length(arr[,1])-1),]))
  ineqmat=as.matrix(arr[,1])%*%nocross+arr[,2]-nextarr[,1]%*%nocross-nextarr[,2]
  as.vector(t(ineqmat))
}

#lines should be aligned with the first startvalue
eqfn = function(arrin)
{  
  arr = cbind(arrin[1:3],arrin[4:6])
  arr[,1]*t[1]+arr[,2]-tlogDT[,1]
}
#start values:
An=c(1,1,1)
bn=tlogDT[,1]
vstart=c(An,bn)

r <- solnp(
  vstart,  gofn,

  eqfun = eqfn, eqB= c(0,0,0),

  ineqfun = ineqfn2,
  ineqLB = rep(0,length(DT[1,])*length(nocross)),
  ineqUB = rep(5000,length(DT[1,])*length(nocross))
)

r$pars[1]
line1 = exp(r$pars[4]+r$pars[1]*t)
line2 = exp(r$pars[5]+r$pars[2]*t)
line3 = exp(r$pars[6]+r$pars[3]*t)

plot(t, DT[,3],log = "y")
points(t, DT[,2],col="green")
points(t, DT[,1],col="blue")

lines(t, line1,lwd=2, col = "blue",  xlab = "Time (s)", ylab = "Counts")
lines(t, line2,lwd=2, col = "green",  xlab = "Time (s)", ylab = "Counts")
lines(t, line3,lwd=2, col = "black",  xlab = "Time (s)", ylab = "Counts")

【讨论】：