0

我面临以下挑战，我似乎没有找到解决方案：

基本上，我想为 StoNED 编写一些代码 - 随机非平滑数据包络。

我从这个开始：

# dim(G) is 35x15
# dim(X) is 5x15
# dim(y) is 5x1
# dim(h) is 15x1

# y is a vector of output, x a matrix of input, 
h a vector of zeros, G a matrix of constraint parameters, betaa shall be estimated

library(CVXR)
betaa <- Variable(15,1)
target <- cvxr_norm(log(y)-log(X%*%betaa))
obj <- Minimize(target)
constr <- list(G%*%betaa>=h)
prob <- Problem(obj,constr)
result <- solve(prob)
print(result$getValue(betaa))

由于这不起作用，我尝试将问题转换为：

library(CVXR)
library(MASS)

#betaa <- Variable((15,1)
test <- Variable(5,1)
target <- cvxr_norm(log(y)-test)
obj <- Minimize(target)
constr <- list(G%*%(ginv(X)%*%exp(test))>=h)
prob <- Problem(obj,constr)
result <- solve(prob)
print(result$getValue(ginv(X)%*%exp(test)))

通过这个，我将 DCP 错误从问题表达式转移到了约束表达式。我真的不知道如何解决这个问题。有没有人有合适的解决方案？提前致谢！

4

1 回答 1

-1

这是用于测试的其余代码

input = 2
xmin = 5
xmax = 15
n=5

X <- cbind(runif(n, xmin, xmax))
if (input>=2){
  for (i in 2:input){
    Xi <-runif(n, xmin, xmax)
    X <- cbind(X,Xi)
  }
}

exponent=1
Y <- cbind(X[,1]^(exponent/ncol(X)))
if (input >=2){
  for (i in 2:input){
    Y <- Y*X[,i]^(exponent/ncol(X))
  }
}

u <- exp(-abs(rnorm(n, 0, 1.0877)))    
v <- exp(rnorm(n, 0, 2*0.1361))

Y_mod <- Y*u*v

COST = 0
MULT = 1
inter = 1
x=X
y=Y_mod
RTS = "VRS"

n <- nrow(as.matrix(x))  # No DMUs
m <- ncol(as.matrix(x))  # No Inputs

x <- as.matrix(x)
estimator <- "MM"

ifelse (RTS=="CRS",inter <- 0,inter <- 1) #
if (is.null(estimator)==TRUE){estimator <- "MM"; print("Moment estimator used
if not specified")}
if (RTS=="CRS") MULT <- 1 
if (RTS=="DRS") MULT <- 1
if (RTS=="IRS") MULT <- 1
if (RTS=="CRS" | RTS=="IRS" | RTS=="DRS") paste("Multiplicative model induced
for this scale assumption") 

if (MULT==1){ 
  if (inter==1){
    X <- diag(n);
    for (i in (1:m)){
      X<- cbind(X,diag(n)*as.matrix(x)[,i])
    }
  }

  if (inter==0){
    X <- diag(n)*as.matrix(x)[,1]
    if (m>1){
      for (i in (2:m)){X<- cbind(X,diag(n)*as.matrix(x)[,i])}}}
  B <- log(y)
}

if (MULT==0){
  if (inter==1){
    X <- diag(n);
    for (i in (1:m)){X<- cbind(X,diag(n)*as.matrix(x)[,i])}}

  if (inter==0){
    X <- diag(n)*as.matrix(x)[,1]
    if (m>1){
      for (i in (2:m)){X<- cbind(X,diag(n)*as.matrix(x)[,i])}}}
  B <- y
} 


###############################################
###### Create non-neg constraints
###############################################


non_neg_beta <- cbind(matrix(0,nrow=m*n,ncol=n),diag(m*n)) 
if(COST==0){ 
  if (inter==1){ 
    if (RTS=="VRS") {sign_cons <- non_neg_beta} 
    if (RTS=="DRS"){RTS_cons <- cbind(1*diag(n),matrix(0,nrow=n,ncol=m*n)); 
    sign_cons <- rbind(non_neg_beta, RTS_cons)}
    if (RTS=="IRS"){RTS_cons <- cbind(-1*diag(n),matrix(0,nrow=n,ncol=m*n)); 
    sign_cons <- rbind(non_neg_beta, RTS_cons)}
  }
}    
if(COST==1){ 
  if (inter==1){
    if (RTS=="VRS") {sign_cons <- non_neg_beta}
    if (RTS=="DRS"){RTS_cons <- cbind(-1*diag(n),matrix(0,nrow=n,ncol=n*m));
    sign_cons <- rbind(non_neg_beta, RTS_cons)}
    if (RTS=="IRS"){RTS_cons <- cbind(1*diag(n),matrix(0,nrow=n,ncol=n*m));
    sign_cons <- rbind(non_neg_beta, RTS_cons)}
  }
}  

if (inter==0){
  sign_cons <- diag(m*n)} 

###############################################
###### Create concavity constraints
###############################################
alpha_matlist <- list()
beta_matlist <- list()

for (i in 1:n){    
  a_matr <- diag(n);
  a_matr[,i] <- a_matr[,i]-1
  alpha_matlist[[i]] <- a_matr 

  # Create beta
  b_matr <-matrix(0,ncol=m*n,nrow=n);

  for (j in 1:m){
    b_matr[,i+n*(j-1)] <- -x[i,j]

    for (k in 1:n){
      b_matr[k,k+(j-1)*n] <- b_matr[k,k+(j-1)*n]+x[i,j]}
  }

  beta_matlist[[i]] <- b_matr}

if (inter==1){sspot_cons <-
cbind(do.call(rbind,alpha_matlist),do.call(rbind,beta_matlist))}
if (inter==0){sspot_cons <- cbind(do.call(rbind,beta_matlist))}  
#####################################
# !!!!!!!!! Cost function !!!!!!!!! 
#####################################
if (COST==1){sspot_cons <- -1*sspot_cons}

#####################################
### Getting the constraints together
#####################################
G <- rbind(sign_cons,sspot_cons) 
h <- rep(0,nrow(sign_cons)+nrow(sspot_cons))
于 2020-05-31T08:01:39.490 回答