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软件包,但是不知道如何在特定问题上使用它们.
I have a 20*10 ruw and 20*10 quw matrix, I now need to generate a yuw matrix which adheres to the constraints. I am coding in R and am familiar with the lpsolve and optimx packages, but don't know how to use them for this particular question.
推荐答案
由于Quw
和ruw
都是数据,因此yuw
决策变量中的所有约束和目标都是线性的.结果,这是一个线性编程问题,lpSolve
程序包可以解决.
Because Quw
and ruw
are both data, all constraints as well as the objective are linear in the yuw
decision variables. As a result, this is a linear programming problem that can be approached by the lpSolve
package.
为了稍微抽象一下,我们假设R=20
和C=10
描述输入矩阵的维数.然后有R*C
个决策变量,我们可以将它们分配为y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC
顺序,读取变量矩阵的列.
To abstract this out a bit, let's assume R=20
and C=10
describe the dimensions of the input matrices. Then there are R*C
decision variables and we can assign them order y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC
, reading down the columns of the matrix of variables.
目标中每个yuw
变量的系数为-Quw
;请注意,求和中的Quw*ruw
项是常量(也不受我们为决策变量选择的值的影响),因此未输入到线性规划求解器中.有趣的是,这意味着ruw
实际上对优化模型解决方案没有影响.
The coefficient of each yuw
variable in the objective is -Quw
; note that the Quw*ruw
terms in the summation are constants (aka not affected by the values we select for the decision variables) and therefore not inputted to the linear programming solver. Interestingly, this means that ruw
actually has no effect on the optimization model solution.
第一个R*(C-1)
约束与yuw >= yu,w+1
约束相对应,接下来的(R-1)*C
约束与yuw >= yu-1,w
约束相对应.最后两个约束对应于y20,1 >= 100
和y1,10 >= 0
约束.
The first R*(C-1)
constraints correspond to the yuw >= yu,w+1
constraints, and the next (R-1)*C
constraints correspond to the yuw >= yu-1,w
constraints. The final two constraints correspond to the y20,1 >= 100
and y1,10 >= 0
constraints.
我们可以使用以下R命令将此模型输入到lpsolve
包中,以简单的Q矩阵作为输入,其中每个条目为-1(结果解决方案应将除左下角以外的所有决策变量都设置为0 ,应该是100):
We can input this model into the lpsolve
package with the following R command, taking as input a simple Q matrix where each entry is -1 (the resulting solution should have all decision variables set to 0 except the bottom left corner, which should be 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))
我们现在可以访问模型解决方案:
We can now access the model solution:
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退出(您可以通过运行mod
或mod$status
看到状态).
Note that your model could easily become infeasible if Quw
changed (for instance if we filled it with 1 instead of -1). In these cases the model will exit with status 3 (you can see this by running mod
or mod$status
).
这篇关于R:等渗回归最小化的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!