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

希望对您有所帮助：

要得到上面的图需要解决两个问题：

使用对数

首先（在您所指的 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) 被硬编码为等同于研究论文带宽，但这不是关键修复。

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")

R 中的威布尔核密度估计

Weibull Kernel Density Estimation in R

r

kernel-density

density-plot