我正面临以下挑战,但似乎没有找到解决方案:
[基本上,我想为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错误从问题表达式转移到约束表达式。我真的不知道如何解决这个问题。有人有合适的解决方案吗?预先感谢!
这里是剩下的测试代码
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))