Need algorithm to compute relative likelihood that data are sample from normal vs lognormal distribution

Question

Let's say you have a set of values, and you want to know if it is more likely that they were sampled from a Gaussian (normal) distribution or sampled from a lognormal distribution?

Of course, ideally you'd know something about the population or about the sources of experimental error, so would have additional information useful to answering the question. But here, assume we only have a set of numbers and no other information. Which is more likely: sampling from a Gaussian or sampling from a lognormal distribution? How much more likely? What I am hoping for is an algorithm to select between the two models, and hopefully quantify the relative likelihood of each.

Answer 1

You could take a best guess at the distribution type by fitting each distribution (normal or lognormal) to the data by maximum likelihood, then comparing the log-likelihood under each model - the model with the highest log-likelihood being the best fit. For example, in R:

# log likelihood of the data given the parameters (par) for 
# a normal or lognormal distribution
logl <- function(par, x, lognorm=F) {
    if(par[2]<0) { return(-Inf) }
    ifelse(lognorm,
    sum(dlnorm(x,par[1],par[2],log=T)),
    sum(dnorm(x,par[1],par[2],log=T))
    )
}

# estimate parameters of distribution of x by ML 
ml <- function(par, x, ...) {
    optim(par, logl, control=list(fnscale=-1), x=x, ...)
}

# best guess for distribution-type
# use mean,sd of x for starting parameters in ML fit of normal
# use mean,sd of log(x) for starting parameters in ML fit of lognormal
# return name of distribution type with highest log ML
best <- function(x) {
    logl_norm <- ml(c(mean(x), sd(x)), x)$value
        logl_lognorm <- ml(c(mean(log(x)), sd(log(x))), x, lognorm=T)$value
    c("Normal","Lognormal")[which.max(c(logl_norm, logl_lognorm))]
}

Now generate numbers from a normal distribution and fit a normal distribution by ML:

set.seed(1)
x = rnorm(100, 10, 2)
ml(c(10,2), x)

Produces:

$par
[1] 10.218083  1.787379

$value
[1] -199.9697
...

Compare log-likelihood for ML fit of normal and lognormal distributions:

ml(c(10,2), x)$value # -199.9697
    ml(c(2,0.2), x, lognorm=T)$value # -203.1891
best(x) # Normal

Try with a lognormal distribution:

best(rlnorm(100, 2.6, 0.2)) # lognormal

Assignment will not be perfect, depending on n, mean and sd:

> table(replicate(1000, best(rnorm(500, 10, 2))))

Lognormal    Normal 
        6       994 
> table(replicate(1000, best(rlnorm(500, 2.6, 0.2))))

Lognormal    Normal 
      999         1 

Answer 2

The Bayesian approach to your problem would be to consider the posterior probability over models $M \in \{ \text{Normal}, \text{Log-normal} \}$ given a set of data points $X = \{ x_1, ..., x_N \}$,

$$P(M \mid X) \propto P(X \mid M) P(M).$$

The difficult part is getting the marginal likelihood,

$$P(X \mid M) = \int P(X \mid \theta, M) P(\theta \mid M) \, d\theta.$$

For certain choices of $p(\theta \mid M)$, the marginal likelihood of a Gaussian can be obtained in closed form. Since saying that $X$ is log-normally distributed is the same as saying that $Y = \{ \log x_1, ..., \log x_N$ } is normally distributed, you should be able to use the same marginal likelihood for the log-normal model as for the Gaussian model, by applying it to $Y$ instead of $X$. Only remember to take into account the Jacobian of the transformation,

$$P(X \mid M = \text{Log-Normal}) = P(Y \mid M=\text{Normal}) \cdot \prod_i \left| \frac{1}{x_i} \right|.$$

For this approach you need to choose a distribution over parameters $P(\theta \mid M)$ – here, presumably $P(\sigma^2, \mu \mid M=\text{Normal})$ – and the prior probabilities $P(M)$.

Example:

For $P(\mu, \sigma^2 \mid M = \text{Normal})$ I choose a normal-inverse-gamma distribution with parameters $m_0 = 0, v_0 = 20, a_0 = 1, b_0 = 100$.

According to Murphy (2007) (Equation 203), the marginal likelihood of the normal distribution is then given by

$$P(X \mid M = \text{Normal}) = \frac{|v_N|^\frac{1}{2}}{|v_0|^\frac{1}{2}} \frac{b_0^{a_0}}{b_n^{a_N}} \frac{\Gamma(a_N)}{\Gamma(a_0)} \frac{1}{\pi^{N/2}2^N}$$

where $a_N, b_N,$ and $v_N$ are the parameters of the posterior $P(\mu, \sigma^2 \mid X, M = \text{Normal})$ (Equations 196 to 200),

\begin{align} v_N &= 1 / (v_0^{-1} + N), \\ m_N &= \left( v_0^{-1}m_0 + \sum_i x_i \right) / v_N, \\ a_N &= a_0 + \frac{N}{2}, \\ b_N &= b_0 + \frac{1}{2} \left( v_0^{-1}m_0^2 - v_N^{-1}m_N^2 + \sum_i x_i^2 \right). \end{align}

I use the same hyperparameters for the log-normal distribution,

$$P(X \mid M = \text{Log-normal}) = P(\{\log x_1, ..., \log x_N \} \mid M = \text{Normal}) \cdot \prod_i \left|\frac{1}{x_i}\right|.$$

For a prior probability of the log-normal of $0.1$, $P(M = \text{Log-normal}) = 0.1$, and data drawn from the following log-normal distribution,

the posterior behaves like this:

The solid line shows the median posterior probability for different draws of $N$ data points. Note that for little to no data, the beliefs are close to the prior beliefs. For around 250 data points, the algorithm is almost always certain that the data was drawn from a log-normal distribution.

When implementing the equations, it would be a good idea to work with log-densities instead of densities. But otherwise it should be pretty straight forward. Here is the code that I used to generate the plots:

https://gist.github.com/lucastheis/6094631

Answer 3

It sounds like you are looking for something quite pragmatic to help analysts who probably aren't professional statisticians and need something to prompt them into doing what should be standard exploratory techniques such as looking at qq plots, density plots, etc.

In which case why not simply do a normality test (Shapiro-Wilk or whatever) on the original data, and one on the log transformed data, and if the second p value is higher raise a flag for the analyst to consider using a log transform? As a bonus, spit out a 2 x 2 graphic of the density line plot and qqnorm plot of the raw and the transformed data.

This won't technically answer your question about the relative likelihood but I wonder if it is all that you need.