Постройте вероятность вейбулла

Я хотел бы построить функцию правдоподобия для образца вейбулла размером 1000 с последовательностью параметров формы тета. Я использовал стандартизированный weibull, поэтому лямбда шкалы равна 1. Однако на выходе получается горизонтальная прямая линия.

n<-1000
lik <- function(theta, x){
  K<- length(theta)
  n<- length(x)
  out<- rep(0,K)
  for(k in 1:K){
    out[k] <- prod(dweibull(x, shape= theta[k], scale=1))   
  }
  return(out)
}
theta<-seq(0.01, 10, by = 0.01)
x <- rweibull(n, shape= 0.5, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)

r plot statistics weibull
person siegfried    schedule 23.05.2020    source источник
comment
Попробуйте построить логарифмическую вероятность; sum(dweibull(x, shape= theta[k], scale=1, log=TRUE)) и держите тета в более правдоподобном диапазоне seq(0.01, 1, by = 0.01)   -  person user20650    schedule 23.05.2020
comment
@ user20650, напишите как ответ?   -  person Ben Bolker    schedule 23.05.2020

Ответы (1)


arrow_upward
2
arrow_downward

В том, что вы сделали, нет ничего на самом деле неправильного, но компьютеры с трудом вычисляют произведение многих маленьких чисел и поэтому могут получить нуль (даже 0.99^1000 = 4^-5). А так проще log трансформировать а потом sum. (Поскольку логарифмическое преобразование является монотонной возрастающей функцией, максимизация логарифмической вероятности аналогична максимизации правдоподобия). Таким образом, измените

prod(dweibull(x, shape= theta[k], scale=1))

to 

sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))

Другое незначительное изменение состоит в том, чтобы построить график вероятности в разумном диапазоне theta, чтобы вы могли видеть кривую.

Рабочий код:

set.seed(1)
n<-1000
lik <- function(theta, x){
  K <- length(theta)
  n <- length(x)
  out <- rep(0,K)
  for(k in 1:K){
    out[k] <- sum(dweibull(x, shape= theta[k], scale=1, log=TRUE))   
  }
  return(out)
}

popTheta = 0.5
theta = seq(0.01, 1.5, by = 0.01)
x = rweibull(n, shape=popTheta, scale= 1)
plot(theta, lik(theta, x), type="l", lwd=2)
abline(v=popTheta)

theta[which.max( lik(theta, x))]
person user20650    schedule 23.05.2020