8
$\begingroup$

I have a variable set of responses that are expressed as an interval such as the sample below.

> head(left)
[1]  860  516  430 1118  860  602
> head(right)
[1]  946  602  516 1204  946  688

where left is the lower bound and right is the upper bound of the response. I want to estimate the parameters according to the lognormal distribution.

For a while when I was trying to calculate the likelihoods directly I was struggling with the fact that since the two bounds are distributed along different set of paramaters, I was getting some negative values like below:

> Pr_high=plnorm(wta_high,meanlog_high,sdlog_high)
> Pr_low=plnorm(wta_low, meanlog_low,sdlog_low)
> Pr=Pr_high-Pr_low
> 
> head(Pr)
[1] -0.0079951419  0.0001207749  0.0008002343 -0.0009705125 -0.0079951419 -0.0022395514

I couldn't really figure out how to resolve it and decided to use the mid-point of the interval instead which is a good compromise until I found mledist function which extracts the loglikelihood of an interval response, this is the summary I get:

> mledist(int, distr="lnorm")
$estimate
meanlog     sdlog 
6.9092257 0.3120138 

$convergence
[1] 0

$loglik
[1] -152.1236

$hessian
         meanlog       sdlog
meanlog 570.760358    7.183723
sdlog     7.183723 1112.098031

$optim.function
[1] "optim"

$fix.arg
NULL

Warning messages:
1: In plnorm(q = c(946L, 602L, 516L, 1204L, 946L, 688L, 1376L, 1376L,  :
NaNs produced
2: In plnorm(q = c(860L, 516L, 430L, 1118L, 860L, 602L, 1290L, 1290L,  :
NaNs produced

The parameter values seem to make sense and the loglikelihood is greater than any other method I have used (mid-point distribution or distribution of either one of the bounds).

There is a warning message which I don't understand so could anyone tell me if I am doing the right thing and what this message means?

Appreciate the help!

$\endgroup$
  • $\begingroup$ Your question amounts to "How do I use a particular R function and what does this Warning message mean?". That's a question for StackOverflow rather than CrossValidated. Further, when you refer to a function from a package, you should mention what package it's from. In this case I presume you mean the function from the package fitdistrplus. $\endgroup$ – Glen_b Apr 13 '13 at 22:19
  • $\begingroup$ Welcome to the site, @ElioDruml. I can't tell if your main question is about how to estimate these parameters, or what the meaning of the warning message is. The former would be a good question for CV, but the latter is really a question for Stack Overflow (see our FAQ). Can you clarify what your primary question is? Would you prefer your Q stay here, or be migrated to SO? (If the latter, flag your Q & we'll migrate it for you, please don't cross-post, though.) $\endgroup$ – gung - Reinstate Monica Apr 13 '13 at 23:16
9
$\begingroup$

It sounds like you might not be computing the likelihood correctly.

When all you know about a value $x$ is that

  1. It is obtained independently from a distribution $F_\theta$ and

  2. It lies between $a$ and $b \gt a$ inclusive (where $b$ and $a$ are independent of $x$),

then (by definition) its likelihood is $${\Pr}_{F_\theta}(a \le x \le b) = F_\theta(b) - F_\theta(a).$$ The likelihood of a set of independent observations therefore is the product of such expressions, one per observation. The log-likelihood, as usual, will be the sum of logarithms of those expressions.


As an example, here is an R implementation where the values of $a$ are in the vector left, the values of $b$ in the vector right, and $F_\theta$ is Lognormal. (This is not a general-purpose solution; in particular, it assumes that $b \gt a$ and $b \ne a$ for all the data.)

#
# Lognormal log-likelihood for interval data.
#
lambda <- function(mu, sigma, left, right) {
  sum(log(pnorm(log(right), mu, sigma) - pnorm(log(left), mu, sigma)))
}

To find the maximum log likelihood we need a reasonable set of starting values for the log mean $\mu$ and log standard deviation $\sigma$. This estimate replaces each interval by the geometric mean of its endpoints:

#
# Create an initial estimate of lognormal parameters for interval data.
#
lambda.init <- function(left, right) {
  mid <- log(left * right)/2
  c(mean(mid), sd(mid))
}

Let's generate some random lognormally distributed data and bin them into intervals:

set.seed(17)
n <- 12                     # Number of data
z <- exp(rnorm(n, 6, .5))   # Mean = 6, SD = 0.5
left <- 100 * floor(z/100)  # Bin into multiples of 100
right <- left + 100

The fitting can be performed by a general-purpose multivariate optimizer. (This one is a minimizer by default, so it must be applied to the negative of the log-likelihood.)

fit <- optim(lambda.init(left,right), 
             fn=function(theta) -lambda(theta[1], theta[2], left, right))
fit$par

6.1188785 0.3957045

The estimate of $\mu$ is $6.12$, not far from the intended value of $6$, and the estimate of $\sigma$ is $0.40$, not far from the intended value of $0.5$: not bad for just $12$ values. To see how good the fit is, let's plot the empirical cumulative distribution function and the fitted distribution function. To construct the ECDF, I just interpolate linearly through each interval:

#
# ECDF of the data.
#
F <- function(x) (1 + mean((abs(x - left) - abs(x - right)) / (right - left)))/2

y <- sapply(x <- seq(min(left) * 0.8, max(right) / 0.8, 1), F)
plot(x, y, type="l", lwd=2, lty=2, ylab="Cumulative probability")
curve(pnorm(log(x), fit$par[1], fit$par[2]), from=min(x), to=max(x), col="Red", lwd=2, 
  add=TRUE)

Plots

Because the vertical deviations are consistently small and vary both up and down, it looks like a good fit.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks a lot for your input @whuber. I did recreate your example and it all makes sense. However, I wasn't able to recreate on my own data of n=56 of which the head is left <- c(860, 516, 430, 1118, 860, 602) and right <- c(946, 602, 516, 1204, 946, 688). I get this warning message: "1: In pnorm(log(right), mu, sigma) : NaNs produced 2: In pnorm(log(left), mu, sigma) : NaNs produced" when fitting with the optimizer to extract the mle estimates. That brings me back to my earlier problem of having negative probabilities when calc. the likelihoods step by step and substracting. $\endgroup$ – Elio Druml Apr 16 '13 at 11:46
  • $\begingroup$ These are the same warning messages given by the mledist function from the fitdistrplus package. However, as you can see above it does give me an output for the mle estimates that look relatively good. Should I trust it and/or what is the issue here? Thanks for the feedback. $\endgroup$ – Elio Druml Apr 16 '13 at 12:01
  • $\begingroup$ Why don't you post your data, Elio, so we can diagnose the problem? Even so, I'm not sure these are critical errors. You might be experiencing the same troubles reported by another user when numerically minimizing a function in Mathematica; the same explanation might apply in your case. $\endgroup$ – whuber Apr 16 '13 at 14:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.