Оценка плотности ядра Вейбулла в R

Question

Я пытаюсь оценить работу Salha et al. (2014), озаглавленный « Оценка функции коэффициента опасности с использованием ядра Вейбулла ». Но мой график плотности ( для реальных данных ) - это просто плоская линия, а не правильный график плотности, подобный статье. Здесь ожидаемая плотность, мой код R и ядро ​​Вейбулла. пожалуйста, помогите мне узнать мою ошибку.

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. Использовать логарифм

Во-первых (в документе , на который вы ссылаетесь) они используют логарифм входных данных - я нашел это в разделе 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) жестко задана, чтобы быть эквивалентной полосе пропускания исследовательских работ, однако это не является решающим исправлением.

2. Для индексации цикла

Цикл 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")