标签: statistics-bootstrap

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

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 …

下面是我的问题的MWE:我已经使用引导程序(通过引导程序包中的引导功能)为某些功能编写了进度条.

只要我不使用并行处理(res_1core下面),这样就可以正常工作.如果我想通过设置parallel = "multicore"和使用并行处理ncpus = 2,则进度条显示不正确(res_2core如下).

library(boot)

rsq <- function(formula, data, R, parallel = c("no", "multicore", "snow"), ncpus = 1) {
  env <- environment()
  counter <- 0
  progbar <- txtProgressBar(min = 0, max = R, style = 3)
  bootfun <- function(formula, data, indices) {
    d <- data[indices,]
    fit <- lm(formula, data = d)
    curVal <- get("counter", envir = env)
    assign("counter", curVal + 1, envir = env)
    setTxtProgressBar(get("progbar", envir = env), curVal + …

我正在尝试向 R 中的引导函数添加一个进度条。我试图使示例函数尽可能简单（因此我在本示例中使用了 mean）。

library(boot)
v1 <- rnorm(1000)
rep_count = 1

m.boot <- function(data, indices) {
  d <- data[indices]
  setWinProgressBar(pb, rep_count)
  rep_count <- rep_count + 1
  Sys.sleep(0.01)
  mean(d, na.rm = T) 
  }

tot_rep <- 200
pb <- winProgressBar(title = "Bootstrap in progress", label = "",
                     min = 0, max = tot_rep, initial = 0, width = 300)
b <- boot(v1, m.boot, R = tot_rep)
close(pb)

引导程序正常运行，但问题是rep_count循环中的值没有增加，并且在此过程中进度条保持冻结状态。

如果我rep_count在引导完成后检查 的值，它仍然是 1。

我究竟做错了什么？也许引导函数不是简单地将m.boot函数插入到循环中，因此其中的变量不会增加？

谢谢你。

这个问题的答案的底部（使用 R 中的约束计算固定效应的 CI）建议人们应该看到user时间 >elapsed并行工作时的时间。尽管parallel = "multicore", ncpus = 4在运行时指定了boot.ci我没有看到那个结果。此外，我在 Mac 的活动监视器运行时只看到大约 30% 的 CPU 负载。这是否意味着我不能与我的 4 核 iMac 进行并行处理？如果没有，关于让它工作的任何建议？

下面是一个例子：

library(car)
library(boot)
set.seed(47)

 y <- rgamma(2000, 2)
 x1 <- 3 * y + rnorm(2000)
 x2 <- y^2 + rnorm(2000)
 x3 <- rnorm(2000)
 MyData <- data.frame(c(y, x1, x2, x3))
 MyModel <- lm(y ~ x1 + x2 + x3, data = MyData)
# Boot doesn't have a parallel option that I …

问题:如何使用boostrap来获得在协方差矩阵的特征值上计算的统计数据的置信区间,分别针对数据框中的每个组(因子级别)？

问题:我无法确定我需要包含适合该boot函数的这些结果的数据结构,或者是一种在组中"映射"引导程序并以适合绘图的形式获得置信区间的方法.

上下文:在heplots包中,boxM计算协方差矩阵相等的Box的M检验.有一种绘图方法可以生成进入此测试的对数决定因素的有用图.该图中的置信区间基于渐近理论近似.

> library(heplots)
> iris.boxm <- boxM(iris[, 1:4], iris[, "Species"])
> iris.boxm

        Box's M-test for Homogeneity of Covariance Matrices

data:  iris[, 1:4]
Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16

> plot(iris.boxm, gplabel="Species")

绘图方法还可以显示特征值的其他函数,但在这种情况下没有可用的理论置信区间.

op <- par(mfrow=c(2,2), mar=c(5,4,1,1))
plot(iris.boxm, gplabel="Species", which="product")
plot(iris.boxm, gplabel="Species", which="sum")
plot(iris.boxm, gplabel="Species", which="precision")
plot(iris.boxm, gplabel="Species", which="max")
par(op)

因此,我希望能够使用boostrap计算这些CI,并将其显示在相应的图中.

我尝试过的:

下面是提升这些统计数据的函数,但对于总样本,不考虑group(Species).

cov_stat_fun <- function(data, indices, 
            stats=c("logdet", "prod", "sum", "precision", "max")
            ) …

我试图从数据集中收集一些摘要统计数据的自举估计值,但我想以不同的速率重新采样部分数据集,这使我依赖于嵌套的for循环.

具体来说,假设我的数据集中有两个组,每个组进一步分为测试和控制.第1组具有75%/ 25%的测试控制比,第2组具有50%/ 50%的测试控制比.

我想重新采样,使得数据集大小相同,但两组的测试控制比率均为90%/ 10%...换句话说,以不同的速率对不同的子组进行重新采样,这让我感觉不同于boot包通常会.

在我的数据集中,我创建了一个group表示组的groupT变量,以及一个表示与test/control连接的组的变量,例如:

    id     group     groupT
     1         1         1T
     2         1         1T
     3         2         2T
     4         1         1C
     5         2         2C

这是我现在正在运行的,nreps任意设置为我的引导复制数:

for (j in 1:nreps){

  bootdat <- datafile[-(1:nrow(datafile)),] ## initialize empty dataset

  for (i in unique(datafile$groups)){

    tstring<-paste0(i,"T") ## e.g. 1T
    cstring<-paste0(i,"C") ## e.g. 1C

    ## Size of test group resample should be ~90% of total group size

    tsize<-round(.90*length(which(datafile$groups==i)),0)

    ## Size of control group resample should be total …

在python中是否有一个等效的boot和boot.ci？在R我会这样做

library(boot)
result <- boot(data,bootfun,10000)
boot.ci(result)

我的目的是编写几个函数，旨在找到两个协方差矩阵之间的整体相似性，方法是将它们与随机向量相乘并关联响应向量，或者通过自举矩阵之一来获得可用于比较的相关系数分布。但在这两种情况下，我都得到了错误的结果。观察到的矩阵间相关性高达 0.93，但分布最多仅达到 0.2。这是函数的代码：

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{
  statSim <- numeric(numR)
  mat1vcv <- cov(mat1)
  mat2vcvT <- cov(mat2)
  ltM1 <- mat1vcv[col(mat1vcv) <= row(mat1vcv)]
  ltM2T <- mat2vcvT[col(mat2vcvT) <= row(mat2vcvT)]
  statObs <- cor(ltM1, ltM2T)                           
  indice <- c(1:length(mat2))
  resamplesIndices <- lapply(1:numR, function(i) sample(indice, replace = F))
  for (i in 1:numR)
  {
    ss <- mat2[sample(resamplesIndices[[i]])]
    ss <- matrix(ss, nrow = dim(mat2)[[1]], ncol = dim(mat2)[[2]])
    mat2ss <- cov(ss)
    ltM2ss <- mat2ss[col(mat2ss) <= row(mat2ss)]
    statSim[i] <- cor(ltM1, ltM2ss)
  }
  if (graph == TRUE) …

我创建了一个基于似然函数和模拟的概率模拟，所有这些都可以用下面的代码复制。

这是似然函数：

probit.ll <- function(par,ytilde,x) {
   a <- par[1] 
   b <- par[2]
 return(  -sum( pnorm(ytilde*(a + b*x),log=TRUE) ))
}

这是进行估计的函数：

my.probit <- function(y,x) {
# use OLS to get start values
par <- lm(y~x)$coefficients
ytilde <- 2*y-1
# Run optim 
res <- optim(par,probit.ll,hessian=TRUE,ytilde=ytilde,x=x)
# Return point estimates and SE based on the inverse of Hessian
names(res$par) <- c('a','b')
se=sqrt(diag(solve(res$hessian)))
names(se) <- c('a','b')
return(list(par=res$par,se=se,cov=solve(res$hessian)))
}

这是生成模拟模型的函数：

probit.data <- function(N=100,a=1,b=1) {
x <- rnorm(N)
y.star <- a + b*x + rnorm(N)
y …

library(boot)
set.seed(1)
x=sample(0:1000,1000)
y=function(u,i) sum(x[i])
o=boot(x,y,1000)
theta1=NULL
theta1=cbind(theta1,o$t)
b=theta1[order(theta1)]
bp1=c(b[25], b[975])
ci=boot.ci(o,type="perc")

我使用两种方法来构建引导百分位数置信区间，但我得到了两个不同的答案。

bp1=c(b[25], b[975]) get (480474,517834)

同时ci=boot.ci(o,type="perc")得到 (480476, 517837 )

boot.ci 如何构建百分位区间？