使用梯度下降优化函数

Question

生长度日是植物物候学中的一个概念，特定作物每天需要积累一定量的热量单位才能从一个阶段移动到另一个阶段。

我有一个给定站点 10 年的每日分辨率的热单位数据，如下所示：

  set.seed(1)
  avg_temp <- data.frame(year_ref = rep(2001:2010, each = 365),
                         doy = rep(1:365, times = 10),
                         thermal.units = sample(0:40, 3650, replace=TRUE))

我在这个网站上也种了一种作物，如果在第 152 天种植，则需要 110 天才能成熟

  planting_date <- 152 
  observed_days_to_mature <- 110

我也有一些初步的随机猜测，关于这种作物在从种植开始到完全成熟的每个阶段通常可能积累多少热单位。例如在下面的例子中，第 1 阶段需要从种植开始积累 50 个热量单位，第 2 阶段需要 120 个热量单位，因为种植，第 3 阶段需要 190 个热量单位，因为种植等等。

  gdd_data <- data.frame(stage_id = 1:4,
                         gdd_required = c(50, 120, 190, 250))

所以根据 gdd 要求，我可以计算每年这种作物需要多少天才能成熟。

  library(dplyr)
  library(data.table)

  days_to_mature_func <- function(gdd_data_df, avg_temp_df, planting_date_d){
    
    gdd.vec <- gdd_data_df$gdd_required 
    
    year_vec <- sort(unique(avg_temp_df$year_ref))
    
    temp_ls <- list()
    
    for(y in seq_along(year_vec)){
      
      year_id <- year_vec[y]
      
      weather_sub <- avg_temp_df %>% 
                     dplyr::filter(year_ref == year_id & 
                                   doy >= planting_date_d)  
      
      stage_vec <- unlist(lapply(1:length(gdd.vec), function(x)  planting_date_d - 1 + which.max(cumsum(weather_sub$thermal.units) >= gdd.vec[x])))
      
      stage_vec[length(stage_vec)] <- ifelse(stage_vec[length(stage_vec)] <= stage_vec[length(stage_vec) - 1], NA, stage_vec[length(stage_vec)])
      
      gdd_doy <- as.data.frame(t(as.data.frame(stage_vec)))
      
      names(gdd_doy) <- paste0('stage_doy', 1:length(stage_vec))
      gdd_doy$year_ref <- year_id
      temp_ls[[y]] <- gdd_doy
  }
    days_to_mature_mod <- rbindlist(temp_ls)
    return(days_to_mature_mod)
  }

  days_to_mature_mod <- days_to_mature_func(gdd_data, avg_temp, planting_date)
  days_to_mature_mod

  stage_doy1 stage_doy2 stage_doy3 stage_doy4 year_ref
  1:        154        160        164        167     2001
  2:        154        157        159        163     2002
  3:        154        157        160        162     2003
  4:        155        157        163        165     2004
  5:        154        156        160        164     2005
  6:        154        161        164        168     2006
  7:        154        156        159        161     2007
  8:        155        158        161        164     2008
  9:        154        156        160        163     2009
  10:       154        158        160        163     2010

由于作物需要 110 天才能成熟，我将误差定义为：

  error_mod <- mean(days_to_mature_mod$stage_doy4 - observed_days_to_mature)^2

我的问题是如何优化 gdd_data 中的 gdd_required 以产生最小的错误。

我实施的一种方法是遍历一系列因素，减少 gdd_required 每一步并计算误差。误差最低的因素是我应用的最终因素 gdd_required 数据。我正在阅读有关梯度下降算法的文章，该算法可能会使此过程更快，但不幸的是，我还没有足够的技术专长来实现这一目标。

来自评论： 我确实有一个不明确的条件 - 我正在优化的函数中的 x 是有序的，即 x[1] < x[2] < x[3] < x[4] 因为它们是累积的.

Answer 1

根据您的示例，您可以定义一个函数，该函数采用任意 gdd_required 和 returns 拟合：

optim_function <- function(x){
  gdd_data <- data.frame(stage_id = 1:4, gdd_required = x)
  days_to_mature_mod <- days_to_mature_func(gdd_data, avg_temp, planting_date)
  error_mod <- mean(days_to_mature_mod$stage_doy4 - observed_days_to_mature)^2
}

函数 optim 允许您找到达到最小值的参数，从您使用的初始设置开始，例如

optim(c(50, 120, 190, 250), optim_function)
#$par
#[1] 266.35738 199.59795 -28.35870  30.21135
# 
#$value
#[1] 1866.24
# 
#$counts
#function gradient 
#      91       NA 
#
#$convergence
#[1] 0
#
#$message
#NULL

因此，使用参数 266.35738、199.59795、-28.35870、30.21135 找到了 1866 年左右的最佳拟合。

帮助页面提供了一些关于进行约束优化的指示，如果它们在特定范围内很重要的话。

鉴于您认为参数应该严格递增的评论，您可以使用 cumsum(exp()) 将任意值转换为递增的值，这样您的代码将变为

optim_function_plus <- function(x){
  gdd_data <- data.frame(stage_id = 1:4, gdd_required = cumsum(exp(x)))
  days_to_mature_mod <- days_to_mature_func(gdd_data, avg_temp, planting_date)
  error_mod <- mean(days_to_mature_mod$stage_doy4 - observed_days_to_mature)^2
}

opt <- optim(log(c(50, 70, 70, 60)), optim_function_plus)
opt
# $par
# [1] 1.578174 2.057647 2.392850 3.241456
# 
# $value
# [1] 1953.64
# 
# $counts
# function gradient 
# 57       NA 
# 
# $convergence
# [1] 0
# 
# $message
# NULL

要使参数恢复到您感兴趣的比例，您需要执行以下操作：

cumsum(exp(opt$par))
# [1]  4.846097 12.673626 23.618263 49.189184

使用梯度下降优化函数

Optimize a function using gradient descent

r

optimization

gradient-descent