Создание MLE для пары дистрибутивов на Python

Итак, в моем текущем коде подбора кривой есть шаг, который использует scipy.stats для определения правильного распределения на основе данных,

distributions = [st.laplace, st.norm, st.expon, st.dweibull, st.invweibull, st.lognorm, st.uniform]
mles = []

for distribution in distributions:
    pars = distribution.fit(data)
    mle = distribution.nnlf(pars, data)
    mles.append(mle)

results = [(distribution.name, mle) for distribution, mle in zip(distributions, mles)]

for dist in sorted(zip(distributions, mles), key=lambda d: d[1]):
    print dist
best_fit = sorted(zip(distributions, mles), key=lambda d: d[1])[0]
print 'Best fit reached using {}, MLE value: {}'.format(best_fit[0].name, best_fit[1])          


print [mod[0].name for mod in sorted(zip(distributions, mles), key=lambda d: d[1])]

Где данные - это список числовых значений. Пока что это отлично работает для подгонки одномодальных распределений, подтвержденных в скрипте, который случайным образом генерирует значения из случайных распределений и использует curve_fit для переопределения параметров.

«Соответствующее

Теперь я хотел бы сделать код способным обрабатывать бимодальные распределения, как в примере ниже:

«Комбинация

Можно ли получить MLE для пары моделей из scipy.stats, чтобы определить, подходит ли конкретная пара распределений для данных?, Что-то вроде

distributions = [st.laplace, st.norm, st.expon, st.dweibull, st.invweibull, st.lognorm, st.uniform]
distributionPairs = [[modelA.name, modelB.name] for modelA in distributions for modelB in distributions]

и использовать эти пары, чтобы получить значение MLE этой пары распределений, соответствующих данным?


python scipy statistics model-fitting distribution
person BruceJohnJennerLawso    schedule 03.02.2017    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Это не полный ответ, но он может помочь вам решить вашу проблему. Допустим, вы знаете, что ваша проблема вызвана двумя плотностями. Решением может быть использование k-среднего или алгоритма EM.

Инициализация. Вы инициализируете свой алгоритм, изменяя каждое наблюдение до той или иной плотности. И вы инициализируете две плотности (вы инициализируете параметры плотности, и один из параметров в вашем случае - «гауссов», «лаплас» и так далее ... Итерация. Затем, итеративно, вы выполняете два следующих шага :

Шаг 1. Оптимизируйте параметры, предполагая, что влияние каждой точки правильное. Теперь вы можете использовать любой решатель оптимизации. Этот шаг дает вам оценку двух лучших плотностей (с заданным параметром), которые соответствуют вашим данным.

Шаг 2. Вы относите каждое наблюдение к той или иной плотности в соответствии с наибольшей вероятностью.

Повторяете до схождения.

Это очень хорошо объясняется на этой веб-странице https://people.duke.edu/~ccc14/sta-663/EMAlgorithm.html.

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

Вот пример кодирования в простом случае: вы знаете, что ваши данные поступают из двух разных гауссиан (вы не знаете, сколько переменных генерируется из каждой плотности). В вашем случае вы можете настроить этот код для цикла на каждой возможной паре плотности (в вычислительном отношении дольше, но я полагаю, что он будет работать эмпирически)

import scipy.stats as st
import numpy as np

#hard coded data generation
data = np.random.normal(-3, 1, size = 1000)
data[600:] = np.random.normal(loc = 3, scale = 2, size=400)

#initialization

mu1 = -1
sigma1 = 1

mu2 = 1
sigma2 = 1

#criterion to stop iteration
epsilon = 0.1
stop = False

while  not stop :  
    #step1
    classification = np.zeros(len(data))
    classification[st.norm.pdf(data, mu1, sigma1) > st.norm.pdf(data, mu2, sigma2)] = 1

    mu1_old, mu2_old, sigma1_old, sigma2_old = mu1, mu2, sigma1, sigma2

    #step2
    pars1 = st.norm.fit(data[classification == 1])
    mu1, sigma1 = pars1
    pars2 = st.norm.fit(data[classification == 0])
    mu2, sigma2 = pars2

    #stopping criterion
    stop = ((mu1_old - mu1)**2 + (mu2_old - mu2)**2 +(sigma1_old - sigma1)**2 +(sigma2_old - sigma2)**2) < epsilon

#result    
print("The first density is gaussian :", mu1, sigma1)
print("The first density is gaussian :", mu2, sigma2)
print("A rate of ", np.mean(classification), "is classified in the first density")

Надеюсь, поможет.

person probaPerception    schedule 09.02.2017
comment
Большое спасибо, похоже, это хорошо работает. Я не уверен, что понимаю, как работает код? Похоже, что он итеративно подгоняет две разные нормальные кривые, сортируя набор данных на два отдельных списка (или, скорее, используя классификацию в качестве массива индикаторов, в какую категорию попадает каждая точка данных? Это потрясающе, я понятия не имел, что вы можете сделать это с массивы numpy). Для случаев, когда дистрибутивы хорошо разделены, кажется, что это хорошо работает: i.imgur.com/8Hrhd0F.png - person BruceJohnJennerLawso; 14.02.2017
comment
Для дистрибутивов, которые не так хорошо разделены, я замечаю, что цикл имеет тенденцию пытаться заставить решение, которое распространяется, например здесь и особенно здесь. Я предполагаю, что это связано с начальными условиями, начинающимися с идентичных сигм и средств распространения, возможно, имеет смысл выполнить несколько прогонов при подборе пары распределений с разными начальными значениями для mu1 / 2 / sigma1 / 2 и сравнить окончательные p ценности. - person BruceJohnJennerLawso; 14.02.2017
comment
Последнее, что я пытаюсь понять, это как совместить мультимодальные перевозки помимо двухрежимных. Я думал о том, чтобы сделать что-то вроде рекурсивной вещи, когда для трех нормальных кривых цикл соответствует одному из распределений, соответствует нормали по оставшимся двум, затем оставшиеся два идентифицируются как действительно плохо подходящие, и цикл выполняется как обычно на них. Но похоже, что не так уж и хороша, даже когда дистрибутивы хорошо разделены. - person BruceJohnJennerLawso; 14.02.2017
comment
Что касается вашего первого комментария, если две плотности не сильно различаются, тогда очень трудно получить хорошие результаты. Но это имеет смысл, поскольку трудно разделить данные от одной плотности к другой. Что касается ваших 3 плотностей, я думаю, что лучший способ решить эту проблему - запустить тот же алгоритм, но с 3 потенциальными плотностями вместо 2. Если вам нужен общий подход с неизвестным числом плотностей, вам нужно оптимизировать более сложную задачу, где вы устанавливаете компромисс между качеством соответствия и количеством плотностей, используемых для представления ваших данных. - person probaPerception; 15.02.2017
comment
То есть это что-то вроде числового применения бритвы Оккама, которое штрафует значение тестовой статистики совпадений с более высоким общим числом распределений? Можете ли вы описать, что именно происходит в задаче классификации штрафов? - person BruceJohnJennerLawso; 18.02.2017
comment
По сути, ваш алгоритм будет пытаться свести к минимуму степень согласия с минимально возможными плотностями. Это немного длинновато для объяснения в комментарии, но если вам интересно, вы можете прочитать этот PDF-файл: statweb.stanford.edu/~jtaylo/courses/stats203/notes/. Речь идет не о классификации, а скорее о регрессии. Однако он очень хорошо объясняет концепцию штрафов. Адаптация к классификации очень похожа (хотя теоретических результатов по ней в литературе меньше [с моей точки зрения]). - person probaPerception; 18.02.2017