Подбор данных методом максимальной вероятности для нового распределения

Question

Я хотел бы знать, как можно приспособить любое распределение к заданному набору данных, используя метод MLE.В качестве конкретного примера, может ли кто-нибудь предложить рабочий код, который даст правильные результаты для MLE для $ \ theta $ и $ \ beta $, когда обобщенное распределение Линдли описано в https://rivista -statistica.unibo.it /article / viewFile / 6836/7039 применяется к данным: 5.1, 6.3, 10.8, 12.1, 18.5, 19.7, 22.2, 23, 30.6, 37.3, 46.3, 53.9, 59.8, 66.2 на стр.156?Как это можно использовать для распределения по гистограмме?

Answer 1

Я собираюсь ответить на ваш теперь удаленный вопрос, который очень похож, основываясь на проблеме в Парари и др. (данные по кондиционированию воздуха от Proschan (Таблица 7.2, результаты в Таблице7.5)

В общем, вам нужно знать разумные начальные значения, чтобы можно было выполнить оценку максимального правдоподобия общего назначения (это общеизвестно, что проблема курицы и яйца, но есть много решений [используйте метод моментов, выбирайте разумные значения, основанные на проблеме, просматривайте график для выбора значений и т. д.]. В этом случае, поскольку авторы дали вам ответы, вы можете использовать эти параметры, чтобы выбрать начальные значения правильного порядка.Чтобы начать, вам нужно знать кое-что о разумных порядках.

Обратите внимание, что существует лотов суетливых числовых проблем с общей оценкой максимального правдоподобия: например,см. главу 7 здесь ...

Данные (x) и функция правдоподобия (dBEPL) являются определенныминиже.Я определяю функцию плотности и использую интерфейс формулы для mle2() ...

dd  <- data.frame(x)
library(bbmle)
## parameters from table
ovec <- list(alpha=.7945,beta=0.1509,omega=6.7278,
              a=0.2035,b=0.2303)
## starting vals -- right order of magnitude
svec <- list(alpha=0.8,beta=0.2,omega=10,a=0.2,b=0.2)
m1 <- mle2( x ~ dBEPL(alpha,beta,omega,a,b),
      data=dd,
      start=svec,
      control=list(parscale=unlist(svec)))
coef(m1)  ## NOT the same as reported, but close
par(las=1,bty="l")
hist(x,col="gray",freq=FALSE,breaks=40,ylim=c(0.,0.014))
with(as.list(coef(m1)),curve(dBEPL(x,alpha,beta,omega,a,b,log=FALSE),add=TRUE,
                             col=2,lwd=2,n=201))

---

x <- scan(textConnection("
194 413 90 74 55 23 97 50 359 50 130 487 57 102 
15 14 10 57 320 261 51 44 9 254 493 33 18 209 41 
58 60 48 56 87 11 102 12 5 14 14 29 37 186 29 104 
7 4 72 270 283 7 61 100 61 502 220 120 141 22 603 35 
98 54 100 11 181 65 49 12 239 14 18 39 3 12 5 32 9 438 
43 134 184 20 386 182 71 80 188 230 152 5 36 79 59 33 
246 1 79 3 27 201 84 27 156 21 16 88 130 14 118 44 15 42 
106 46 230 26 59 153 104 20 206 566 34 29 26 35 5 82 31 
118 326 12 54 36 34 18 25 120 31 22 18 216 139 67 310 
346 210 57 76 14 111 97 62 39 30 7 44 11 63 23 22 23 14 
18 13 34 16 18 130 90 163 208 124 70 16 101 52 
208 95 62 11 191 14 71"))


dBEPL <- function(x,alpha,beta,omega,a,b,log=TRUE) {
    r <- log(alpha*beta^2*omega/(beta(a,b)*(beta+1)))+
           log(1+x^alpha)+(alpha-1)*log( x)-beta* x^alpha+(omega*a-1) *
           log(1-(1+beta* x^alpha/(beta+1))*exp(-beta* x^alpha))+
           (b-1)*log(1-(1-(1+beta* x^alpha/(beta+1))*exp(-beta* x^alpha))^omega)
    if (log) return(r) else return(exp(r))
}