R中的威布尔核密度估计

Question

我正在尝试评估 Salha 等人的工作。 (2014)，题为“使用 Weibull 内核的危险率函数估计”。但我的密度图（真实数据）只是一条平线，而不是类似于文章的适当密度图。这是预期的密度、我的 R 代码和 Weibull 内核。请帮助我找出我的错误。

R 代码：

k<-200
yy<-c(1,1,1,5,7,8,8,13,14,14,17,18,21,21,22,25,27,27,30,30,31,31,32,34,35,36,37,38,39,39,40,49,49,54,56,56,62,63,65,65,67,75,76,79,82,83,84,84,84,90,91,92,93,93,103,103,111,112,119,122,123,126,129,134,144,147,153,163,167,175,228,231,235,242,256,256,257,311,314,322,369,415,573,609,640,737)
y<-log(yy)
n<-length(yy)

 h<-0.79 * IQR(y) * length(y) ^ (-1/5)
 x <- seq(min(yy) + 0.05, max(yy), length = k)

 KWeibull <- matrix(rep(0, k * n), ncol = k)
 fhat <- rep(0, k)
###########weibull###########
for (j in 1:k) {
    for (i in 1:n) {
        fn <- gamma(1 + h)
        KWeibull[i, j] <- (fn/(h * x[i])) * ((yy[i] * fn)/x[i])^((1/h) - 1) * exp(-((yy[i] * 
            fn)/x[i])^(1/h))
    }
    fhat[j] <- 1/n * (sum(KWeibull[, j]))
}

plot(x,fhat, type = "l")

Answer 1

希望这会有所帮助：

要实现上述情节，需要解决两个问题：

1. 使用对数

首先（在您所指的paper 中）他们使用输入数据的对数 - 我在论文的第 5.2 节中找到了这一点 --- 下面是这个修复：

k<-200
yy<-c(1,1,1,5,7,8,8,13,14,14,17,18,21,21,22,25,27,27,30,30,31,31,32,34,35,36,37,38,39,39,40,49,49,54,56,56,62,63,65,65,67,75,76,79,82,83,84,84,84,90,91,92,93,93,103,103,111,112,119,122,123,126,129,134,144,147,153,163,167,175,228,231,235,242,256,256,257,311,314,322,369,415,573,609,640,737)
y<-log(yy)
n<-length(yy)

#h<-0.79 * IQR(y) * length(y) ^ (-1/5)
x <- seq(min(y) + 0.05, max(y), length = k)
h <- 0.480411

KWeibull <- matrix(rep(0, k * n), ncol = k)
fhat <- rep(0, k)

请注意，带宽 (h) 被硬编码为等同于研究论文的带宽，但这不是关键的修复。

1. For 循环索引

for 循环 - 您正在迭代您的 yy （我认为这是您在内核密度估计器中的时间变量），但您的 x 随机样本每次都在相同的集合上进行迭代。并且还使用 y 而不是 yy，因为这是对数转换后的数据。

请参阅以下修复：（这包括对数修复）

###########weibull###########
for (j in 1:k) {
  for (i in 1:n) {
    fn <- gamma(1 + h)
    KWeibull[i, j] <- (fn/(h * x[j])) * ((y[i] * fn)/x[j])^((1/h) - 1) * exp(-((y[i] * 
                                                                                   fn)/x[j])^(1/h))
  }
  fhat[j] <- 1/n * (sum(KWeibull[, j]))
}

plot(yy,KWeibull[,86], type = "l")
plot(x,fhat, type = "l")