来自算法的奇怪错误

Question

我编写了这个函数,它根据Newton-Raphson算法以数字方式计算Cauchy分布中的MLE:

mlec <- function(x,theta0=median(x),numstp=100,eps=0.01){
    numfin <- numstp
    ic <- 0
    istop <- 0
    while(istop==0){
        ic <- ic+1
        ltheta <- -2*sum((x-theta0)/(1+(x-theta0)^2))
        lprimetheta <- -2*(sum(2*(x-theta0)^2/
                            (1+(x-theta0)^2)^2-1/(1+(x-theta0)^2)^2))
        theta1 <- theta0-(ltheta/lprimetheta)
        check <- abs((theta1-theta0)/theta1)
        if(check < eps ) { istop <- 1 }
        theta0 <- theta1
    }
    list(theta1=theta1,check=check,realnumstps=ic)
}

然后,目标是从具有比例参数2的Cauchy分布生成观测值,并查看MLE如何执行.问题在于,对于某些样本,MLE对其他人运行得非常好,我得到了奇怪的错误

Error in if (check < eps) { : missing value where TRUE/FALSE needed

这里发生了什么？我已经定义了什么"检查"是不应该发生的.

谢谢.

Answer 1

我添加了一些仪器(参见cat()中间的语句),并修复了二阶导数表达式(fixed见下文):

mlec <- function(x,theta0=median(x),numstp=100,eps=0.01,
                 debug=TRUE,fixed=FALSE){
    numfin <- numstp
    ic <- 0
    istop <- 0
    while(istop==0){
        ic <- ic+1
        ltheta <- -2*sum((x-theta0)/(1+(x-theta0)^2))
        lprimetheta <- -2*(sum(2*(x-theta0)^2/
                            (1+(x-theta0)^2)^2-1/(1+(x-theta0)^2)^2))
        if (!fixed) {
            theta1 <- theta0-(ltheta/lprimetheta)
        } else theta1 <- theta0-ltheta/ff(theta0)
        check <- abs((theta1-theta0)/theta1)
        if (debug) cat(ic,ltheta,lprimetheta,theta0,theta1,check,"\n")
        if(check < eps ) { istop <- 1 }
        theta0 <- theta1
    }
    list(theta1=theta1,check=check,realnumstps=ic)
}

set.seed(1)
x <- rcauchy(100,2)

mlec(x)

这是输出的尾端:

## ic  ltheta        lprimetheta      theta0        theta1     check
## 427 -4.48838e-75 -2.014555e-151 -4.455951e+76 -6.683926e+76 0.3333333 
## 428 -2.992253e-75 -8.953579e-152 -6.683926e+76 -1.002589e+77 0.3333333 
## 429 -1.994835e-75 -3.979368e-152 -1.002589e+77 -1.503883e+77 0.3333333 
## 430 -1.32989e-75 0 -1.503883e+77 -Inf NaN 

现在,它为什么会发生？要么你某处有错误,要么算法不稳定. 文艺青年最爱的事实证明,答案居然是"既"; 你的二阶导数表达似乎是错误的,但即使它是正确的,NR算法对于这个问题极其不稳定.

这是你的衍生函数和二阶导数函数(Vectorize()为了方便起见,我将它们包装起来,以便我可以theta同时评估多个值):

lthetafun <- Vectorize(function(theta) {
    -2*sum((x-theta)/(1+(x-theta)^2))
})
lprimethetafun <- Vectorize(function(theta) {
    -2*(sum(2*(x-theta)^2/
                (1+(x-theta)^2)^2-1/(1+(x-theta)^2)^2))
})

基于内置函数的负对数似然dcauchy函数:

thetafun <- Vectorize(function(theta) -sum(dcauchy(x,theta,log=TRUE)))

检查区别(在任意选择的点):

library("numDeriv")
all.equal(grad(thetafun,2),lthetafun(2))  ## TRUE

检查二阶导数:

hessian(thetafun,2) ## 36.13297
lprimethetafun(2)   ## 8.609859: ???

我认为你的二阶导数表达式是错误的.

以下替代二阶导数函数基于与Wolfram Alpha的懒惰作弊,区分您的渐变函数(与有限差分近似匹配):

ff <- Vectorize(function(theta)
    2*sum(((x-theta)^2-1)/((x-theta)^2+1)^2))
ff(2)  ## matches hessian() above.

但看起来你可能还有其他问题.

负对数似然表面看起来正常:

curve(thetafun, from=-10,to=10,n=501)

但麻烦即将出现:

curve(lthetafun, from=-10,to=10, n=501)

这看起来不规则,向上一级到二阶导数表明它是:

curve(ff, from=-10, to=10, n=501)

以下是NR更新的曲线:

ff2 <- function(x) x-lthetafun(x)/ff(x)
curve(ff2, from=-10, to=10, n=501,ylim=c(-100,100))

哎呀!这表明为什么Newton-Raphson方法可能出错,除非你开始接近最小值(任何时候可能性表面都有一个拐点,NR更新曲线将有一个极点......).对问题的进一步分析可能会告诉你为什么Cauchy的二阶导数是如此可怕.

如果你只是想找到MLE,你可以通过一些更强大的1-D方法来实现:

library("bbmle")
mle2(x~dcauchy(location=m),
     data=data.frame(x),
     start=list(m=median(x)),
     method="Brent",
     lower=-100,upper=100)
## 
## Call:
## mle2(minuslogl = x ~ dcauchy(location = m), start = list(m = median(x)), 
##     method = "Brent", data = data.frame(x), lower = -100, upper = 100)
## 
## Coefficients:
##       m 
## 1.90179 
## 
## Log-likelihood: -262.96 
## 

如果你开始足够近,NR似乎工作正常:

  mlec(x,1.85,debug=FALSE,fixed=TRUE,eps=0.0001)
 ## $theta1
 ## [1] 1.901592
 ## 
 ## $check
 ## [1] 5.214763e-05
 ## 
 ## $realnumstps
 ## [1] 37079