r - 有约束的优化

Question

我正在处理一个模型的输出，其中的参数估计可能不符合先验预期。我想编写一个函数来强制这些效用估计符合这些期望。为此，该函数应最小化起始值与新估计值之间的平方偏差之和。由于我们有先验期望，优化应该受到以下约束:

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

例如，下面的原始参数估计值针对 B2 和 B3 进行了翻转。栏目 Delta和 Delta^2显示原始参数估计与新系数之间的偏差。我正在尝试最小化列 Delta^2 .我已经在 Excel 中对此进行了编码，并展示了 Excel 的求解器如何通过提供一组约束来优化此问题:

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 中设置它。我确定我只是有点密集，但可以在正确的方向使用一些指针!

2012 年 3 月 24 日 - 添加了赏金，因为我不够聪明，无法翻译第一个答案。

这是一些应该在正确路径上的 R 代码。假设 Beta 开始于:

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]

IE。，

-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]).

有几个例子
here
或 there .