置换重采样

Question

我正在尝试获取我的模型的(引导)输入数据。

源文件：2010.csv?dl=0

library("dplyr")
library("readr")
library("reshape2")
library("ggplot2")

sub <- read_csv("EddyData_2010.csv", 
                col_types = list(col_integer(), col_integer(), col_double(),
                                 col_double(), col_double(), col_double(),
                                 col_double(), col_double(), col_double(),
                                 col_double(), col_double(), col_double()),
                col_names = c("Year", "DoY", "Hour", "NEE", "LE", "H", "Rg",
                              "Tair", "Tsoil", "rH", "Ustar", "VPD")) %>%
  filter(DoY == 170) %>%
  mutate(hour = 1:48) %>%
  select(NEE:hour)

# Number of resampling 
n_resempling <- 1000 
result_resampling <- NULL

# Do the resampling
for (i in 1:n_resempling) {
  result_resampling <- sample(48, 48, replace = T) %>%
    slice(sub, .) %>%
    mutate(resempling_number = i) %>%
    bind_rows(. , result_resampling)
}

这将生成替换的重采样，如

​

​

输出结果显示，每天有48个半小时的替换，重新放大了1000个引导。

我的问题是：

重新采样和替换是随机混合半个小时的一天，没有任何控制。例如，我不想把半个小时的晚上和白天的半小时混为一谈。结果导致后来的奇怪计算。是否有可能以这样的方式对此进行编码:我定义了什么是允许的，什么是不允许的？

示例：

• 只允许从晚上10点到下午5点重新采样。
• 夜间时间不能用白天时间重新安排(反之亦然)

Answer 1

我为CRD设计做了天真的引导，但没有为时间数据。那是时间序列数据吗？据我所知，你想要的是下午2点而不是下午3点。理解采样对于做正确的分析是很重要的，因为在R中很容易出错。

我注意到您使用了一个循环来引导，而不是一个包。我用'boot‘软件包来进行简单的引导，所以我求助于Google来查看时间数据。下面是我发现的，我很抱歉，这是我所拥有的(由于缺乏rep而不能发表评论)使用引导包--我敢打赌，任何事情都比使用循环更快，您可以使用system.time( )来检查，特别是如果您有大量数据的话。

https://stat.ethz.ch/R-manual/R-devel/library/boot/html/tsboot.html

下面是我如何处理我天真的引导：

my.boot.fnx<-function(var, ind,alpha=0.95){
  newdf      <-na.omit(var[ind])
  mymean     <-mean(newdf,na.rm=TRUE)
  mytrimmean <-mean(newdf, trim=0.1, na.rm=TRUE)
  mymedian   <-median(newdf,na.rm=TRUE)
  mysd       <-sd(newdf,na.rm=TRUE)
  nonmiss    <- length(na.omit(newdf))
  semean     <- mysd/sqrt(nonmiss)
  lcl        <- mymean - qt(1-alpha/2,nonmiss-1)*semean
  ucl        <- mymean + qt(1-alpha/2,nonmiss-1)*semean
  mygini     <-
   sum(abs(outer(newdf,newdf,FUN="")))/
   length(newdf)/(length(newdf)-1)*sqrt(pi)/2
   c(mean=mymean,median=mymedian,se=semean, 
   lcl=lcl,ucl=ucl,sd=mysd,gsd=mygini)
#gini coef is a robust measure of SE
}
strap.df <- boot(df$var,statistic=my.boot.fnx, R=1000)

我还找到了时间数据http://eranraviv.com/bootstrapping-time-series-r-code/的源代码。

为不同设计提供适当的分析方法也是很好的资源：

http://people.stat.sfu.ca/~cschwarz/CourseNotes/

不管怎样，对不起，我没有得到太多的帮助，但我想分享一些想法。