mle2 из пакета bbmle дает бессмысленные ответы с пользовательской функцией

Я хочу сопоставить некоторые полевые данные со специальной отрицательной экспоненциальной функцией плотности вероятности. (Мотивация: в конечном итоге я хочу подогнать полевые данные ко многим распределениям в таблице 3 Буллок Ши и Скарпаас 2006).

Сначала я определил функцию dnegexp в соответствии с этим сообщением: Ошибка с определением пользовательской функции плотности для вызова формулы mle2

dnegexp <- function(x, mya, myb, log=FALSE){
   logne <- log(mya)-myb*x
   if(log) return(logne) else return(exp(logne))
}

Затем я сделал функцию rnegexp, которая генерирует набор данных из этого распределения с любыми заданными параметрами. (Пожалуйста, простите за неэлегантность этой функции - это не тема этого поста, но идеи по ее улучшению приветствуются...)

rnegexp <- function(n, thisa, thisb){
  rnums <- array(data=NA,dim=n)
  counter=0
  while(counter<=n){
    candidate <- runif(1,0,100)
    if(runif(1,0,1)<dnegexp(candidate,thisa, thisb)) {
      rnums[counter] <- candidate
      counter <- counter+1
    }
  }
  return(rnums)
}

Имея в руках эти две функции, я тестирую свой рабочий процесс, создавая набор данных, который следует отрицательному экспоненциальному распределению с известными параметрами:

set.seed(501)
mynegexpvals <- rnegexp(100,0.08, 0.09)

hist(mynegexpvals, freq=FALSE)
mydists <- seq(0,100, by=1)
lines(mydists,dnegexp(mydists, 0.08, 0.09), col="blue", lwd=2)

гистограмма сгенерированных значений со строкой генерации pdf

Когда я пытаюсь использовать mle2 из пакета bbmle для поиска параметров, он дает бессмысленные значения, хотя я дал ему точные параметры генерации в качестве начальных значений:

> library(bbmle)
> mle2(mynegexpvals ~ dnegexp(a,b), start = list(a=0.08, b=0.09), data=data.frame(mynegexpvals))

Call:
mle2(minuslogl = mynegexpvals ~ dnegexp(a, b), start = list(a = 0.08, 
    b = 0.09), data = data.frame(mynegexpvals))

Coefficients:
            a             b 
 2.421577e+12 -6.849330e+12 

Log-likelihood: 6.148807e+15

Попытки ограничить пространство поиска дают оценки параметров на границе:

> mle2(mynegexpvals ~ dnegexp(a,b), start = list(a=0.08, b=0.09), data=data.frame(mynegexpvals), method="L-BFGS-B", lower=c(a=0.04, b=0.0001), upper=c(a=1000, b=1000))

Call:
mle2(minuslogl = mynegexpvals ~ dnegexp(a, b), start = list(a = 0.08, 
    b = 0.09), method = "L-BFGS-B", data = data.frame(mynegexpvals), 
    lower = c(a = 0.04, b = 1e-04), upper = c(a = 1000, b = 1000))

Coefficients:
    a     b 
1e+03 1e-04 

Log-likelihood: 690.69 
Warning message:
In mle2(mynegexpvals ~ dnegexp(a, b), start = list(a = 0.08, b = 0.09),  :
  some parameters are on the boundary: variance-covariance calculations based on Hessian may be unreliable

Я пропустил что-то важное здесь? Я неправильно определил свою функцию dnegexp?


person emudrak    schedule 16.02.2016    source источник
comment
Я думаю, что фундаментальная проблема здесь заключается в том, что ваша функция dnegexp не нормализована (т.е. интеграл от 0 до бесконечности dnegexp(x,...) не всегда равен 1). a*exp(-b*x), данное в Баллоке et al., представляет собой плотность семян на расстоянии x, а не плотность вероятности того, что семена преодолеют расстояние x . .. Вы должны выяснить, соответствуют ли ваши данные исходному распределению, исходной тени или дисперсионному распределению (все немного отличается).   -  person Ben Bolker    schedule 18.02.2016