Bernardo Santos
2014-Dec-08 20:04 UTC
[R] MLE with parameters restricted to a defined range using bbmle
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]]
Rolf Turner
2014-Dec-08 20:35 UTC
[R] MLE with parameters restricted to a defined range using bbmle
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
Ben Bolker
2014-Dec-09 01:05 UTC
[R] MLE with parameters restricted to a defined range using bbmle
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
Bernardo Santos
2014-Dec-10 13:35 UTC
[R] MLE with parameters restricted to a defined range using bbmle
Thanks, Rolf and Ben, Both solutions worked, but I finished by contraining parameters values using optim(). Best,Bernardo Niebuhr Em Segunda-feira, 8 de Dezembro de 2014 18:36, Rolf Turner <r.turner at auckland.ac.nz> escreveu: 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 [[alternative HTML version deleted]]