有约束时如何在 R 中进行优化

Question

我需要优化以下object 函数。 qgpd 来自包fExtremes。

var.asym <- function(alpha1, alpha2, xi, beta, n){
  term11 <- alpha1*(1-alpha1)^(2*xi-1)
  term12 <- alpha1*(1-alpha1)^(xi-1)*(1-alpha2)^xi
  term22 <- alpha2*(1-alpha2)^(2*xi-1)
  Sigma <- matrix(c(term11, term12, term12, term22), nrow=2, byrow=TRUE)
  Sigma*beta^2/n
}

mop.jacob.inv <- function(alpha1, alpha2, xi, beta){
  term11 <- -qgpd(alpha1, xi, beta)/xi - beta*(1-alpha1)^xi*log(1-alpha1)/xi
  term12 <- qgpd(alpha1, xi, beta)/beta
  term21 <- -qgpd(alpha2, xi, beta)/xi - beta*(1-alpha2)^xi*log(1-alpha2)/xi
  term22 <- qgpd(alpha2, xi, beta)/beta
  jacob <- matrix(c(term11, term12, term21, term22), nrow=2, byrow=TRUE)
  jacob.inv <- solve(jacob)
  jacob.inv
}

var.asym2 <- function(alpha1, alpha2) var.asym(alpha1, alpha2, 0.2, 1, 1000)

mop.jacob.inv2 <- function(alpha1, alpha2) mop.jacob.inv(alpha1, alpha2, 0.2, 1)

# Function to be optimised:

object <- function(alpha1, alpha2){
  term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))
  sum(diag(term1))
}

为了最小化 object，我有一个额外的约束 0 < alpha1 < alpha2 < 1。我的问题是我是否可以使用R 中的通用optim 函数来做到这一点。如果是这样，语法是什么，即如何在R 中设置问题？还有其他和/或更好的方法吗？如果没有，有没有办法在R 中做到这一点？谢谢。

更新：

在评论的帮助下，我有以下几点：

object <- function(alpha1, alpha2){
  term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))
  1/sum(diag(term1))*(alpha1>0)*(alpha2>alpha1)*(alpha2<1)
}

optim(c(0.01, 0.75), object)

然后我收到错误Error in stopifnot(min(p, na.rm = TRUE) >= 0) : argument "alpha2" is missing, with no default。出了什么问题？

Answer 1

您可以为此使用constrOptim(...)。我们需要稍微改变object的定义。

object <- function(alpha){
  alpha1 <- alpha[1]
  alpha2 <- alpha[2]
  term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))
  sum(diag(term1))
}
ui <- matrix(c(1,0,-1,0,-1,1),nc=2)
ci <- c(0,-1,0)
result <- constrOptim(th=c(0.4,0.6),object, grad=NULL, ui=ui, ci=ci)
result$par
# [1] 1.962097e-10 7.962686e-01

使用ui=... 和ci=... 参数应用约束。 ui 是 k x p 矩阵，其中 p 是参数的数量（在您的情况下为 2），k 是约束的数量（在您的情况下为 3），ci 是长度为 @ 987654331@。必须指定约束，以便：

ui × alpha - ci ≥ 0

所以在你的情况下，约束是：

α₁ ≥ 0

-α₂ + 1 ≥ 0

-α₁ + α₂ ≥ 0

以上代码中ui 和ci 的定义强制执行这些约束。

我们可以使用网格搜索来检查结果。

# check using grid search
x <- seq(0,1,by=0.1)
m <- expand.grid(a1=x,a2=x)
m <- m[m$a1<m$a2,]
grid <- apply(m,1,object)
m[which.min(grid),]
#    a1  a2
# 89  0 0.8

这会产生非常接近优化解决方案的结果。