Минимизация вероятности отрицательного журнала: ошибка в optim (), вызванная начальными значениями

Я хочу получить параметры максимального правдоподобия (MLE) для кумулятивной нормальной кривой, подогнанной к некоторым данным о пропорциях, путем прямой минимизации отрицательного логарифмического правдоподобия (без использования glm). Для некоторых начальных значений, введенных в optim (), проблем нет:

x <- c(-0.250, -0.056,  0.137,  0.331,  0.525,  0.719,  0.912,  1.100,  1.300)
k <- c(0, 0, 5, 11, 12, 12, 12, 12, 12)
n <- c(12, 12, 12, 12, 12, 12, 12, 12, 12)


nll <- function(p) {
  phi <- pnorm(x, p[1], p[2])
  -sum(k * log(phi) + (n - k) * log(1 - phi))
}

para<- optim(c(0.5, 0.1), nll)$par

xseq <- seq(-.5, 1.5, len = 100)
yseq <- pnorm(xseq, para[1],para[2])
curve <- data.frame(xseq, yseq)

dat <- data.frame(x, k, n)

library(ggplot2)
ggplot(dat,aes(x = x, y = k / n)) + 
  geom_point()+
  geom_line(data = curve, aes(x = xseq, y = yseq))

Но, если я использую начальные значения, которые на самом деле ближе к параметрам MLE

 para<- optim(c(0.1, 0.1), nll)$par

Я получил следующую ошибку:

Error in optim(c(0.1, 0.1), nll) : function cannot be evaluated at initial parameters

Кажется, что ошибка вызвана некоторой бесконечностью в оценке отрицательного логарифмического правдоподобия. Я обнаружил, что если я увеличиваю точность, используя параметр log.p параметра pnorm, я не обнаруживаю ошибку.

nll <- function(p) {
  logphi1 <- pnorm(x, p[1], p[2], lower.tail = T, log.p = T)
  logphi2 <- pnorm(x, p[1], p[2], lower.tail = F, log.p = T)
  -sum(k * logphi1 + (n - k) * logphi2)
}
para<- optim(c(0.1, 0.1), nll)$par

но проблема в том, что в дополнение к pnorm я также хочу подогнать a + b * pnorm кривые с константами a и b, и в этих случаях я не могу использовать log.p для увеличения точности.


r optimization mle
person danilinares    schedule 11.12.2014    source источник

Ответы (2)


arrow_upward
1
arrow_downward

Кажется, что замена очень маленьких чисел машинным эпсилон и чисел, очень близких к 1, на 1 - (машинный эпсилон), ошибки не возникает, и соответствие кажется разумным.

x <- c(-0.250, -0.056,  0.137,  0.331,  0.525,  0.719,  0.912,  1.100,  1.300)
k <- c(0, 0, 5, 11, 12, 12, 12, 12, 12)
n <- c(12, 12, 12, 12, 12, 12, 12, 12, 12)

nll <- function(p) {
  phi <- pnorm(x, p[1], p[2])
  phi[phi < .Machine$double.eps] <- .Machine$double.eps
  phi[phi > (1 - .Machine$double.eps)] <- 1 - .Machine$double.eps
  -sum(k * log(phi) + (n - k) * log(1 - phi))
}

para<- optim(c(0.1, 0.1), nll)$par

xseq <- seq(-.5, 1.5, len = 100)
yseq <- pnorm(xseq, para[1],para[2])
curve <- data.frame(xseq, yseq)

dat <- data.frame(x, k, n)

library(ggplot2)
ggplot(dat,aes(x = x, y = k / n)) +
  geom_point()+
  geom_line(data = curve, aes(x = xseq, y = yseq))
person danilinares    schedule 12.12.2014

arrow_upward
0
arrow_downward

Проблема заключается в 8-й точке данных и значениях параметров; они вызывают NaN при оценке вероятности, потому что pnorm оценивается как 1 (численно):

p <- c(0.1,0.1)
pnorm(x[8], p[1], p[2])
## 1
1-pnorm(x[8], p[1], p[2])
## 0
pnorm(x[8], p[1], p[2], lower.tail=FALSE)
## 7.6e-24

Последнее значение находится ниже машинного эпсилон, поэтому, даже если вы с большой вероятностью напишете 1 - pnorm(x[8], p[1], p[2], lower.tail=FALSE), это не предотвратит потери значимости.

person renato vitolo    schedule 11.12.2014