Jackknife en R para obtener intervalos de estimación

06 Aug, 2023 Programación 0

Tengo una pregunta sobre cómo utilizar el jackknife usando el paquete bootstrap. Quiero obtener el intervalo de estimación para el método del jackknife.

He intentado ejecutar el código a continuación, pero no obtengo resultados para mi estimación de parámetros.

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

set.seed(20)
lambda <- 0.02
beta <- 0.5
alpha <- 0.10
n <- 40
N <- 1000
lambdahat <- NULL
betahat <- NULL
cp <- NULL
jacklambda <- matrix(NA, nrow = N, ncol = 2)
jackbeta <- matrix(NA, nrow = N, ncol = 2)

Agrupar todos los marcos de datos generados por el bucle en una lista de marcos de datos.

datafull <- list()
for(i in 1:N){
  u <- runif(n)
  ci <- rexp(n, 0.0001)
  ti <- (log(1 - (1 / lambda) * log(1 - u))) ^ (1 / beta)
  si <- 1 * (ti < ci)
  t <- pmin(ti, ci)
  datafull[[i]] <- data.frame(u, ti, ci, si, t)
}

Función estadística para jackknife().

estjack <- function(data, j) {
  data <- data[j, ]
  data0 <- data[which(data$si == 0), ] # datos no censurados
  data1 <- data[which(data$si == 1), ] # datos censurados a la derecha

  LLF <- function(para) {
    t1 <- data$ti

    # Código omitido: funciones y cálculos adicionales

    f <- sum(e + r)
    return(f)
  }

  mle <- maxLik(LLF, start = c(para = c(0.02, 0.5)))
  lambdahat[i] <- mle$estimate[1]
  betahat[i] <- mle$estimate[2]
  return(c(lambdahat[i], betahat[i]))
}

Crear muestras de jackknife.

jackkniferesample <- list()
for(i in 1:N) {
  jackkniferesample[[i]] <- datafull[[i]][-i]
  results <- jackknife(jackkniferesample, estjack, R=1000)
  jacklambda[i,] <- lambdahat[i] + c(-1,1)*qt(alpha/2,n-1,lower.tail = FALSE)*results$jack.se
  jackbeta[i,] <- betahat[i] + c(-1,1)*qt(alpha/2,n-1,lower.tail = FALSE)*results$jack.se
}

No pude obtener la estimación de parámetros que se ejecuta en MLE y, por lo tanto, no pude continuar con el siguiente paso. ¿Alguien puede ayudar?

davy.ai

August 6, 2023 at 3:17 am

Basado en el código proporcionado, parece que puede haber algunos problemas con la función estjack.

En primer lugar, la línea data debe ser eliminada ya que no tiene ningún efecto en la función.

En segundo lugar, parece que hay un error tipográfico en la línea que define r. Debería ser r <- (1-data$s_i)*log(exp(lambda*(data$c_i-t1^beta))) para calcular la verosimilitud logarítmica de las observaciones censuradas a la derecha.

Además, el uso de s_i debe ser reemplazado por data$s_i para hacer referencia correctamente al indicador de censura en los cálculos de verosimilitud logarítmica.

Por último, la salida de la función debe ser un vector, no una matriz. La última línea debe modificarse a return(c(lambda_hat[i], beta_hat[i])).

La función estjack corregida se muestra a continuación:

estjack <- function(data, j) {
  data <- data[j, ]
  data0 <- data[which(data$s_i == 0), ] # datos no censurados
  data1 <- data[which(data$s_i == 1), ] # datos censurados a la derecha

  LLF <- function(para) {
    t1 <- data$t_i

    lambda <- para[1]
    beta <- para[2]
    e <- data$s_i*log(lambda*t1^(beta-1)*beta*exp(t1^beta)*exp(lambda*(1-exp(t1^beta))))
    r <- (1-data$s_i)*log(exp(lambda*(data$c_i-t1^beta)))
    f <- sum(e + r)
    return(f)
  }
  mle <- maxLik(LLF, start = 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]))
}

Una vez corregido esto, puedes usar la función bootstrap() del paquete bootstrap para obtener una estimación del intervalo para el método de jackknife. El siguiente código muestra un ejemplo de cómo hacerlo:

jackknife_resample <- list()
for (i in 1:N) {
  jackknife_resample[[i]] <- data_full[[i]][-i]
}

estimates <- lapply(jackknife_resample, estjack)

lambda_hat <- sapply(estimates, "[", 1)
beta_hat <- sapply(estimates, "[", 2)

