What is an appropriate method for providing bounds when performing maximum likelihood parameter estimation?

Question

I have successfully implemented a maximum likelihood estimation of model parameters with bounds by creating a likelihood function that returns NA or Inf values when the function is out of bounds. I optimize the function using optim in R.

detailed example available on github

Quick example:

likelihood.fun<-function(par, ...){

  likelihood<- -sum(dnorm(..., log=T))

  if(any(c(par[1]<0,
         par[1]>5, 
         par[2]>5,
         par[2]<0)){likelihood<-NA}

  return(likelihood)

}

Is this equivalent to box optimization or deprecated compare to box optimization?

If this is not equivalent:

How can I implement this using

optim(..., method="L-BFGS-B", lower=c(...), upper=c(...))


from example, this does not seem to work:
optim(..., method="L-BFGS-B", lower=c(0,0), upper=c(5,5))

or

constrOptim()

This is linked to this question on constrOptim.

Answer 1

What you are doing in your first code block is indeed equivalent to box-constrained optimisation. Here's some sample code, with some unnecessary output removed to save space:

> foo.unconstr <- function(par, x) -sum(dnorm(x, par[1], par[2], log=TRUE))
> 
> foo.constr <- function(par, x)
+ {
+   ll <- NA
+   if (par[1] > 0 && par[1] < 5 && par[2] > 0 && par[2] < 5)
+   {
+     ll <- -sum(dnorm(x, par[1], par[2], log=TRUE))
+   }
+   ll
+ }
> 
> x <- rnorm(100,1,1)
> par <- c(1,1)
> optim(par, foo.constr, x=x)
$par
[1] 1.147690 1.077712

$value
[1] 149.3724

> 
> par <- c(1,1)
> optim(par, foo.unconstr, lower=c(0,0), upper=c(5,5), method="L-BFGS-B", x=x)
$par
[1] 1.147652 1.077654

$value
[1] 149.3724

They won't give quite the same answers, because they are different algorithms.

I'll answer your constrOptim question over there, so other people who might be interested will see it.