[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
> 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
More information about the R-help
mailing list