jack_lambda <- bootstrap(lambda_hat, estfun = mean, R = 1000, type = "jackknife")
jack_beta <- bootstrap(beta_hat, estfun = mean, R = 1000, type = "jackknife")

alpha <- 0.10
quantiles <- c((1-alpha/2), alpha/2)

lambda_interval <- quantile(jack_lambda$t, quantiles)
beta_interval <- quantile(jack_beta$t, quantiles)

print(paste0("Intervalo de estimación de jackknife para lambda: [", lambda_interval[1], ", ", lambda_interval[2], "]"))
print(paste0("Intervalo de estimación de jackknife para beta: [", beta_interval[1], ", ", beta_interval[2], "]"))

Este código primero calcula las estimaciones de lambda y beta utilizando la función estjack para cada muestra jackknife. Estas estimaciones se utilizan luego para calcular los errores estándar del jackknife utilizando la función bootstrap() con type = "jackknife". Por último, los intervalos de estimación se obtienen tomando los cuantiles de las estadísticas bootstrap. Tenga en cuenta que la variable alpha se utiliza para especificar el nivel de confianza del intervalo de estimación.

Al ejecutar este código, deberían obtenerse los intervalos estimados para lambda y beta.

Answer

davy.ai

August 6, 2023 at 3:17 am
Basado en el código proporcionado, parece que puede haber algunos problemas con la función estjack.

En primer lugar, la línea data debe ser eliminada ya que no tiene ningún efecto en la función.

En segundo lugar, parece que hay un error tipográfico en la línea que define r. Debería ser r <- (1-data$s_i)*log(exp(lambda*(data$c_i-t1^beta))) para calcular la verosimilitud logarítmica de las observaciones censuradas a la derecha.

Además, el uso de s_i debe ser reemplazado por data$s_i para hacer referencia correctamente al indicador de censura en los cálculos de verosimilitud logarítmica.

Por último, la salida de la función debe ser un vector, no una matriz. La última línea debe modificarse a return(c(lambda_hat[i], beta_hat[i])).

La función estjack corregida se muestra a continuación:
```
estjack <- function(data, j) {
  data <- data[j, ]
  data0 <- data[which(data$s_i == 0), ] # datos no censurados
  data1 <- data[which(data$s_i == 1), ] # datos censurados a la derecha

  LLF <- function(para) {
    t1 <- data$t_i

    lambda <- para[1]
    beta <- para[2]
    e <- data$s_i*log(lambda*t1^(beta-1)*beta*exp(t1^beta)*exp(lambda*(1-exp(t1^beta))))
    r <- (1-data$s_i)*log(exp(lambda*(data$c_i-t1^beta)))
    f <- sum(e + r)
    return(f)
  }
  mle <- maxLik(LLF, start = 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]))
}
```
Una vez corregido esto, puedes usar la función bootstrap() del paquete bootstrap para obtener una estimación del intervalo para el método de jackknife. El siguiente código muestra un ejemplo de cómo hacerlo:
```
jackknife_resample <- list()
for (i in 1:N) {
  jackknife_resample[[i]] <- data_full[[i]][-i]
}

estimates <- lapply(jackknife_resample, estjack)

lambda_hat <- sapply(estimates, "[", 1)
beta_hat <- sapply(estimates, "[", 2)

jack_lambda <- bootstrap(lambda_hat, estfun = mean, R = 1000, type = "jackknife")
jack_beta <- bootstrap(beta_hat, estfun = mean, R = 1000, type = "jackknife")

alpha <- 0.10
quantiles <- c((1-alpha/2), alpha/2)

lambda_interval <- quantile(jack_lambda$t, quantiles)
beta_interval <- quantile(jack_beta$t, quantiles)

print(paste0("Intervalo de estimación de jackknife para lambda: [", lambda_interval[1], ", ", lambda_interval[2], "]"))
print(paste0("Intervalo de estimación de jackknife para beta: [", beta_interval[1], ", ", beta_interval[2], "]"))
```
Este código primero calcula las estimaciones de lambda y beta utilizando la función estjack para cada muestra jackknife. Estas estimaciones se utilizan luego para calcular los errores estándar del jackknife utilizando la función bootstrap() con type = "jackknife". Por último, los intervalos de estimación se obtienen tomando los cuantiles de las estadísticas bootstrap. Tenga en cuenta que la variable alpha se utiliza para especificar el nivel de confianza del intervalo de estimación.

Al ejecutar este código, deberían obtenerse los intervalos estimados para lambda y beta.

Comments are closed.

es.davy.ai

¿Tienes una pregunta?

Jackknife en R para obtener intervalos de estimación

Answer

davy.ai