R help - Sep 2015 - maximum-likelihood-estimation with mle()

Hi everyone,

I have a problem with maximum-likelihood-estimation in the following
situation:

Assume a functional relation y = f(x) (the specific form of f should be
irrelevant). For my observations I assume (for simplicity) white noise,
such that hat(y_i) = f(x_i) + epsilon_i, with the epsilon_i iid N(0,
sigma^2). Y_i should then be N(f(x_i), sigma^2)-distributed and due to
the iid assumption the density of Y = (Y_1, ..., Y_n) is simply the
product of the individual densities, taking the log gives the the sum of
the log of individual densities.

I tried coding this in R with a simple example: f(x) = a*x + b (simple
linear regression). This way I wanted to compare the results from my
ml-estimation (specifying the log-likelihood manually and estimating
with mle()) with the results from using lm(y~x). In my example however
it doesn't work:

x <- 1:10
y <- 3*x - 1 + rnorm(length(x), mean=0, sd=0.5)

library("stats4")
nLL <- function(a, b, sigma)
{
  -sum(dnorm(y, mean=a*x+b, sd=sigma, log=TRUE))
}

fit <- mle(nLL1, start=list(a=0, b=0, sigma=1), nobs=length(y))

summary(lm(y~x))
summary(fit)

These should be the same but the aren't. I must have made some mistake
specifying the (negative) log-likehood (but I just don't see it). I also
actually don't care much (at the moment) for estimating sigma but I
don't know of a way to specify (and estimate) the (negative)
log-likelihood without estimating sigma.

Thanks a lot
and kind regards

Ronald Koelpin

You are being over-optimistic with your starting values, and/or with constrains
on the parameter space.
Your fit is diverging in sigma for some reason known only to nonlinear-optimizer
gurus...

For me, it works either to put in an explicit constraint or to reparametrize
with log(sigma).

E.g.
> fit <- mle(nLL, start=list(a=0, b=0, sigma=1),
method="L-BFGS-B", lower=c(-Inf, -Inf, 0.001))
> summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = nLL, start = list(a = 0, b = 0, sigma = 1), method =
"L-BFGS-B",
    lower = c(-Inf, -Inf, 0.001))

Coefficients:
        Estimate Std. Error
a      3.0293645 0.01904026
b     -1.1516406 0.11814174
sigma  0.1729418 0.03866673

-2 log L: -6.717779 

or 
> nLL <- function(a, b, lsigma)
+  -sum(dnorm(y, mean=a*x+b, sd=exp(lsigma), log=TRUE))
> fit <- mle(nLL, start=list(a=0, b=0, lsigma=0))
> summary(fit)
Maximum likelihood estimation

Call:
mle(minuslogl = nLL, start = list(a = 0, b = 0, lsigma = 0))

Coefficients:
        Estimate Std. Error
a       3.029365 0.01903975
b      -1.151642 0.11813856
lsigma -1.754827 0.22360675

-2 log L: -6.717779 

both of which reproduce lm() except for DF issues.

To fix sigma, use fixed=list(sigma=0.19) in mle(). This also stabilizes the
convergence (which it blinking well should, since it is now a purely quadratic
problem).

-pd


On 11 Sep 2015, at 14:20 , Ronald K?lpin <ronald.koelpin at gmail.com>
wrote:
> Hi everyone,
> 
> I have a problem with maximum-likelihood-estimation in the following
> situation:
> 
> Assume a functional relation y = f(x) (the specific form of f should be
> irrelevant). For my observations I assume (for simplicity) white noise,
> such that hat(y_i) = f(x_i) + epsilon_i, with the epsilon_i iid N(0,
> sigma^2). Y_i should then be N(f(x_i), sigma^2)-distributed and due to
> the iid assumption the density of Y = (Y_1, ..., Y_n) is simply the
> product of the individual densities, taking the log gives the the sum of
> the log of individual densities.
> 
> I tried coding this in R with a simple example: f(x) = a*x + b (simple
> linear regression). This way I wanted to compare the results from my
> ml-estimation (specifying the log-likelihood manually and estimating
> with mle()) with the results from using lm(y~x). In my example however
> it doesn't work:
> 
> x <- 1:10
> y <- 3*x - 1 + rnorm(length(x), mean=0, sd=0.5)
> 
> library("stats4")
> nLL <- function(a, b, sigma)
> {
>  -sum(dnorm(y, mean=a*x+b, sd=sigma, log=TRUE))
> }
> 
> fit <- mle(nLL1, start=list(a=0, b=0, sigma=1), nobs=length(y))
> 
> summary(lm(y~x))
> summary(fit)
> 
> These should be the same but the aren't. I must have made some mistake
> specifying the (negative) log-likehood (but I just don't see it). I
also
> actually don't care much (at the moment) for estimating sigma but I
> don't know of a way to specify (and estimate) the (negative)
> log-likelihood without estimating sigma.
> 
> Thanks a lot
> and kind regards
> 
> Ronald Koelpin
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com

peter dalgaard <pdalgd <at> gmail.com> writes:
> 
> You are being over-optimistic with your starting values, and/or 
> with constrains on the parameter space. 
> Your fit is diverging in sigma for some reason known
>  only to nonlinear-optimizer gurus...
> 
> For me, it works either to put in an explicit
>  constraint or to reparametrize with log(sigma).
  [snip]

  Another few strategies:

set.seed(101)
x <- 1:10
y <- 3*x - 1 + rnorm(length(x), mean=0, sd=0.5)

library("stats4")
nLL <- function(a, b, sigma) {
  -sum(dnorm(y, mean=a*x+b, sd=sigma, log=TRUE))
}

fit <- mle(nLL, start=list(a=0, b=0, sigma=1), nobs=length(y))


   Nelder-Mead is generally more robust than BFGS, although slower:


fit2 <- mle(nLL, start=list(a=0, b=0, sigma=1),
     method="Nelder-Mead",nobs=length(y))

   bbmle can be used to simplify things slightly (this repeats
Peter Dalgaard's transformation of sigma):

library("bbmle")
fit3 <- mle2(y~dnorm(mean=mu,sd=exp(logsigma)),
   parameters=list(mu~x),
   data=data.frame(x,y),
   start=list(mu=0, logsigma=0))

  You can also profile out the sigma:

## dnorm with sd profiled out
dnorm2 <- function(x,mean,log=FALSE) {
  ssq <- sum((x-mean)^2)
  dnorm(x,mean,sd=sqrt(ssq/length(x)),log=log)
}
fit4 <- mle2(y~dnorm2(mean=mu),
   parameters=list(mu~x),
   data=data.frame(x,y),
   start=list(mu=0))