Нужен алгоритм для вычисления относительной вероятности того, что данные являются выборкой из нормального и логнормального распределения

Допустим, у вас есть набор значений, и вы хотите знать, более ли вероятно, что они были выбраны из гауссова (нормального) распределения или из логнормального распределения?

Конечно, в идеале вы должны были бы что-то знать о населении или об источниках экспериментальной ошибки, поэтому имели бы дополнительную информацию, полезную для ответа на вопрос. Но здесь, предположим, у нас есть только набор чисел и никакой другой информации. Что является более вероятным: выборка по Гауссу или выборка по логнормальному распределению? Насколько более вероятно? Я надеюсь на алгоритм выбора между двумя моделями и, надеюсь, количественную оценку относительной вероятности каждой из них.

Вы можете сделать наиболее точную оценку типа распределения, подгоняя каждое распределение (нормальное или логнормальное) к данным с максимальной вероятностью, а затем сравнивая логарифмическое правдоподобие для каждой модели - модель с наибольшим логарифмическим правдоподобием, которая наилучшим образом подходит. Например, в 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))]
}

Теперь сгенерируйте числа из нормального распределения и подгоните нормальное распределение по ML:

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

Производит:

$par
[1] 10.218083  1.787379

$value
[1] -199.9697
...

Сравните правдоподобие для соответствия ML нормального и логнормального распределений:

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

Попробуйте с логнормальным дистрибутивом:

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

Назначение не будет идеальным, в зависимости от n, среднего и сд:

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

Трудная часть заключается в получении предельной вероятности ,

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

Пример:

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

Согласно Мерфи (2007) (уравнение 203), предельная вероятность нормального распределения определяется как

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

Для предварительной вероятности нормального логарифма $0.1$, $P(M = \text{Log-normal}) = 0.1$и данные, взятые из следующего лог-нормального распределения,

задний ведет себя так:

Сплошная линия показывает среднюю апостериорную вероятность для разных розыгрышей $N$Точки данных. Обратите внимание на то, что по небольшим или никаким данным убеждения близки к предыдущим убеждениям. Приблизительно для 250 точек данных алгоритм почти всегда уверен, что данные были получены из лог-нормального распределения.

При реализации уравнений было бы неплохо работать с логарифмическими плотностями вместо плотностей. Но в остальном все должно быть довольно просто. Вот код, который я использовал для создания графиков:

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

Похоже, вы ищете что-то весьма прагматичное, чтобы помочь аналитикам, которые, вероятно, не являются профессиональными статистиками и нуждаются в чем-то, что побуждает их делать то, что должно быть стандартными исследовательскими методами, такими как просмотр графиков qq, графиков плотности и т. Д.

В таком случае, почему бы просто не выполнить тест нормальности (Shapiro-Wilk или любой другой) для исходных данных и один для данных, преобразованных в журнал, и, если второе значение p выше, поднимите флаг для аналитика, чтобы рассмотреть возможность использования преобразования журнала ? В качестве бонуса выложите график 2 x 2 графика линейной плотности и график qqnorm необработанных и преобразованных данных.

Технически это не ответит на ваш вопрос об относительной вероятности, но мне интересно, если это все, что вам нужно.