约束优化

Question

我正在使用模型的输出,其中存在可能不符合先验预期的参数估计.我想写一个函数,迫使这些效用估计值符合这些预期.为此,该函数应最小化起始值和新估计之间的平方偏差之和.由于我们有先验预测,因此优化应遵循以下约束:

B0 < B1
B1 < B2
...
Bj < Bj+1 

例如,下面的原始参数估计是针对B2和B3的翻转翻转.列Delta和Delta^2显示原始参数估计和新系数之间的偏差.我想尽量减少列Delta^2.我在Excel中对此进行了编码,并显示了Excel的Solver如何在提供约束条件的情况下优化此问题:

Beta    BetaRaw    Delta    Delta^2    BetaNew
B0       1.2       0        0          1.2
B1       1.3       0        0          1.3
B2       1.6       -0.2     0.04       1.4
B3       1.4       0        0          1.4
B4       2.2       0        0          2.2

通过阅读后?optim和?constrOptim,我不能神交如何设置这在R.我敢肯定,我只是作为一个有点迟钝,但可以在正确的方向上使用一些指针!

3/24/2012 - 添加了赏金,因为我不够聪明,无法翻译第一个答案.

这里有一些R代码应该在正确的路径上.假设测试版开始于:

betas <- c(1.2,1.3,1.6,1.4,2.2)

我想尽量减少以下功能 b0 <= b1 <= b2 <= b3 <= b4

f <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  x3 <- x[3]
  x4 <- x[4]
  x5 <- x[5]

 loss <- (x1 - betas[1]) ^ 2 + 
         (x2 - betas[2]) ^ 2 + 
         (x3 - betas[3]) ^ 2 + 
         (x4 - betas[4]) ^ 2 +
         (x5 - betas[5]) ^ 2    

  return(loss)
}

为了证明该函数有效,如果我们将原始beta版传递给:

> f(betas)
[1] 0

并且随机输入相对较大:

> set.seed(42)
> f(rnorm(5))
[1] 8.849329

并在我能够在Excel中计算的值最小化:

> f(c(1.2,1.3,1.4,1.4,2.2))
[1] 0.04

Answer 1

1. 由于目标是二次方且约束是线性的,因此可以使用solve.QP.

它发现b最小化

(1/2) * t(b) %*% Dmat %*% b - t(dvec) %*% b 

在约束下

t(Amat) %*% b >= bvec. 

在这里,我们希望b最小化

sum( (b-betas)^2 ) = sum(b^2) - 2 * sum(b*betas) + sum(beta^2)
                   = t(b) %*% t(b) - 2 * t(b) %*% betas + sum(beta^2).

从上一个词开始,sum(beta^2)是不变的,我们可以放弃它,我们可以设置

Dmat = diag(n)
dvec = betas.

约束是

b[1] <= b[2]
b[2] <= b[3]
...
b[n-1] <= b[n]

即

-b[1] + b[2]                       >= 0
      - b[2] + b[3]                >= 0
               ...
                   - b[n-1] + b[n] >= 0

这t(Amat)就是

[ -1  1                ]
[    -1  1             ]
[       -1  1          ]
[             ...      ]
[                -1  1 ]

并且bvec为零.

这导致以下代码.

# Sample data
betas <- c(1.2, 1.3, 1.6, 1.4, 2.2)

# Optimization
n <- length(betas)
Dmat <- diag(n)
dvec <- betas
Amat <- matrix(0,nr=n,nc=n-1)
Amat[cbind(1:(n-1), 1:(n-1))] <- -1
Amat[cbind(2:n,     1:(n-1))] <-  1
t(Amat)  # Check that it looks as it should
bvec <- rep(0,n-1)
library(quadprog)
r <- solve.QP(Dmat, dvec, Amat, bvec)

# Check the result, graphically
plot(betas)
points(r$solution, pch=16)

2. 您可以constrOptim以相同的方式使用(目标函数可以是任意的,但约束必须是线性的).

3. 更一般地说,例如,optim如果将问题重新参数化为非约束优化问题,则可以使用

b[1] = exp(x[1])
b[2] = b[1] + exp(x[2])
...
b[n] = b[n-1] + exp(x[n-1]).

这里 或那里有一些例子 .