在 Matrix 中获得 Bootstrap 个结果

Question

我写了下面的代码。

library(quantreg)

# return the g function:
G = function(m, N, gamma) {
  Tm = m * N
  k = 1:Tm
  Gvalue = sqrt(m) * (1 + k/m) * (k/(m + k))^gamma
  return(Gvalue)
}



sqroot <- function(A) {
  e = eigen(A)
  v = e$vectors
  val = e$values
  sq = v %*% diag(sqrt(val)) %*% solve(v)
  return(t(sq))
}

fa = function(m, N, a) {
  Tm = m * N
  k = 1:Tm
  t = (m + k)/m
  f_value = (t - 1) * t * (a^2 + log(t/(t - 1)))
  return(sqrt(f_value))
}

m = 50
N = 2
n= 50*3
x1 = matrix(runif(n, 0, 1), ncol = 1)
x = cbind(1, x1)
beta = c(1, 1)
xb = x %*% beta
pr = 1/(1+exp(-xb))  
y = rbinom(n,1,pr)

# calculate statistic:
stat = function(y, x, m, N, a) {
  y_train = y[1:m]
  x_train = x[(1:m),]
  y_test = y[-(1:m)]
  x_test = x[-(1:m),]
  
  fit = glm(y ~ 0 + x, family="binomial")
  coef = coef(fit)
  log_predict = predict(fit, type="response")
  sigma = sqrt(1/(m-1)* sum((y_train - log_predict)^2))
  
  Jvalue = t(x_train) %*% x_train/m * sigma^2
  Jsroot = sqroot(Jvalue)
  
  fvalue = fa(m, N, a)
  score1 = apply((x_test * as.vector((y_test - x_test  %*% coef))), 2, cumsum)
  statvalue1 = t(solve(Jsroot) %*% t(score1))/fvalue/sqrt(m)
  statmax1 = pmax(abs(statvalue1[, 1]), abs(statvalue1[, 2]))
  
  result = list(stat = statmax1)
  return(result)
}
m =50
N = 2
a = 2.795
value = stat(y, x, m, N, a)
value

我要执行bootstrap得到B=999个数的统计。我使用以下 r 代码。但它会产生一个错误，说“统计错误（数据，原始，...）：缺少参数“m”，没有默认值

library(boot)
data1 = data.frame(y = y, x = x1, m = m , N = N, a = a)
head(data1)
boot_value = boot(data1, statistic = stat, R = 999)

谁能给我提示？另外，我能否以矩阵格式获得 bootstrap 结果？由于 stat function 给出 100 值。

Answer 1

有多种 bootstrapping。如果您想从您的数据中提取 999 个具有相同大小的数据副本的样本，您可以只使用 replicate，不需要包。

我们将要重采样的数据放到一个数据框中。在我看来 m, N, a 保持不变，所以我们只是将其作为向量提供。

data2 <- data.frame(y=y, x=x)

stat 函数需要修改以解包 y 和 x-矩阵。在底部，我们删除了 list 调用以仅返回一个向量。 unnameing 只会给我们数字。

stat2 <- function(data, m, N, a) {
  y_train <- data[1:m, 1]
  x_train <- as.matrix(data[1:m, 2:3])
  y_test <- data[-(1:m), 1]
  x_test <- as.matrix(data[-(1:m), 2:3])
  y <- data[, "y"]
  x <- as.matrix(data[, 2:3])
  fit <- glm(y ~ 0 + x, family="binomial")
  coef <- coef(fit)
  log_predict <- predict(fit, type="response")
  sigma <- sqrt(1/(m-1) * sum((y_train - log_predict)^2))
  Jvalue <- t(x_train) %*% x_train/m * sigma^2
  Jsroot <- sqroot(Jvalue)
  fvalue <- fa(m, N, a)
  score1 <- apply((x_test * as.vector((y_test - x_test %*% coef))), 2, cumsum)
  statvalue1 <- t(solve(Jsroot) %*% t(score1))/fvalue/sqrt(m)
  statmax1 <- pmax(abs(statvalue1[, 1]), abs(statvalue1[, 2]))
  result <- unname(statmax1)
  return(result)
}

replicate 是 sapply 的表亲，专为重复计算而设计。在调用中，我们只是 sample 行 999 次，并且已经得到一个矩阵。与 sapply 一样，我们需要 t 转换我们的结果。

res <- t(replicate(999, stat2(data2[sample(1:nrow(data2), nrow(data2), replace=TRUE), ], m, N, a)))

结果

结果我们在行中得到 999 bootstrap 个复制，在列中得到 100 个属性。

str(res)
# num [1:999, 1:100] 0.00205 0.38486 0.10146 0.12726 0.47056 ...

代码运行速度也相当快。

user  system elapsed 
3.46    0.01    3.49

注意，有不同种类的 bootstrapping。例如。有时只对样本的一部分进行重新采样，使用权重，应用聚类等。由于您尝试使用 boot，因此显示的方法应该是默认方法。

在 Matrix 中获得 Bootstrap 个结果

Obtain Bootstrap Results in Matrix

r

resampling

statistics-bootstrap

结果