Necesita un algoritmo para calcular la probabilidad relativa de que los datos sean muestras de distribución normal versus distribución lognormal

Supongamos que tiene un conjunto de valores y desea saber si es más probable que se muestrearon de una distribución gaussiana (normal) o de una distribución lognormal.

Por supuesto, idealmente sabría algo sobre la población o sobre las fuentes de error experimental, por lo que tendría información adicional útil para responder la pregunta. Pero aquí, supongamos que solo tenemos un conjunto de números y ninguna otra información. ¿Qué es más probable: muestreo de un gaussiano o muestreo de una distribución lognormal? ¿Cuánto más probable? Lo que espero es un algoritmo para seleccionar entre los dos modelos y, con suerte, cuantificar la probabilidad relativa de cada uno.

Puede adivinar mejor el tipo de distribución ajustando cada distribución (normal o logarítmica normal) a los datos por la máxima probabilidad, y luego comparando la probabilidad logarítmica en cada modelo: el modelo con la mayor probabilidad logarítmica es el mejor ajuste. Por ejemplo, en 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))]
}

Ahora genera números a partir de una distribución normal y ajusta una distribución normal por ML:

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

Produce:

$par
[1] 10.218083  1.787379

$value
[1] -199.9697
...

Compare la probabilidad logarítmica para el ajuste de ML de distribuciones normales y lognormales:

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

Pruebe con una distribución lognormal:

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

La asignación no será perfecta, dependiendo de n, mean y 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 

$M \in \{ \text{Normal}, \text{Log-normal} \}$$X = \{ x_1, ..., x_N \}$

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

La parte difícil es obtener la probabilidad marginal ,

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

$p(\theta \mid M)$$X$$Y = \{ \log x_1, ..., \log x_N$$Y$$X$,

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

$P(\theta \mid M)$$P(\sigma^2, \mu \mid M=\text{Normal})$$P(M)$

Ejemplo:

$P(\mu, \sigma^2 \mid M = \text{Normal})$$m_0 = 0, v_0 = 20, a_0 = 1, b_0 = 100$

Según Murphy (2007) (Ecuación 203), la probabilidad marginal de la distribución normal viene dada por

$$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}$$

$a_N, b_N,$$v_N$$P(\mu, \sigma^2 \mid X, M = \text{Normal})$

$$\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}$$

Yo uso los mismos hiperparámetros para la distribución log-normal,

$$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|.$$

$0.1$$P(M = \text{Log-normal}) = 0.1$

el posterior se comporta así:

$N$

Al implementar las ecuaciones, sería una buena idea trabajar con densidades logarítmicas en lugar de densidades. Pero de lo contrario debería ser bastante sencillo. Aquí está el código que usé para generar las tramas:

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

Parece que está buscando algo bastante pragmático para ayudar a los analistas que probablemente no sean estadísticos profesionales y que necesiten algo que los impulse a hacer lo que deberían ser técnicas exploratorias estándar, como observar gráficos qq, gráficos de densidad, etc.

En cuyo caso, ¿por qué no simplemente hacer una prueba de normalidad (Shapiro-Wilk o lo que sea) en los datos originales, y una en el registro de datos transformados, y si el segundo valor p es más alto, active un indicador para que el analista considere usar una transformación de registro ? Como beneficio adicional, escupe un gráfico de 2 x 2 de la gráfica de línea de densidad y la gráfica de qqnorm de los datos sin procesar y transformados.

Esto técnicamente no responderá a su pregunta sobre la probabilidad relativa, pero me pregunto si es todo lo que necesita.