Potrzebujesz algorytmu, aby obliczyć względne prawdopodobieństwo, że dane są próbkami z rozkładu normalnego w porównaniu z rozkładem logarytmicznym

Powiedzmy, że masz zestaw wartości i chcesz wiedzieć, czy bardziej prawdopodobne jest, że próbkowano z rozkładu Gaussa (normalnego) lub próbkowano z rozkładu logarytmicznego?

Oczywiście idealnie byłoby wiedzieć coś o populacji lub o źródłach błędów eksperymentalnych, więc mielibyśmy dodatkowe informacje przydatne w odpowiedzi na pytanie. Ale tutaj załóżmy, że mamy tylko zestaw liczb i żadnych innych informacji. Co jest bardziej prawdopodobne: pobieranie próbek z Gaussa lub pobieranie próbek z rozkładu logarytmicznego? O ile bardziej prawdopodobne? Mam nadzieję na algorytm wyboru między dwoma modelami i, mam nadzieję, ilościowe oszacowanie względnego prawdopodobieństwa każdego z nich.

Można odgadnąć typ rozkładu, dopasowując każdy rozkład (normalny lub logarytmiczny) do danych według maksymalnego prawdopodobieństwa, a następnie porównując prawdopodobieństwo dziennika dla każdego modelu - model o najwyższym prawdopodobieństwie dziennika jest najlepiej dopasowany. Na przykład w 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))]
}

Teraz generuj liczby z rozkładu normalnego i dopasuj rozkład normalny do ML:

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

Produkuje:

$par
[1] 10.218083  1.787379

$value
[1] -199.9697
...

Porównaj prawdopodobieństwo dziennika dla dopasowania ML rozkładów normalnych i logarytmicznych:

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

Spróbuj z lognormalną dystrybucją:

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

Przypisanie nie będzie idealne, w zależności od n, średniej i 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).$$

Trudną częścią jest uzyskanie marginalnego prawdopodobieństwa ,

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

Przykład:

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

Według Murphy'ego (2007) (równanie 203) krańcowe prawdopodobieństwo rozkładu normalnego podaje następnie

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

Używam tych samych hiperparametrów do rozkładu log-normalnego,

$$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$ oraz danych pochodzących z następującego rozkładu log-normalnego,

tylny zachowuje się tak:

Linia ciągła pokazuje środkowe prawdopodobieństwo z tyłu dla różnych losowań $N$

Podczas implementacji równań dobrze byłoby pracować z gęstością logarytmiczną zamiast gęstości. Ale w przeciwnym razie powinno być całkiem prosto. Oto kod, którego użyłem do wygenerowania wykresów:

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

Wygląda na to, że szukasz czegoś dość pragmatycznego, aby pomóc analitykom, którzy prawdopodobnie nie są zawodowymi statystykami i potrzebują czegoś, co zachęci ich do zrobienia czegoś, co powinno być standardowymi technikami eksploracyjnymi, takimi jak wykresy qq, wykresy gęstości itp.

W takim przypadku dlaczego nie po prostu wykonać testu normalności (Shapiro-Wilk lub cokolwiek innego) na oryginalnych danych i jednego na danych przekształconych w dzienniku, a jeśli druga wartość p jest wyższa, podnieś flagę, aby analityk rozważył użycie transformacji dziennika ? Jako bonus wypluj grafikę 2 x 2 wykresu linii gęstości i wykres qqnorm surowych i przetworzonych danych.

To technicznie nie odpowie na twoje pytanie dotyczące względnego prawdopodobieństwa, ale zastanawiam się, czy to wszystko, czego potrzebujesz.