标签: statistics-bootstrap

我试图使用boot.ciR的boot包来计算参数自举的偏差和偏斜校正自举置信区间.从我阅读的手册页和实验中,我得出结论,我必须自己计算折刀估计值并将其输入boot.ci,但这并未在任何地方明确说明.我一直无法找到其他文档,虽然公平地说我还没有查看代码所基于的原始Davison和Hinkley书籍......

如果我天真地跑b1 <- boot(...,sim="parametric"),然后boot.ci(b1),我得到错误influence values cannot be found from a parametric bootstrap.当且仅当我指定type="all"或时,才会发生此错误type="bca"; boot.ci(b1,type="bca")给出了同样的错误.那样做empinf(b1).我能得到的东西的工作的唯一办法是明确地计算刀切估计(使用empinf()与data参数)和饲料到这些boot.ci.

构建数据:

set.seed(101)
d <- data.frame(x=1:20,y=runif(20))
m1 <- lm(y~x,data=d)

引导:

b1 <- boot(d$y,
           statistic=function(yb,...) {
             coef(update(m1,data=transform(d,y=yb)))
           },
           R=1000,
           ran.gen=function(d,m) {
             unlist(simulate(m))
           },
           mle=m1,
           sim="parametric")

好到目前为止.

boot.ci(b1)
boot.ci(b1,type="bca")
empinf(b1)

都给出了上述错误.

这有效:

L <- empinf(data=d$y,type="jack",
            stype="i",
            statistic=function(y,f) {
              coef(update(m1,data=d[f,]))
            })

boot.ci(b1,type="bca",L=L)

有谁知道这是我应该这样做的方式吗？

更新:boot …

这是我第一次发帖到这个论坛,我想从一开始就说我不是一个熟练的程序员.如果问题或代码不清楚,请告诉我!

我试图通过引导来获得交互的95%置信区间(CI)(这是我的测试统计).我正在使用"boot"包.我的问题是,对于每个重新采样,我希望在受试者内进行随机化,以便来自不同受试者的观察结果不会混合.这是生成类似于我的数据帧的代码.如您所见,我有两个主体内因素("Num"和"Gram",我对两者之间的相互作用感兴趣):

Subject = rep(c("S1","S2","S3","S4"),4)
Num     = rep(c("singular","plural"),8) 
Gram    = rep(c("gram","gram","ungram","ungram"),4)
RT      = c(657,775,678,895,887,235,645,916,930,768,890,1016,590,978,450,920)
data    = data.frame(Subject,Num,Gram,RT) 

这是我用来获得经验交互值的代码:

summary(lm(RT ~ Num*Gram, data=data))

如您所见,我的两个因素之间的相互作用是-348.我想得到这个统计信息的bootstrap置信区间,我可以使用"boot"包生成:

# You need the following packages
install.packages("car") 
install.packages("MASS")
install.packages("boot")
library("car")
library("MASS")
library("boot")

#Function to create the statistic to be boostrapped
boot.huber <- function(data, indices) {
data <- data[indices, ] #select obs. in bootstrap sample
mod <- lm(RT ~ Num*Gram, data=data)
coefficients(mod)       #return coefficient vector
}

#Generate bootstrap estimate
data.boot <- boot(data, boot.huber, 1999)

