[R] MLE with parameters restricted to a defined range using bbmle

Ben Bolker bbolker at gmail.com
Tue Dec 9 02:05:39 CET 2014

```Rolf Turner <r.turner <at> auckland.ac.nz> writes:

>
>
> I know nothing about the bbmle package and its mle2() function, but it
> is a general truth that if you need to constrain a parameter to be
> positive in an optimisation procedure a simple and effective approach is
> to reparameterize using exp().

mle2() is a wrapper for R's built-in optim() function, so
you can use method="L-BFGS-B" and set the minimum values via
the lower= argument.  The only potentially tricky part is that
you may want to set the lower bound slightly above the desired
bound, as L-BFGS-B uses as finite difference approximation to
compute the gradients, and I'm not 100% sure that the finite
difference computation always respects the bounds automatically.
(The finite-difference stepsize is set by the 'ndeps' parameter
and is 0.001 by default.)

>
> I.e. represent xmin as exp(lxmin) (say) and use lxmin as the argument
>
> This strategy rarely if ever fails to work.
>
> cheers,
>
> Rolf Turner
>
> On 09/12/14 09:04, Bernardo Santos wrote:
> > Dear all,
> > I am fitting models to data with mle2 function of the bbmle package.
> In specific, I want to fit a power-law
> distribution model, as defined here
> (http://arxiv.org/pdf/cond-mat/0412004v3.pdf), to data.
> > However, one of the parameters - xmin -, must be necessarily greater
> than zero. What can I do to restrict the
> possible values of a parameter that are passed to the optimizer?
> > Here there is a sample of my code:
library(bbmle)

# Creating data
set.seed(1234)
data <- rexp(1000, rate = 0.1) #
## The fit will not be too good, but it is just to test

# Creating the power-law distribution density function
dpowlaw <- function(x, alfa, xmin, log=FALSE){
c <- (alfa-1)*xmin^(alfa-1)
if(log) ifelse(x < xmin, 0, log(c*x^(-alfa)))
else ifelse(x < xmin, 0, c*x^(-alfa))
}
# Testing the function
integrate(dpowlaw, -Inf, Inf, alfa=2, xmin=1)
curve(dpowlaw(x, alfa=2.5, xmin=10), from=0, to=100, log="")
curve(dpowlaw(x, alfa=2.5, xmin=1), from=1, to=100, log="xy")

# Negative log-likelihood function
LLlevy <- function(mu, xmin){
-sum(dpowlaw(data, alfa=mu, xmin=xmin, log=T))
}

# Fitting model to data
mlevy <- mle2(LLlevy, start=list(mu=2, xmin=1))
The result of model fitting here is Coefficients:
mu      xmin
-916.4043  890.4248

but this does not make sense!xmin must be > 0,
and mu must be > 1.What should I do?