Poor performance of "Optim"

Ben Bolker sent me a private email rightfully correcting me that was
factually wrong when I wrote that ML /is/ a numerical method (I had
written sloppily and under time pressure). He is of course right to point
out that not all maximum likelihood estimators require numerical methods
to solve. Further, only numerical optimization will show the behavior
discussed in this post for the given reasons. (I hope this post isn't yet
another blooper of mine at 5 a.m. in the morning). 

Best,
Daniel


Daniel Malter wrote:

With respect, your statement that R's optim does not give you a reliable
estimator is bogus. As pointed out before, this would depend on when
optim believes it's good enough and stops optimizing. In particular if
you stretch out x, then it is plausible that the likelihood function will
become flat enough "earlier," so that the numerical optimization will
stop earlier (i.e., optim will "think" that the slope of the likelihood
function is flat enough to be considered zero and stop earlier than it
will for more condensed data). After all, maximum likelihood is a
numerical method and thus an approximation. I would venture to say that
what you describe lies in the nature of this method. You could also
follow the good advice given earlier, by increasing the number of
iterations or decreasing the tolerance. 

However, check the example below: for all purposes it's really close
enough and has nothing to do with optim being "unreliable."

n<-1000
x<-rnorm(n)
y<-0.5*x+rnorm(n)
z<-ifelse(y>0,1,0)

X<-cbind(1,x)
b<-matrix(c(0,0),nrow=2)

#Probit
reg<-glm(z~x,family=binomial("probit"))

#Optim reproducing probit (with minor deviations due to difference in
method)
LL<-function(b){-sum(z*log(pnorm(X%*%b))+(1-z)*log(1-pnorm(X%*%b)))}
optim(c(0,0),LL)

#Multiply x by 2 and repeat optim
X[,2]=2*X[,2]
optim(c(0,0),LL)

HTH,
Daniel



yehengxin wrote:

What I tried is just a simple binary probit model.   Create a random
data and use "optim" to maximize the log-likelihood function to estimate
the coefficients.   (e.g. u = 0.1+0.2*x + e, e is standard normal.  And
y = (u > 0),  y indicating a binary choice variable)

If I estimate coefficient of "x", I should be able to get a value close
to 0.2 if sample is large enough.  Say I got 0.18. 

If I expand x by twice and reestimate the model, which coefficient
should I get?  0.09, right?

But with "optim", I got something different.  When I do the same thing
in both Gauss and Matlab, I can exactly get 0.09, evidencing that the
coefficient estimator is reliable.  But R's "optim" does not give me a
reliable estimator.