R语言 数值优化但使用向量

9wbgstp7  于 2023-06-27  发布在  其他
关注(0)|答案(1)|浏览(113)

我是新的解决一些数字,所以我要求这得到一个入门的方法,以一个问题,我真的很清楚。
假设你有这个优化问题:

其中您知道\gamma_c\gamma_\ell\tau\Bar{\alpha}的值
我用拉格朗日乘子手工求解,得到了一个封闭形式的解。因此,我对消费(c)、休闲(\ell)和劳动力供给(h)给出了这些答案。

所以问题是我可以计算最优值(c\ellh)如下:(我在R中这样做了,但Python或Julia中的过程可能非常相似)

library(tidyverse)

w_par = c(4, 0.4)
i_par = c(3, 0.04)
e_par = c(0, 0.01^2)
gamma_l = 8; gamma_c = 50; tau = 0.08; Time = 24; alpha_bar = 0.7;N = 10000
gamma_h = Time - gamma_l
theta_true = c(gamma_h, gamma_c, alpha_bar, sqrt(e_par[2]))

set.seed(1)
df <- data.frame(w = exp(rnorm(n = N, mean = w_par[1], sd = sqrt(w_par[2]))),
                 I = exp(rnorm(n = N, mean = i_par[2], sd = sqrt(i_par[2]))),
                 e = rnorm(n = N, mean = e_par[1], sd = sqrt(e_par[2]))) %>% 
  mutate(a = alpha_bar + e,
         h = a*gamma_h - (((1-a)*(I-gamma_c))/((1-tau)*w)),
         L = Time - h,
         C = (1-tau)*w*h+I,
         U = a*log(C-gamma_c) + (1-a)*log(L-gamma_l))

好了,现在取 Dataframe 的第一部分,它只包含这4个变量(wIea)和参数。
有没有一种方法可以通过优化器获得hLC(optima)?我应该遵循哪些步骤来找到最佳列?使用优化器获得的列是否与我使用封闭式解决方案获得的列具有相同的值?
我不需要一个超级明确的答案,但一些开始弄清楚如何做到这一点。
我陈述这个小模型是因为我知道有一个封闭形式的解。但是为了工作,我必须得到一个没有封闭形式解的模型的最优解,他们告诉我的只是用数值求解(我不知道怎么做,但我愿意学习)
先谢谢你了!
编辑:在我的符号中有一个错字,而不是T_i只是T

gg0vcinb

gg0vcinb1#

不知道为什么会投反对票,但这里有一个实现,其中x是一个向量--只需要在函数的前几行分配元素。

U_func_3 <- function(x, w=NULL, I=NULL, e=NULL, a=NULL){
  
  h <- x[1]
  L <- 24-x[1]
  C <-(1-tau)*w*h+I
  

  ## the U to be maximized
  U <- a*log(C-gamma_c) + (1-a)*log(L-gamma_l)
  -U
  
}

  
use_apply_and_optim <-
    apply(df,
          1,
          function(DAT){
        -optim(c(Time/2, Time/2, 600), fn=U_func_3,
                
                w = DAT[1],
                I = DAT[2],
                e = DAT[3],
                a = DAT[4],
                
                
          method="L-BFGS-B",
          upper=c(Time, Time, Inf),
          lower=c(0,gamma_l, gamma_c ))$value}
    )

我能够将问题简化为一个具有动态下限的单变量优化问题。本例中的解决方案使用apply()逐行遍历数据集,然后使用optim()获取数据值以通知动态下限。如果我没看错模型和约束,

  • 对于给定的h,如果我们知道Time,则确定\ell,并且
  • 对于给定的h,如果我们知道tauwI,则确定C

自然地,自然对数中的量需要大于0,因此

  • 求解c > gamma_c使得h的动态下限和
  • 求解l > gamma_c使得h的上限为静态Time - gamma_l
## define U function
U_func <- function(x, w=NULL, I=NULL, e=NULL, a=NULL){
  
  h <- x[1]
  L <- Time - x[1]
  C <- (1-tau)*w*h+I
  
  ## the U to be maximized
  U <- a*log(C-gamma_c) + (1-a)*log(L-gamma_l)
  -U
  
}

## store results of optim() via row-wise apply
use_apply_and_optim <-
    apply(df,
          1,
          function(DAT){
        -optim(Time/2, fn=U_func,
                
                w = DAT[1],
                I = DAT[2],
                e = DAT[3],
                a = DAT[4],
                       
          method="L-BFGS-B",
          upper=c(Time - gamma_l),
          lower=c( (gamma_c - DAT[2])/((1-tau)*DAT[1]) ))$value}
    )

head(cbind(closed_solution = df$U, use_apply_and_optim))

结果如下:

> head(cbind(closed_solution = df$U, use_apply_and_optim))
     closed_solution use_apply_and_optim
[1,]        4.541093            4.541093
[2,]        4.940625            4.940625
[3,]        4.396769            4.396769
[4,]        5.476242            5.476242
[5,]        5.050314            5.050314
[6,]        4.419881            4.419881

直线:

plot(df$U, use_apply_and_optim)
abline(a=0,b=1,col="blue")

不完全相同,但MSE非常低,可能存在舍入误差(?):

> identical(df$U, use_apply_and_optim)
[1] FALSE
> mean((df$U-use_apply_and_optim)^2)
[1] 1.127322e-17

相关问题