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().

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

--
Rolf Turner
Technical Editor ANZJS

______________________________________________
R-help@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.

Reply via email to