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

[Ben Bolker]

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
> to your objective function.
> 
> 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:
# Loading library
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?
Thanks in advance!Bernardo Niebuhr