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
[[alternative HTML version deleted]]
______________________________________________
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.
