如何将方程内的这个最大化问题转化为 R?
How do I translate this maximization problem inside an equation to R?
各位程序员。我正在研究一本关于经济学数值解的书(Judd 1998)。我正在尝试在 R 中重现同一本书中的问题,因此我可以使用 optim
包来查看是否可以获得类似的结果。
作者建立的问题是这个one: and his results were these.
我试图将这个问题转录到 R,结果是这个代码块:
DisutilityJudd <- function(L){
if(L == 0){
return(0)
}else{
return(0.1)
}
}
AgentUtilityJudd <- function(w, L){
(-exp(-2*w) + 1) - DisutilityJudd(L)
}
reservation.utility.judd <- AgentUtilityJudd(1, 1)
MaxEffortUtility <- function(w1, w2, L = 1){
0.8 * AgentUtilityJudd(w1, L) + 0.2 * AgentUtilityJudd(w2, L)
}
LeastEffortUtility <- function(w1, w2, L = 0){
0.4 * AgentUtilityJudd(w1, L) + 0.6 * AgentUtilityJudd(w2, L)
}
UtilityDifferenceJudd <- function(w1, w2){
MaxEffortUtility(w1, w2) - LeastEffortUtility(w1, w2)
}
PenaltyFunctionJudd <- function(w1, w2, P = 100000){
if(length(w1) == 2){
y <- -1 * (0.8 * (2 - w1[1]) - 0.2 * w1[2] - P *
(pmax(0, -MaxEffortUtility(w1[1], w1[1]) - reservation.utility.judd))^2 -
P * (pmax(0, -UtilityDifferenceJudd(w1[1], w1[1])))^2)
}else{
y <- -1 * (0.8 * (2 - w1) - 0.2 * w2 - P *
(pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd))^2 -
P * (pmax(0, -UtilityDifferenceJudd(w1, w2)))^2)
}
return(y)
}
没有错误,但我的代码生成的结果与我的预期相去甚远:
optim(c(1.1, 0.5), PenaltyFunctionJudd)
$par
[1] 1.343909e+49 -2.370681e+51
$value
[1] -4.633849e+50
$counts
function gradient
501 NA
$convergence
[1] 1
$message
NULL
可能是我的惩罚函数有问题。我假设这是由于 pmax
函数造成的。有人可以帮我鉴定一下吗?谢谢,感谢您的关注。
编辑:打字错误。
我相信你的意思是 w1[2]
在 if(length(w1) == 2)
为真时。
我已经修改了你的代码,没有触及你如何定义以前的功能。目前尚不清楚它是否是预期的结果:IV(-1) 是什么意思,结果是负 1 吗? 10 的幂 ?
PenaltyFunctionJudd <- function(w1, w2, P = 1e5){
if(length(w1) > 1){
w2 <- w1[2]
w1 <- w1[1]
}
# cat("length is 2 \n")
y <- 0.8 * (2 - w1) - 0.2 * w2 - P *
( pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd) )^2 -
P * ( pmax(0, -UtilityDifferenceJudd(w1, w2)) )^2
# cat("pmax1 :", pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd), "\n")
# cat("pmax2 :", pmax(0, -UtilityDifferenceJudd(w1, w2)), "\n")
return(y)
}
optim(c(1.1, 0.5), PenaltyFunctionJudd, control = list(fnscale = -1) )
optim(c(11, 5), PenaltyFunctionJudd, method = "BFGS", control = list(fnscale = -1, maxit = 100) )
您可以使用 cat
或 print
来检查您的值(这里我注意到一些 Inf 和 0 导致我注意到代码错误)。
友情提示:如果你正确定义了之前的功能,优化中存在很多不稳定因素(问题设置不当?需要更多惩罚?)。的确当运行两倍以上算法参数波动很大...
各位程序员。我正在研究一本关于经济学数值解的书(Judd 1998)。我正在尝试在 R 中重现同一本书中的问题,因此我可以使用 optim
包来查看是否可以获得类似的结果。
作者建立的问题是这个one: and his results were these.
我试图将这个问题转录到 R,结果是这个代码块:
DisutilityJudd <- function(L){
if(L == 0){
return(0)
}else{
return(0.1)
}
}
AgentUtilityJudd <- function(w, L){
(-exp(-2*w) + 1) - DisutilityJudd(L)
}
reservation.utility.judd <- AgentUtilityJudd(1, 1)
MaxEffortUtility <- function(w1, w2, L = 1){
0.8 * AgentUtilityJudd(w1, L) + 0.2 * AgentUtilityJudd(w2, L)
}
LeastEffortUtility <- function(w1, w2, L = 0){
0.4 * AgentUtilityJudd(w1, L) + 0.6 * AgentUtilityJudd(w2, L)
}
UtilityDifferenceJudd <- function(w1, w2){
MaxEffortUtility(w1, w2) - LeastEffortUtility(w1, w2)
}
PenaltyFunctionJudd <- function(w1, w2, P = 100000){
if(length(w1) == 2){
y <- -1 * (0.8 * (2 - w1[1]) - 0.2 * w1[2] - P *
(pmax(0, -MaxEffortUtility(w1[1], w1[1]) - reservation.utility.judd))^2 -
P * (pmax(0, -UtilityDifferenceJudd(w1[1], w1[1])))^2)
}else{
y <- -1 * (0.8 * (2 - w1) - 0.2 * w2 - P *
(pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd))^2 -
P * (pmax(0, -UtilityDifferenceJudd(w1, w2)))^2)
}
return(y)
}
没有错误,但我的代码生成的结果与我的预期相去甚远:
optim(c(1.1, 0.5), PenaltyFunctionJudd)
$par
[1] 1.343909e+49 -2.370681e+51
$value
[1] -4.633849e+50
$counts
function gradient
501 NA
$convergence
[1] 1
$message
NULL
可能是我的惩罚函数有问题。我假设这是由于 pmax
函数造成的。有人可以帮我鉴定一下吗?谢谢,感谢您的关注。
编辑:打字错误。
我相信你的意思是 w1[2]
在 if(length(w1) == 2)
为真时。
我已经修改了你的代码,没有触及你如何定义以前的功能。目前尚不清楚它是否是预期的结果:IV(-1) 是什么意思,结果是负 1 吗? 10 的幂 ?
PenaltyFunctionJudd <- function(w1, w2, P = 1e5){
if(length(w1) > 1){
w2 <- w1[2]
w1 <- w1[1]
}
# cat("length is 2 \n")
y <- 0.8 * (2 - w1) - 0.2 * w2 - P *
( pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd) )^2 -
P * ( pmax(0, -UtilityDifferenceJudd(w1, w2)) )^2
# cat("pmax1 :", pmax(0, -MaxEffortUtility(w1, w2) - reservation.utility.judd), "\n")
# cat("pmax2 :", pmax(0, -UtilityDifferenceJudd(w1, w2)), "\n")
return(y)
}
optim(c(1.1, 0.5), PenaltyFunctionJudd, control = list(fnscale = -1) )
optim(c(11, 5), PenaltyFunctionJudd, method = "BFGS", control = list(fnscale = -1, maxit = 100) )
您可以使用 cat
或 print
来检查您的值(这里我注意到一些 Inf 和 0 导致我注意到代码错误)。
友情提示:如果你正确定义了之前的功能,优化中存在很多不稳定因素(问题设置不当?需要更多惩罚?)。的确当运行两倍以上算法参数波动很大...