替换重采样
Resampling with replacement
我正在尝试为我的模型获取(引导的)输入数据。
源文件:https://www.dropbox.com/s/dudzxhozr50uhr7/EddyData_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)
}
这会生成带有替换的重采样,例如
输出显示 1000 个 bootstraps 在一天的 48 个半小时内通过替换重新采样。
这是我的问题:
带替换的重采样在一天中随机混合半小时,没有任何控制。例如,我不想把晚上的半小时和白天的半小时混在一起。结果导致事后进行奇怪的计算。有没有可能以我定义允许什么和不允许什么的方式对此进行编码?
示例:
- 只允许在晚上 10 点到下午 5 点之间重新采样
- 夜间时间不能用白天时间重新采样(反之亦然)
我对 CRD 设计做了幼稚的 bootstrap,但对时间数据没有做过。那是时间序列数据吗?据我了解,您希望仅在下午 2 点而不是下午 3 点对下午 2 点进行采样。了解抽样对于进行正确的分析很重要,因为在 R 中很容易出错。
我注意到您对 bootstrap 使用了循环而不是包。
我使用 'boot' 包用于天真的 bootstrap,所以我转向 Google 查看时间数据。这是我发现的,我很抱歉这就是我所拥有的(由于缺乏代表我无法发表评论)
使用启动包我敢打赌任何东西都比使用循环更快
您可以使用 system.time( ) 检查,尤其是当您有大量数据时。
https://stat.ethz.ch/R-manual/R-devel/library/boot/html/tsboot.html
这是我如何处理我的天真 bootstrap:
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/
无论如何,抱歉,我没有提供太多帮助,但想分享一些想法。
我正在尝试为我的模型获取(引导的)输入数据。
源文件:https://www.dropbox.com/s/dudzxhozr50uhr7/EddyData_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)
}
这会生成带有替换的重采样,例如
输出显示 1000 个 bootstraps 在一天的 48 个半小时内通过替换重新采样。
这是我的问题:
带替换的重采样在一天中随机混合半小时,没有任何控制。例如,我不想把晚上的半小时和白天的半小时混在一起。结果导致事后进行奇怪的计算。有没有可能以我定义允许什么和不允许什么的方式对此进行编码?
示例:
- 只允许在晚上 10 点到下午 5 点之间重新采样
- 夜间时间不能用白天时间重新采样(反之亦然)
我对 CRD 设计做了幼稚的 bootstrap,但对时间数据没有做过。那是时间序列数据吗?据我了解,您希望仅在下午 2 点而不是下午 3 点对下午 2 点进行采样。了解抽样对于进行正确的分析很重要,因为在 R 中很容易出错。
我注意到您对 bootstrap 使用了循环而不是包。 我使用 'boot' 包用于天真的 bootstrap,所以我转向 Google 查看时间数据。这是我发现的,我很抱歉这就是我所拥有的(由于缺乏代表我无法发表评论) 使用启动包我敢打赌任何东西都比使用循环更快 您可以使用 system.time( ) 检查,尤其是当您有大量数据时。
https://stat.ethz.ch/R-manual/R-devel/library/boot/html/tsboot.html
这是我如何处理我的天真 bootstrap:
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/
无论如何,抱歉,我没有提供太多帮助,但想分享一些想法。