Search code examples
rfunctionstatisticsmathematical-optimizationmle

R Optim() function: Truncated Log-normal Maximum Likelihood Estimation (solve for mu and sd)

I'm attempting to fit a truncated log-normal distribution in R. In the below x is list of integers. I believe the function Optim() is likely to be the best way to do this, as below. 

log.lklh.lnorm <- function(x, mu, sd, input_min, input_max){-sum(log(dlnorm(x, meanlog = mu, sdlog = sd, log = FALSE)/((plnorm(input_max, meanlog = mu, sdlog = sd, lower.tail = TRUE, log.p = FALSE)) - (plnorm(input_min, meanlog = mu, sdlog = sd, lower.tail = TRUE, log.p = FALSE)))))}


optim(par = c(0,1,1,100000), log.lklh.lnorm)

I used excel solver to solve (i.e. find the maximum value of this sum) for mu and sd (subject to the input max and min). However, i cant seem to replicate this in R. I've tried various versions of the above code, including:

  • Set input min and input max to certain values, such as 1 and 100,000 respectively.
  • Order of Optim() calls, i.e. input x, par and fn order)

Any help is greatly appreciated, thanks!

Dave

Solution

  • I think you want to optimize a function of just two parameters. So, we could do:

    x <- c(1,2,3,3,3,4,4,5,5,6,7)

log.lklh.lnorm <- function(x, mu, sd, input_min, input_max){-sum(log(dlnorm(x, meanlog = mu, sdlog = sd, log = FALSE)/((plnorm(input_max, meanlog = mu, sdlog = sd, lower.tail = TRUE, log.p = FALSE)) - (plnorm(input_min, meanlog = mu, sdlog = sd, lower.tail = TRUE, log.p = FALSE)))))}


f <- function(y) {
  log.lklh.lnorm(x,y[1],y[2],1,100000)    
}

optim(par = c(3,1), f)

    Notes in response to further questions in the comments:

    1. What is the function of par = c(3,1)? This has two functions. (1) It is the initial point (guess). (2) It determines the number of parameters to optimize for. In this case, this is 2. It is noted that you can (and should) calculate very good initial values. You can get very good estimates for μ and σ just by calculating the average and sample standard deviation. (In statistics this is sometimes called method of moments) This will make the task of optim much easier. So instead of par=c(3,1) we really should do: par = c(mean(x),sd(x)).

    2. What does y[1]/y[2] do?. The optim function organizes all parameters to optimize for, as a single vector. So I represented (μ,σ) by a single vector y (of length 2). In the function f, this array is unpacked and passed on to log.lklh.lnorm as individual arguments.