Блог

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

Мои данные структурированы как серия записей. Каждая запись содержит пару наблюдений, которые относятся к фиксированному периоду в один час. Одно наблюдение — это общее количество событий, которые происходят в течение данного часа. Другое наблюдение — это количество успехов за этот период времени. Так, например, в точке данных может быть указано, что за данный 1 часовой период всего произошло 1000 событий, и из этих 1000 100 были успешными. В другой период времени может быть всего 1000000 событий, из которых 120000 являются успешными. Дисперсия наблюдений не является постоянной и зависит от общего количества событий, и отчасти этот эффект я хотел бы контролировать и моделировать.

Мой первый шаг в этом, чтобы оценить базовый показатель успеха. Я подготовил приведенный ниже код, который призван имитировать ситуацию, предоставляя два набора «наблюдаемых» данных, генерируя их с помощью scipy. Однако он не работает должным образом.
Я ожидаю, что он найдет:

• loss_lambda_factor равен примерно 0,1
• total_lambda (и total_lambda_mu) составляет примерно 120.

Вместо этого модель сходится очень быстро, но к неожиданному ответу.

• total_lambda и total_lambda_mu соответственно являются резкими пиками около 5e5.
• loss_lambda_factor равен примерно 0.

График трассировки (который я не могу опубликовать из-за репутации ниже 10) довольно неинтересен — быстрая конвергенция и резкие пики при числах, которые не соответствуют входным данным. Мне любопытно, есть ли что-то принципиально неправильное в подходе, который я использую. Как следует изменить следующий код, чтобы получить правильный / ожидаемый результат?

 from pymc import Model, Uniform, Normal, Poisson, Metropolis, traceplot 
from pymc import sample 
import scipy.stats

totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000)
totalCounts = scipy.stats.poisson.rvs(mu=totalRates) 
successRate = 0.1*totalRates 
successCounts = scipy.stats.poisson.rvs(mu=successRate) 

with Model() as success_model: 
    total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000)
    total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000)
    total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
    total = Poisson('total', mu=total_lambda, observed=totalCounts) 

    loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
    success_rate = Poisson('success_rate', mu=total_lambda*loss_lambda_factor, observed=successCounts) 

with success_model: 
    step =  Metropolis() 
    success_samples = sample(20000, step) #, start)


plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples)
  

В вашем подходе нет ничего принципиально неправильного, за исключением подводных камней любого байесовского анализа MCMC: (1) несходимость, (2) приоры, (3) модель.

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

Это не очень хорошо, и чтобы более четко понять, почему, я бы изменил код схемы трассировки, чтобы показывать только вторую половину трассировки, traceplot(success_samples[10000:]) :

Предыдущее: одной из основных проблем для конвергенции является ваше предыдущее включение total_lambda_tau , которое является типичной ошибкой в байесовском моделировании. Хотя использование prior может показаться довольно неинформативным Uniform('total_lambda_tau', lower=0, upper=100000) , результатом этого является утверждение, что вы совершенно уверены, что total_lambda_tau оно велико. Например, вероятность того, что оно меньше 10, равна .0001. Изменение предшествующего

 total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100)
total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000)
  

в результате получается схема трассировки, которая является более многообещающей:

Однако это все еще не то, что я ищу на графике трассировки, и чтобы получить что-то более удовлетворительное, я предлагаю использовать шаг «последовательное сканирование мегаполиса» (который по умолчанию используется в PyMC2 для аналогичной модели). Вы можете указать это следующим образом:

 step =  pm.CompoundStep([pm.Metropolis([total_lambda_mu]),
                         pm.Metropolis([total_lambda_tau]),
                         pm.Metropolis([total_lambda]),
                         pm.Metropolis([loss_lambda_factor]),
                         ]) 
  

Это создает график трассировки, который кажется приемлемым:

Модель: Как ответил @KaiLondenberg, подход, который вы использовали с приоритетами total_lambda_tau и total_lambda_mu не является стандартным подходом. Вы описываете сильно различающиеся итоговые значения событий (1000 в один час и 1 000 000 в следующий), но ваша модель утверждает, что она распределена нормально. В пространственной эпидемиологии подход, который я видел для аналогичных данных, представляет собой модель, более похожую на эту:

 import pymc as pm, theano.tensor as T
with Model() as success_model: 
    loss_lambda_rate = pm.Flat('loss_lambda_rate')
    error = Poisson('error', mu=totalCounts*T.exp(loss_lambda_rate), 
            observed=successCounts)
  

Я уверен, что есть и другие способы, которые покажутся более знакомыми и в других исследовательских сообществах.

Вот записная книжка, в которой собраны эти комментарии.

Я вижу несколько потенциальных проблем с моделью.

1.) Я бы подумал, что подсчет успеха (называемый ошибкой?) Должен следовать биномиальному (n = total, p = loss_lambda_factor) распределению, а не Пуассону.

2.) С чего начинается цепочка? Начинать с конфигурации MAP или MLE имеет смысл, если вы не используете чистую выборку Гиббса. В противном случае для записи цепочки может потребоваться много времени, что может быть тем, что здесь происходит.

3.) Ваш выбор иерархического предшествующего значения для total_lambda (т.Е. Нормального с двумя одинаковыми значениями для этих параметров) гарантирует, что цепочке потребуется много времени для сходимости, если вы не выберете разумное начало (как в пункте 2.).). По сути, вы вводите множество ненужных степеней свободы, в которых цепочка MCMC может потеряться. Учитывая, что total_lambda должен быть неотрицательным, я бы выбрал единообразное значение для total_lambda в подходящем диапазоне (например, от 0 до наблюдаемого максимума).

4.) Вы используете пробоотборник Metropolis. для этого может быть недостаточно 20000 выборок. Попробуйте 60000 и отбросьте первые 20000 как записанные. Пробоотборнику Metropolis может потребоваться некоторое время для настройки размера шага, поэтому вполне может быть, что первые 20000 выборок были потрачены в основном на отклонение предложений и настройку. Попробуйте другие сэмплеры, такие как NUTS.