R: Isotonic regression Minimisation -
i want minimize following equation:
f=sum{u 1:20}sum{w 1:10} quw(ruw-yuw) with following constraints:
yuw >= yu,w+1 yuw >= yu-1,w y20,0 >= 100 y0,10 >= 0 i have 20*10 ruw , 20*10 quw matrix, need generate yuw matrix adheres constraints. coding in r , familiar lpsolve , optimx packages, don't know how use them particular question.
because quw , ruw both data, constraints objective linear in yuw decision variables. result, linear programming problem can approached lpsolve package.
to abstract out bit, let's assume r=20 , c=10 describe dimensions of input matrices. there r*c decision variables , can assign them order y11, y21, ... yr1, y12, y22, ... yr2, ..., y1c, y2c, ..., yrc, reading down columns of matrix of variables.
the coefficient of each yuw variable in objective -quw; note quw*ruw terms in summation constants (aka not affected values select decision variables) , therefore not inputted linear programming solver. interestingly, means ruw has no effect on optimization model solution.
the first r*(c-1) constraints correspond yuw >= yu,w+1 constraints, , next (r-1)*c constraints correspond yuw >= yu-1,w constraints. final 2 constraints correspond y20,1 >= 100 , y1,10 >= 0 constraints.
we can input model lpsolve package following r command, taking input simple q matrix each entry -1 (the resulting solution should have decision variables set 0 except bottom left corner, should 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 access model solution:
mod # success: objective function 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 note model become infeasible if quw changed (for instance if filled 1 instead of -1). in these cases model exit status 3 (you can see running mod or mod$status).
Comments
Post a Comment