如何计算核密度及其累积分布函数的积分？

Question

我正在尝试在 R 中计算以下积分： 在此输入图片描述 在哪里 在此输入图片描述 表示累积分布函数的核密度函数和核估计，h 是带宽参数。我想计算上面的估计量。我在 R 中编写了以下代码：

JXE <- function(t, n, a) {

z <- runif(n)

f  <- approxfun(density(z)$x, density(z)$y)

fhat <- function(t) ifelse(is.na(f(t)), 0, f(t))

Fhat <-function(t) integrate(fhat, lower = 0, upper = t)

fhat_kn<-function(x,t){

(-log(Fhat(x)/Fhat(t)))^(n-1)*(fhat(x)/Fhat(t))

}

Hhatkn<-function(t,a){

integrate(function(x) (fhat_kn(x,t))^a,lower=0,upper=t)$value

}

return(Hhatkn(t,a))

}

n=5

a=0.2

t=0.5

JXE(t,n,a)

但是，我有以下错误：

Error in integrate(fhat, lower = 0, upper = t) : 
  length(upper) == 1 is not TRUE

Is there anyone to help me what is the problem?
Thanks a lot for your valuable comments!

我有一个问题，如果可能的话需要收到我的答复。

Answer 1

从评论中添加，我们可以像这样继续，用一个最小的例子：

# Install and load the kerdiest package if not already installed
if (!require(kerdiest)) { remotes::install_github("cran/kerdiest") }
library(kerdiest)

JXE <- function(t, n, a, bw, seed = 123) {
  # simulate data
  set.seed(seed) # to replicate results
  z <- runif(n, 0, t)

  # kernel density and cdf estimation using kerdiest package
  kde_result <- kde(vec_data = z, bw = bw)
  f <- approxfun(kde_result$grid, kde_result$Estimated_values)

  # cumulative distribution function estimation
  Fhat <- function(t) {
    result <- tryCatch(
      integrate(function(x) f(x), lower = 0, upper = t)$value,
      error = function(e) NA # in case of integration errors
    )
    return(result)
  }

  # density estimation function
  fhat <- function(t) ifelse(is.na(f(t)), 0, f(t))

  # density and cdf estimate
  fhat_kn <- function(x, t) {
    Fhat_t <- Fhat(t) # avoid repeated computation
    Fhat_x <- Fhat(x)
    if (is.na(Fhat_t) || is.na(Fhat_x) || Fhat_t == 0 || Fhat_x == 0) {
      return(0)
    }
    value <- (-log(Fhat_x / Fhat_t))^(n - 1) * (fhat(x) / Fhat_t)
    if (is.infinite(value) || is.nan(value)) {
      return(0)
    }
    return(value)
  }

  # final integral
  Hhatkn <- function(t, a) {
    tryCatch(
      integrate(function(x) (fhat_kn(x, t))^a, lower = 0, upper = t)$value,
      error = function(e) 0 # handle integration errors
    )
  }

  # debug
  cat("Fhat(t):", Fhat(t), "\n")
  cat("fhat_kn values:\n")
  for (x in seq(0, t, length.out = 10)) {
    cat("x:", x, "fhat_kn:", fhat_kn(x, t), "\n")
  }

  Hhatkn(t, a)
}

# parameters
n <- 1000
a <- 0.2
t <- 0.5
bw <- 0.1

result <- JXE(t, n, a, bw, seed = 123)

cat("Result:", result, "\n")

如果清楚请告诉我。

我编辑了问题以添加评论。