2

我想通过直接最小化负对数似然（不使用 glm）来获得拟合到某些比例数据的累积正态曲线的最大似然参数（MLE）。对于引入 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.ppnorm 选项提高精度，我不会得到错误

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提高精度。

4

2 回答 2

1

似乎用机器 epsilon 替换非常小的数字和非常接近 1 的数字 1 - (机器 epsilon)，不会发生错误并且拟合似乎是明智的。

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))
于 2014-12-12T14:10:16.087 回答
0

问题出在第 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

后一个值低于机器epsilon，因此即使您写出1 - pnorm(x[8], p[1], p[2], lower.tail=FALSE)您的可能性，这也无法避免下溢。

于 2014-12-11T22:12:02.900 回答