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

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.

• It could be a fun exercise to try and characterize the distribution over distributions in nature/published literature. Then again- it will never be more than a fun exercise. For a serious treatment, you can either look for a theory justifying your choice, or given enough data- visualize and test the goodness of fit of each candidate distribution. – JohnRos Jul 21 '13 at 19:18
• If it's a matter of generalising from experience I would say that positively skewed distributions are the most common type, especially for response variables which are of central interest, and that lognormals are more common than normals. A 1962 volume The scientist speculates edited by the famous statistician I.J. Good included an anonymous piece "Bloggins's working rules", containing the assertion "The log normal distribution is more normal than the normal". (Several of the other rules are strongly statistical.) – Nick Cox Jul 21 '13 at 21:23
• I seem to interpret your question differently from JohnRos and anxoestevez. To me, your question sounds like one about plain model selection, that is, a matter of computing $P(M \mid D)$, where $M$ is either the normal or log-normal distribution and $D$ is your data. If model selection is not what you're after, can you clarify? – Lucas Jul 21 '13 at 21:28
• @lucas I think your interpretation is not so much different from mine. In either case you need to do apriori assumptions. – anxoestevez Jul 21 '13 at 21:33
• Why not just calculate the generalized likelihood ratio & alert the user when it favours the log-normal? – Scortchi - Reinstate Monica Jul 24 '13 at 1:46

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

• You don't need to find the maximum-likelihood parameter estimates numerically for either the normal or log-normal (though it shows how you'd generalize the idea to comparison of other distributions). Apart from that, very sensible approach. – Scortchi - Reinstate Monica Jul 25 '13 at 12:46
• I've barely used R or the concept of maximum likelihood, so here is a basic question. I know we can't compare the AIC (or BIC) from fitting a normal distribution to the data vs. to the logs of the data, because the AIC or BIC would not be comparable. One needs to fit two models to one set of data (with no transforms; no outlier exclusions etc.), and transforming the data will change AIC or BIC regardless making the comparison bogus. What about ML? Is this comparison legit? – Harvey Motulsky Jul 26 '13 at 18:44
• We find the best fitting normal and lognormal distributions to the data, then calculate the probability of observing the data assuming they were from those distributions (the likelihood or p(X|\theta)). We're not transforming the data. We print out the distribution for which the probability of observing the data is highest. This approach is legit but has the disadvantage that we don't infer the probability of the model given the data p(M|X), i.e. the probability that the data is from a normal vs lognormal distribution (e.g. p(normal) = 0.1, p(lognormal) = 0.9) unlike the Bayesian approach. – waferthin Jul 27 '13 at 8:48
• @Harvey True enough, but irrelevant - you asked about fitting normal vs log-normal distributions to the same data, & this is what whannymahoots is answering. Because the number of free parameters is the same for both models, comparing AICs or BICs reduces to comparing log-likelihoods. – Scortchi - Reinstate Monica Jul 27 '13 at 19:02
• @wannymahoots Any reasonable prior for a Bayesian approach in this context - relying on estimating the relative probabilities that a software user is trying to fit normal or log-normal data - is going to be so uninformative that it'll give similar results to an approach based just on likelihood. – Scortchi - Reinstate Monica Jul 27 '13 at 19:41

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

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.

• Clever. Maybe this is enough, and avoids the need to explain likelihood calculations.... Thanks. – Harvey Motulsky Jul 26 '13 at 18:44