Necesita el algoritmo para calcular la probabilidad relativa de que los datos son una muestra de lo normal vs distribución lognormal

Question

Digamos que usted tiene un conjunto de valores, y desea saber si es más probable que ellos se tomaron muestras de una distribución Gaussiana (normal) o de muestreo de una distribución logarítmico-normal?

Por supuesto, lo ideal sería que usted sabría algo acerca de la población o sobre las fuentes de error experimental, por lo que habría información adicional útil para responder a la pregunta. Pero aquí, supongamos que sólo tenemos un conjunto de números y no de otra información. Lo que es más probable: el muestreo de una Gaussiana o el muestreo de una distribución logarítmico-normal? Cuánto más probable? Lo que estoy esperando es un algoritmo para seleccionar entre los dos modelos, y esperemos que cuantificar la probabilidad relativa de cada uno.

Answer 1

El enfoque Bayesiano para su problema sería considerar la probabilidad posterior sobre los modelos de $M \in \{ \text{Normal}, \text{Log-normal} \}$ dado un conjunto de puntos de datos $X = \{ x_1, ..., x_N \}$,

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

La parte difícil es conseguir que los marginales de la probabilidad,

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

Para ciertas opciones de $p(\theta \mid M)$, la probabilidad marginal de una Gaussiana puede ser obtenida en forma cerrada. Ya diciendo que $X$ es de registro-normalmente distribuida es lo mismo que decir que $Y = \{ \log x_1, ..., \log x_N$ } se distribuye normalmente, usted debería ser capaz de utilizar el mismo marginales de la probabilidad de la log-normal como modelo para el modelo de Gauss, mediante la aplicación a $Y$ en lugar de $X$. Sólo recuerde tomar en cuenta el Jacobiano de la transformación,

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

Para este método se necesita elegir una distribución a través de los parámetros de $P(\theta \mid M)$ – aquí, presumiblemente $P(\sigma^2, \mu \mid M=\text{Normal})$ – y las probabilidades previas $P(M)$.

Ejemplo:

Para $P(\mu, \sigma^2 \mid M = \text{Normal})$ I elegir una normal inverso-distribución gamma con parámetros a $m_0 = 0, v_0 = 20, a_0 = 1, b_0 = 100$.

Según Murphy (2007) (Ecuación 203), los marginales de la probabilidad de la distribución normal está dado 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}$$

donde $a_N, b_N,$ $v_N$ son los parámetros de la parte posterior de la $P(\mu, \sigma^2 \mid X, M = \text{Normal})$ (Ecuaciones 196 a 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}

Yo uso el mismo hyperparameters para la log-normal de distribución,

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

Para una probabilidad anterior de la log-normal de $0.1$, $P(M = \text{Log-normal}) = 0.1$, y los datos extraídos de la siguiente log-normal de distribución,

la parte posterior se comporta como esta:

La línea sólida muestra la mediana de probabilidad posterior para los diferentes sorteos de $N$ puntos de datos. Tenga en cuenta que por muy poco no los datos, las creencias son cerca de las creencias anteriores. De alrededor de 250 puntos de datos, el algoritmo es casi siempre la certeza de que los datos fueron extraídos de una log-normal de distribución.

Cuando la aplicación de las ecuaciones, sería una buena idea trabajar con el registro de densidades en lugar de densidades. Pero de lo contrario, debería ser bastante sencillo. Aquí está el código que he utilizado para generar las parcelas:

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

Answer 2

Usted podría tomar una mejor respuesta en el tipo de distribución por ajuste de cada distribución (normal o lognormal) a los datos de máxima verosimilitud, a continuación, comparar el logaritmo de la probabilidad de cada modelo - el modelo con el mayor registro de la probabilidad de ser 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 generar números a partir de una distribución normal, y el ajuste a 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
...

Comparar la log-verosimilitud para ML ajuste de normal y lognormal distribuciones:

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 logarítmico-normal:

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

Asignación no será perfecto, dependiendo de n, media 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 

Answer 3

Suena como que usted está buscando algo muy pragmático para ayudar a los analistas que probablemente no son los estadísticos profesionales y necesitan algo para obligarlos a hacer lo que deberían ser la norma de técnicas de exploración como mirar gráficos qq, la densidad de terrenos, etc.

En cuyo caso, ¿por qué no simplemente hacer un test de normalidad (Shapiro-Wilk o lo que sea) en los datos originales, y uno en el registro de datos transformados, y si el segundo valor de p es mayor levantar una bandera para el analista considere el uso de un registro de transformar? Como un bono, escupir un 2 x 2 gráfica de la densidad de línea de parcela y qqnorm parcela de la materia prima y los datos transformados.

Esto no técnicamente responder a su pregunta acerca de la probabilidad relativa pero me pregunto si eso es todo lo que usted necesita.