Оценка плотности ядра Вейбулла в R
Я пытаюсь оценить работу 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")
Можете добавить ожидаемый результат? Мы не знаем, как должен выглядеть сюжет.
Ответы 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) жестко запрограммирована, чтобы быть эквивалентной полосе пропускания исследовательских работ, однако это не критическое исправление.
- Для индексации цикла
Цикл 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")
но они не похожи на требуемый. Даже масштаб по оси x и оси y совершенно разные.
К сожалению, в вопросе Q неясно, что касается входных данных - в исследовательской статье обнаружено использование журнала yy
Я просто использовал пропускную способность (0,480411), которая указана в документе, чтобы обеспечить точно такие же результаты.
Отлично :) --- @Angel, можешь ли ты отметить ответ?
Почему вы проголосовали против? объясните пожалуйста причину.