r - 生态时间序列模型中参数蛮力估计的替代方法

Question

我正在模拟一个水文过程（湖泊中的水位 [阶段] 以毫米为单位测量），可以描述为：

其中是从不同的模型估计的，并在该模型中用作常数。是未知数，该值预计在 (-0.001,0.001) 之间。S的起始值只要大于10m（10000mm）即可。该模型以每日时间步长运行。我从多个不同的湖泊观察了舞台，并独立地适应了每个湖泊。

目前，我通过以下方式暴力识别参数值：

1. 创建跨越 (-0.001,0.001) 的参数值的 100 值序列
2. 使用上述方程进行预测阶段并估计建模数据和观测数据之间的 RMSE（观测值明显少于建模数据点）
3. 识别具有最低 RMSE的B并在任一侧选择B值以创建新的参数值序列以进行搜索
4. 重复步骤 2 和 3，直到 RMSE 减少小于 0.01 或增加。

下面是我一直在使用的蛮力方法的代码以及与单个湖相关的数据。

考虑到上述模型以及我只观察到有限天数的数据这一事实，是否有另一种方法来估计未知参数Beta2 ？

谢谢！

library(tidyverse)
library(lubridate)
library(Metrics)

#The Data
dat <- read_csv("https://www.dropbox.com/s/skg8wfpu9274npb/driver_data.csv?dl=1")
observeddata <-    read_csv("https://www.dropbox.com/s/bhh27g5rupoqps3/observeddata.csv?dl=1") %>% select(Date,Value)

#Setup initial values and vectors
S = matrix(nrow = nrow(dat),ncol = 1) #create an empty matrix for predicted values
S[1,1] = 10*1000 #set initial value (mm)
rmse.diff <- 10^100 #random high value for difference between min RMSE between successive
                    #parameter searches
b.levels <- seq(from = -0.001,to = 0.001,length.out=100) #random starting parameter that should contain
                                                         #the final value being estimated 
n = 0 # counter

#Loop to bruteforce search for best parameter estimate
while(rmse.diff > 0.01 ) {
  rmse.vec = rep(NA,length(b.levels))
  for(t in 1:length(b.levels)){
    for(z in 2:nrow(S)){
      S[z,1] <- S[(z-1),1] + (1.071663*(dat$X[z])) + (b.levels[t]*(S[(z-1),1])) #-1.532236
    } #end of time series loop
    extrap_level <- data.frame(Date= dat$Date, level = S) # predicted lake levels

    #calculate an offset to center observed data on extrapolated data 
    dat.offset = observeddata %>% left_join(extrap_level) %>% 
      mutate(offset = level-Value) %>% drop_na()
    offset <- mean(dat.offset$offset)
    dat.compare <- observeddata %>% left_join(extrap_level) %>% 
      mutate(Value = Value + offset) %>% drop_na()
    #calculate RMSE between observed and extrapolated values
    rmse.vec[t] <- rmse(actual = dat.compare$Value,predicted = dat.compare$level)
    #plot the data to watch how parameter choice influences fit while looping
    #plots have a hard time keeping up
    if(t ==1 | t==50 | t==100) {
    plot(extrap_level$Date,extrap_level$level,type="l")
    lines(dat.compare$Date,dat.compare$Value,col="red")
    }
  }
  #find minimum RMSE value
  min.rmse <- which(rmse.vec==min(rmse.vec))
  if(n == 0) rmse.best <- rmse.vec[min.rmse] else rmse.best = c(rmse.best,rmse.vec[min.rmse])
  if (n >= 1) rmse.diff <- (rmse.best[n]-rmse.best[n+1])
  if(rmse.diff < 0) break()
  best.b <- b.levels[min.rmse]
  #take the parameter values on either side of the best prior RMSE and use those as search area
  b.levels <- seq(from = (b.levels[min.rmse-1]),to = (b.levels[min.rmse+1]),length.out=100)
  n = n + 1
}

rmse.best #vector of RMSE for each parameter search 
best.b #Last identified best parameter value

标签： rtime-series