R 中的 Jackknife 获取区间估计值

Question

我有一个关于如何使用

jackknife

包使用

bootstrap

的问题。我想获得折刀法的区间估计。

我尝试运行下面的代码，但没有得到参数估计的结果。

rm(list=ls())
library(bootstrap)
library(maxLik)

set.seed(20)
lambda <- 0.02
beta <- 0.5
alpha <-  0.10
n <- 40
N <- 1000
lambda_hat <- NULL
beta_hat <- NULL
cp <- NULL
jack_lambda <- matrix(NA, nrow = N, ncol = 2)
jack_beta <- matrix(NA, nrow = N, ncol = 2)

### group all data frame generated from for loop into a list of data frame
data_full <- list()
for(i in 1:N){
  u <- runif(n)
  c_i <- rexp(n, 0.0001)
  t_i <- (log(1 - (1 / lambda) * log(1 - u))) ^ (1 / beta)
  s_i <- 1 * (t_i < c_i)
  t <- pmin(t_i, c_i)
  data_full[[i]] <- data.frame(u, t_i, c_i, s_i, t)
}

### statistic function for jackknife()
estjack <- function(data, j) {
  data <- data[j, ]
  data0 <- data[which(data$s_i == 0), ] #uncensored data
  data1 <- data[which(data$s_i == 1), ] #right censored data
  data
 
  LLF <- function(para) {
    t1 <- data$t_i
   
    lambda <- para[1]
    beta <- para[2]
    e <- s_i*log(lambda*t1^(beta-1)*beta*exp(t1^beta)*exp(lambda*(1-exp(t1^beta))))
    r <- (1-s_i)*log(exp(lambda*(1-exp(t1^beta))))
    f <- sum(e + r)
    return(f)
  }
  mle <- maxLik(LLF, start = c(para = c(0.02, 0.5)))
  lambda_hat[i] <- mle$estimate[1]
  beta_hat[i] <- mle$estimate[2]
  return(c(lambda_hat[i], beta_hat[i]))
}

jackknife_resample<-list()
for(i in 1:N) {
  jackknife_resample[[i]]<-data_full[[i]][-i]
  results <- jackknife(jackknife_resample, estjack,R=1000)
  jack_lambda[i,]<-lambda_hat[i]+c(-1,1)*qt(alpha/2,n-1,lower.tail = FALSE)*results$jack.se
  jack_beta[i,]<-beta_hat[i]+c(-1,1)*qt(alpha/2,n-1,lower.tail = FALSE)*results$jack.se
}```

I couldn't get the parameter estimate that run in MLE and hence couldn't proceed to the next step. Can anyone help?

Answer 1

看来您尝试自己实施留一法。这是可能的，但在这种情况下，您不需要

bootstrap::jackknife

:

jk <- numeric(N)
for (i in 1:N) {
  jk[i] <- your_statistic(data_full[-i,])
}
jk.m <- mean(jk)
jk.se <- sqrt((n-1)/n * sum((jk-jk.m)^2))

函数

bootstrap::jackknife

为您提供了留一法，但是当统计数据需要的不仅仅是一个简单的向量作为输入时，定义计算统计数据的函数有点棘手。在这种情况下，请查看

bootstrap::jackknife

文档中的最后一个示例：

# To jackknife functions of more  complex data structures, 
# write theta so that its argument x
#  is the set of observation numbers  
#  and simply  pass as data to jackknife the vector 1,2,..n. 
# For example, to jackknife
# the correlation coefficient from a set of 15 data pairs:      
                         
xdata <- matrix(rnorm(30),ncol=2)
n <- 15
theta <- function(x,xdata) {
  cor(xdata[x,1],xdata[x,2])
}
results <- jackknife(1:n,theta,xdata)