Perlu algoritma untuk menghitung kemungkinan relatif bahwa data adalah sampel dari distribusi normal vs lognormal

Katakanlah Anda memiliki seperangkat nilai, dan Anda ingin tahu apakah lebih besar kemungkinannya diambil sampelnya dari distribusi Gaussian (normal) atau sampel dari distribusi lognormal?

Tentu saja, idealnya Anda akan tahu sesuatu tentang populasi atau tentang sumber kesalahan eksperimental, sehingga akan memiliki informasi tambahan yang berguna untuk menjawab pertanyaan. Tapi di sini, anggaplah kita hanya memiliki satu set angka dan tidak ada informasi lain. Mana yang lebih mungkin: pengambilan sampel dari Gaussian atau pengambilan sampel dari distribusi lognormal? Berapa besar kemungkinannya? Apa yang saya harapkan adalah algoritma untuk memilih antara dua model, dan mudah-mudahan menghitung kemungkinan masing-masing model.

Anda dapat mengambil tebakan terbaik pada tipe distribusi dengan memasang setiap distribusi (normal atau lognormal) ke data dengan kemungkinan maksimum, kemudian membandingkan kemungkinan log di setiap model - model dengan kemungkinan log tertinggi yang paling cocok. Misalnya, dalam 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))]
}

Sekarang hasilkan angka dari distribusi normal dan paskan distribusi normal dengan ML:

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

Menghasilkan:

$par
[1] 10.218083  1.787379

$value
[1] -199.9697
...

Bandingkan kemungkinan log untuk ML fit dari distribusi normal dan lognormal:

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

Coba dengan distribusi lognormal:

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

Tugas tidak akan sempurna, tergantung pada n, mean dan 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).$$

Bagian yang sulit adalah mendapatkan kemungkinan marjinal ,

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

Contoh:

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

Menurut Murphy (2007) (Persamaan 203), kemungkinan marginal dari distribusi normal kemudian diberikan oleh

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

Saya menggunakan hyperparameters yang sama untuk distribusi 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$

posterior berperilaku seperti ini:

$N$

Ketika menerapkan persamaan, itu akan menjadi ide yang baik untuk bekerja dengan kepadatan log, bukan kepadatan. Tetapi sebaliknya itu harus lurus ke depan. Berikut adalah kode yang saya gunakan untuk membuat plot:

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

Sepertinya Anda mencari sesuatu yang cukup pragmatis untuk membantu analis yang mungkin bukan ahli statistik profesional dan membutuhkan sesuatu untuk mendorong mereka melakukan apa yang seharusnya menjadi teknik eksplorasi standar seperti melihat plot qq, plot kepadatan, dll.

Dalam hal ini mengapa tidak hanya melakukan tes normalitas (Shapiro-Wilk atau apa pun) pada data asli, dan satu pada log mengubah data, dan jika nilai p kedua lebih tinggi menaikkan bendera untuk analis untuk mempertimbangkan menggunakan log transformasi ? Sebagai bonus, keluarkan grafik 2 x 2 dari plot garis kerapatan dan plot qqnorm dari data mentah dan yang diubah.

Ini tidak akan secara teknis menjawab pertanyaan Anda tentang kemungkinan relatif tetapi saya ingin tahu apakah itu yang Anda butuhkan.