R:等渗回归最小化

R: Isotonic regression Minimisation

我想最小化以下等式:

F=SUM{u 1:20}sum{w 1:10}   Quw(ruw-yuw)

具有以下限制条件:

yuw >= yu,w+1
yuw >= yu-1,w
y20,0 >= 100
y0,10 >= 0

我有一个 20*10 ruw 和 20*10 quw 矩阵,我现在需要生成一个遵守约束的 yuw 矩阵。我在 R 中编码并且熟悉 lpsolve 和 optimx 包,但不知道如何将它们用于这个特定问题。

因为Quwruw都是数据,所以所有约束以及objective在yuw决策变量中都是线性的。因此,这是一个可以通过 lpSolve 包解决的线性规划问题。

为了抽象一点,我们假设 R=20C=10 描述了输入矩阵的维度。然后有 R*C 个决策变量,我们可以为它们分配顺序 y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC,读取变量矩阵的列。

objective中每个yuw变量的系数为-Quw;请注意,求和中的 Quw*ruw 项是常数(也不受我们 select 的决策变量值的影响),因此不会输入到线性规划求解器。有趣的是,这意味着 ruw 实际上对优化模型解决方案没有影响。

R*(C-1)个约束对应yuw >= yu,w+1个约束,接下来的(R-1)*C个约束对应yuw >= yu-1,w个约束。最后两个约束对应于 y20,1 >= 100y1,10 >= 0 约束。

我们可以使用以下 R 命令将此模型输入到 lpsolve 包中,将一个简单的 Q 矩阵作为输入,其中每个条目为 -1(生成的解决方案应将所有决策变量设置为 0,除了左下角,应该是100):

# Sample data
Quw <- matrix(-1, nrow=20, ncol=10)
R <- nrow(Quw)
C <- ncol(Quw)

# Build constraint matrix
part1 <- matrix(0, nrow=R*(C-1), ncol=R*C)
part1[cbind(1:(R*C-R), 1:(R*C-R))] <- 1
part1[cbind(1:(R*C-R), (R+1):(R*C))] <- -1
part2 <- matrix(0, nrow=(R-1)*C, ncol=R*C)
pos2 <- as.vector(sapply(2:R, function(r) r+seq(0, R*(C-1), by=R)))
part2[cbind(1:nrow(part2), pos2)] <- 1
part2[cbind(1:nrow(part2), pos2-1)] <- -1
part3 <- rep(0, R*C)
part3[R] <- 1
part4 <- rep(0, R*C)
part4[(C-1)*R + 1] <- 1
const.mat <- rbind(part1, part2, part3, part4)

library(lpSolve)
mod <- lp(direction = "min",
          objective.in = as.vector(-Quw),
          const.mat = const.mat,
          const.dir = rep(">=", nrow(const.mat)),
          const.rhs = c(rep(0, nrow(const.mat)-2), 100, 0))

我们现在可以访问模型解决方案:

mod
# Success: the objective function is 100
matrix(mod$solution, nrow=R)
#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#  [1,]    0    0    0    0    0    0    0    0    0     0
#  [2,]    0    0    0    0    0    0    0    0    0     0
#  [3,]    0    0    0    0    0    0    0    0    0     0
#  [4,]    0    0    0    0    0    0    0    0    0     0
#  [5,]    0    0    0    0    0    0    0    0    0     0
#  [6,]    0    0    0    0    0    0    0    0    0     0
#  [7,]    0    0    0    0    0    0    0    0    0     0
#  [8,]    0    0    0    0    0    0    0    0    0     0
#  [9,]    0    0    0    0    0    0    0    0    0     0
# [10,]    0    0    0    0    0    0    0    0    0     0
# [11,]    0    0    0    0    0    0    0    0    0     0
# [12,]    0    0    0    0    0    0    0    0    0     0
# [13,]    0    0    0    0    0    0    0    0    0     0
# [14,]    0    0    0    0    0    0    0    0    0     0
# [15,]    0    0    0    0    0    0    0    0    0     0
# [16,]    0    0    0    0    0    0    0    0    0     0
# [17,]    0    0    0    0    0    0    0    0    0     0
# [18,]    0    0    0    0    0    0    0    0    0     0
# [19,]    0    0    0    0    0    0    0    0    0     0
# [20,]  100    0    0    0    0    0    0    0    0     0

请注意,如果 Quw 更改(例如,如果我们用 1 而不是 -1 填充它),您的模型很容易变得不可行。在这些情况下,模型将以状态 3 退出(您可以通过 运行 modmod$status 查看)。