如何在约束优化中将参数总和设置为1

6
这里是代码(如果太长,我很抱歉,但这是我手头唯一的例子);我正在使用A. Wittmann的CreditMetrics包中的CVaR示例以及DEoptim求解器进行优化:
library(CreditMetrics)
library(DEoptim)

N <- 3
n <- 100000
r <- 0.003
ead <- rep(1/N,N)
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
lgd <- 0.99
rating <- c("BBB", "AA", "B")   
firmnames <- c("firm 1", "firm 2", "firm 3")
alpha <- 0.99

# correlation matrix
rho <- matrix(c(  1, 0.4, 0.6,
                  0.4,   1, 0.5,
                  0.6, 0.5,   1), 3, 3, dimnames = list(firmnames, firmnames),
              byrow = TRUE)

# one year empirical migration matrix from standard&poors website
rc <- c("AAA", "AA", "A", "BBB", "BB", "B", "CCC", "D")
M <- matrix(c(90.81,  8.33,  0.68,  0.06,  0.08,  0.02,  0.01,   0.01,
              0.70, 90.65,  7.79,  0.64,  0.06,  0.13,  0.02,   0.01,
              0.09,  2.27, 91.05,  5.52,  0.74,  0.26,  0.01,   0.06,
              0.02,  0.33,  5.95, 85.93,  5.30,  1.17,  1.12,   0.18,
              0.03,  0.14,  0.67,  7.73, 80.53,  8.84,  1.00,   1.06,
              0.01,  0.11,  0.24,  0.43,  6.48, 83.46,  4.07,   5.20,
              0.21,     0,  0.22,  1.30,  2.38, 11.24, 64.86,  19.79,
              0,     0,     0,     0,     0,     0,     0, 100
)/100, 8, 8, dimnames = list(rc, rc), byrow = TRUE)

cm.CVaR(M, lgd, ead, N, n, r, rho, alpha, rating)

y <- cm.cs(M, lgd)[which(names(cm.cs(M, lgd)) == rating)]

现在我写下了我的函数...
fun <- function(w) {
  # ... 
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, 
                           rho, alpha, rating)
}

我想对它进行优化:

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
        control = DEoptim.control())

请告诉我,在优化的过程中,我需要在 # ... 中插入什么内容才能使得 sum(w) = 1

以下是根据flodel的提示所展示的优化结果:

# The first trick is to include B as large number to force the algorithm to put sum(w) = 1

fun <- function(w) {
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) + 
    abs(10000 * (sum(w) - 1))
}

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
        control = DEoptim.control())

$optim$bestval
[1] -0.05326055

$optim$bestmem
par1        par2        par3 
0.005046258 0.000201286 0.994752456

parsB <- c(0.005046258, 0.000201286, 0.994752456)

> fun(parsB)
            [,1]
[1,] -0.05326089

如您所见,第一种技巧的效果更好,因为他找到了比第二种结果更小的结果。不幸的是,似乎他需要更长的时间。

# The second trick needs you use w <- w / sum(w) in the function itself

fun <- function(w) {
  w <- w / sum(w)
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) #+ 
    #abs(10000 * (sum(w) - 1))
}

DEoptim(fn = fun, lower = rep(0, N), upper = rep(1, N), 
        control = DEoptim.control())

$optim$bestval
[1] -0.0532794

$optim$bestmem
par1         par2         par3 
1.306302e-15 2.586823e-15 9.307001e-01

parsC <- c(1.306302e-15, 2.586823e-15, 9.307001e-01)
parC <- parsC / sum(parsC)

> fun(parC)
           [,1]
[1,] -0.0532794

有什么意见吗?

如果要优化的函数“过于随机”,我是否应该增加迭代次数?


在上面的“第二个技巧”的代码中,您忘记在“fun”的主体中添加w <- w / sum(w)。您能否更新您的代码和结果? - flodel
谢谢,我刚刚更新了。根据这些结果,最好的方法是什么? - Lisa Ann
你说得对,我删错了。我忘记了将变量的上限设为1,如果在函数体中不插入 w <- w / sum(w),将会得到错误的结果。 - Lisa Ann
关于 w <- c(w, 1-sum(w)),它在 t(w) %*% y 中返回了一个数组不兼容的参数错误。 - Lisa Ann
根据这种最后的方法,最后一个变量(通过1-sum(w)获得的变量)可能是负数吗?这将违反我的限制。 - Lisa Ann
显示剩余2条评论
3个回答

9

尝试:

w <- w / sum(w)

如果DEoptim给出一个最优解w*,使得sum(w*) != 1,那么w*/sum(w*)应该是您的最优解。

另一种方法是在除一个变量之外的所有变量上求解。我们知道最后一个变量的值必须为1 - sum(w),因此在函数体中写入:

w <- c(w, 1-sum(w))

并且对 DEoptim 返回的最优解做相同处理:w* <- c(w*, 1-sum(w*)) 这两种解决方案都要求您将问题重新表述为无约束优化(不考虑变量限制),以便可以使用 DEoptim。这强制您在 DEoptim 之外进行一些额外的工作,以恢复原始问题的解决方案。
回答您的评论,如果您想让 DEoptim 立即给出正确答案(即不需要后续转换),您也可以尝试在目标函数中包含惩罚成本:例如添加 B * abs(sum(w)-1),其中 B 是一些任意大的数字,因此 sum(w) 将被强制为 1

谢谢,flodel,那就是我习惯使用的解决方案。但我希望优化输出已经总和为1。几个月前,我记得有一个函数示例,它在函数本身中嵌入了这个约束条件,但我无法找回它 :( - Lisa Ann
感谢您的建议,flodel。不幸的是,我不擅长约束优化,所以我将向您展示我使用您两个建议得到的结果。 “大B技巧”在收敛速度上较慢:实际上,它试图最小化B,然后进入了我的方程的第一部分;参数归一化可以加快收敛速度,但是,在归一化之后,我得到了次优解。我将编辑我的问题,以便向您展示我的代码... - Lisa Ann

1
我认为你应该对任何偏差都加以惩罚。 在最小化问题中添加一个项+(sum(weights) - 1)^2 * 1e10。你会发现这个巨大的惩罚将强制权重总和为1!

0

通过你所应用的技巧:

fun <- function(w) {
  w <- w / sum(w)
  - (t(w) %*% y - r) / cm.CVaR(M, lgd, ead = w, N, n, r, rho, alpha, rating) #+ 
    #abs(10000 * (sum(w) - 1))
}

为什么你不在这种情况下使用优化?我认为它会更快。

这并不是一个真正的答案。请在评论部分提出问题和澄清请求,而不是在回答部分。 - Roger-123

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接