#Get confidence interval
boot.ci(data.boot, …

我正在生成一个脚本,用于从cats数据集(从-MASS-包中)创建bootstrap样本.

根据Davidson和Hinkley教科书[1],我进行了一个简单的线性回归,并采用了一个基本的非参数程序,用于从iid观察引导,即对重新采样.

原始样本的形式如下:

Bwt   Hwt

2.0   7.0
2.1   7.2

...

1.9    6.8

通过单变量线性模型,我们想通过他们的大脑重量来解释猫的重量.

代码是:

library(MASS)
library(boot)


##################
#   CATS MODEL   #
##################

cats.lm <- glm(Hwt ~ Bwt, data=cats)
cats.diag <- glm.diag.plots(cats.lm, ret=T)


#######################
#   CASE resampling   #
#######################

cats.fit <- function(data) coef(glm(data$Hwt ~ data$Bwt)) 
statistic.coef <- function(data, i) cats.fit(data[i,]) 

bootl <- boot(data=cats, statistic=statistic.coef, R=999)

现在假设存在一个聚类变量cluster = 1, 2,..., 24(例如,每只猫属于一个给定的垃圾).为简单起见,假设数据是平衡的:我们对每个簇有6个观察值.因此,24窝中的每一只由6只猫组成(即n_cluster = 6和n = 144).

可以通过以下方式创建伪cluster变量:

q <- …

我想加快我的引导功能,它本身运行得非常好.我读过,因为R 2.14有一个名为的软件包parallel,但我觉得这对某人很难.对计算机科学知识不足以真正实现它.也许有人可以提供帮助.

所以这里有一个引导程序:

n<-1000
boot<-1000
x<-rnorm(n,0,1)
y<-rnorm(n,1+2*x,2)
data<-data.frame(x,y)
boot_b<-numeric()
for(i in 1:boot){
  bootstrap_data<-data[sample(nrow(data),nrow(data),replace=T),]
  boot_b[i]<-lm(y~x,bootstrap_data)$coef[2]
  print(paste('Run',i,sep=" "))
}

目标是使用并行处理/利用我的PC的多个核心.我正在Windows下运行R. 谢谢!

编辑(诺亚回复后)

以下语法可用于测试:

library(foreach)
library(parallel)
library(doParallel)
registerDoParallel(cores=detectCores(all.tests=TRUE))
n<-1000
boot<-1000
x<-rnorm(n,0,1)
y<-rnorm(n,1+2*x,2)
data<-data.frame(x,y)
start1<-Sys.time()
boot_b <- foreach(i=1:boot, .combine=c) %dopar% {
  bootstrap_data<-data[sample(nrow(data),nrow(data),replace=T),]
  unname(lm(y~x,bootstrap_data)$coef[2])
}
end1<-Sys.time()
boot_b<-numeric()
start2<-Sys.time()
for(i in 1:boot){
  bootstrap_data<-data[sample(nrow(data),nrow(data),replace=T),]
  boot_b[i]<-lm(y~x,bootstrap_data)$coef[2]
}
end2<-Sys.time()
start1-end1
start2-end2
as.numeric(start1-end1)/as.numeric(start2-end2)

但是,在我的机器上,简单的R代码更快.这是并行处理的已知副作用之一,即它会导致开销过程,这会增加像这样的"简单任务"中的时间吗？

编辑:在我的机器上,parallel代码比"简单"代码长约5倍.当我增加任务的复杂性(例如增加boot或n)时,这个因素显然不会改变.那么代码或我的机器可能存在问题(基于Windows的处理？).

我正在尝试有效地实现块引导技术来获得回归系数的分布.主要内容如下.

我有一个面板数据集,并说公司和年份是指数.对于bootstrap的每次迭代,我希望对n个主题进行替换.从这个样本中,我需要构建一个新的数据框,它是rbind()每个采样主题的所有观察的堆栈,运行回归,并拉出系数.重复一堆迭代,比如说100.

• 每个公司都可能被多次选中,因此我需要在每个迭代的数据集中多次包含它.
• 使用循环和子集方法,如下所示,似乎计算繁琐.
• 请注意,对于我的实际数据帧,n和数字迭代比下面的示例大得多.

我的想法最初是使用split()命令按主题将现有数据框分成列表.从那里,使用

sample(unique(df1$subject),n,replace=TRUE)

获取新列表,然后可能quickdf从plyr包中实现构造新的数据框.

示例慢代码:

require(plm)
data("Grunfeld", package="plm")

firms = unique(Grunfeld$firm)
n = 10
iterations = 100
mybootresults=list()

for(j in 1:iterations){

  v = sample(length(firms),n,replace=TRUE)
  newdata = NULL

  for(i in 1:n){
    newdata = rbind(newdata,subset(Grunfeld, firm == v[i]))
  }

  reg1 = lm(value ~ inv + capital, data = newdata)
  mybootresults[[j]] = coefficients(reg1)

}

mybootresults = as.data.frame(t(matrix(unlist(mybootresults),ncol=iterations)))
names(mybootresults) = names(reg1$coefficients)
mybootresults

  (Intercept)      inv    capital
1    373.8591 6.981309 -0.9801547
2    370.6743 6.633642 …

请我出去!我感谢任何帮助!谢谢!

重复进行1000次重复采样我遇到了麻烦.我尝试使用replicate()来做到这一点,但它不起作用.有没有其他方法可以做到这一点？任何人都可以告诉我,如果这可能通过使用lapply完成？以下是我的代码:

#sampling 1000 betas0 & 1 (coefficients) from the data
get.beta=function(data,indices){ 
  data=data[indices,] #let boot to select sample
  lm.out=lm(y ~ x,data=data)
  return(lm.out$coefficients)
}
n=nrow(data)
get.beta(data,1:n)

bootcoe=boot(data,get.beta,R=1000) #generate 1000 random samples
head(bootcoe$t) #look at the betas

从上面的代码我可以通过随机抽样数据得到1000 betas0&1.而且我想做1000次以获得不同的测试版.除了replicate()之外我该怎么做呢？

我正在阅读有关交叉验证以及如何使用它来选择最佳模型和估计参数的内容，但我并没有真正理解它的含义。

假设我建立一个线性回归模型并进行 10 倍交叉验证，我认为这 10 个中的每一个都会有不同的系数值，现在我应该从 10 个不同的系数中选择哪个作为我的最终模型或估计参数。

或者我们使用交叉验证只是为了找到平均误差（在我们的例子中是 10 个模型的平均值）并与另一个模型进行比较？

我有一个包含自变量和一组因变量的数据集.我想使用自举非线性最小二乘过程为每组自变量拟合一个函数.在某些情况下,自变量是"良好的质量",即合理地适合函数.在其他情况下,他们很吵.

在所有情况下,我都可以用来nls()估计参数.但是,当数据有噪声时,引导程序会抛出错误Error in nls(...) : singular gradient.我可以理解为什么nls适合噪声数据会失败,例如在太多次迭代后无法收敛,但我不明白为什么它是一个奇异的梯度误差,以及为什么我只得到质量差的重采样数据集.

码:

require(ggplot2)
require(plyr)
require(boot)

# Data are in long form: columns are 'enzyme', 'x', and 'y'
enz <- read.table("http://dl.dropbox.com/s/ts3ruh91kpr47sj/SE.txt", header=TRUE)

# Nonlinear formula to fit to data
mmFormula <- formula(y ~ (x*Vmax) / (x + Km))

nls完全能够拟合数据(即使在某些情况下a,我怀疑模型是否适合数据.

# Use nls to fit mmFormula to the data - this works well enough
fitDf <- ddply(enz, .(enzyme), function(x) coefficients(nls(mmFormula, x, start=list(Km=100, Vmax=0.5))))

# Create points to …

我有这个值数组:

> df
[1] 2 0 0 2 2 0 0 1 0 1 2 1 0 1 3 0 0 1 1 0 0 0 2 1 2 1 3 1 0 0 0 1 1 2 0 1 3
[38] 1 0 2 1 1 2 2 1 2 2 2 1 1 1 2 1 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 0 0 0 0 0
[75] …

最重要的是,我正在寻找一种快速(呃)方式对矩阵进行子集化/索引很多次:

for (i in 1:99000) {
  subset.data <- data[index[, i], ]
}

背景:
我正在实现一个涉及R中引导程序的顺序测试程序.想要复制一些模拟结果,我遇到了需要完成大量索引的瓶颈.为了实现块引导,我创建了一个索引矩阵,我用它来对原始数据矩阵进行子集以绘制数据的重采样.

# The basic setup

B <- 1000 # no. of bootstrap replications
n <- 250  # no. of observations
m <- 100  # no. of models/data series

# Create index matrix with B columns and n rows.
# Each column represents a resampling of the data.
# (actually block resamples, but doesn't matter here).

boot.index <- matrix(sample(1:n, n * B, replace=T), nrow=n, ncol=B)

# Make matrix …