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 ______________________________________________ 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.
