Использование pymc3 для соответствия распределению Стьюдента

Не уверен, что я делаю что-то глупое или у pymc3 есть ошибка, но пытаясь подогнать T-распределение к нормальному, я получаю количество степеней свободы (от 0,18 до 0,25, я ожидаю чего-то высокого, 4-5 по крайней мере). Конечно, я получаю ту же ошибку, если пробую T-распределение с разумным количеством степеней свободы, например, с 3 или 5.

import pymc3 as pm
Nsample = 200000
tst = np.random.normal(loc = 1e4, scale = 5e4, size = 250)
with pm.Model() as m:
    mean = pm.Normal('mean',mu=0,sd = 1e5)
    sigma = pm.Flat('sigma') # I tried uniform, gamma, exponential
    df = pm.Flat("df") # the same
    v = pm.T("pl",nu=df,mu = mean, lam = 1.0/sigma, observed = tst)

    start = {'df':5,'mean': 1e4, 'sigma':5e4} #start = pm.find_MAP()

    step = pm.Metropolis()
    trace = pm.sample(Nsample, step,start=start, progressbar=True)

pm.traceplot(trace[100000:],vars = ['df', 'sigma', 'mean']);

Не могли бы вы предложить какое-то исправление (изменение априоров, метод выборки)?


bayesian mcmc pymc3 pymc
person Alexey Goldin    schedule 21.10.2014    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Почему вы ожидаете увидеть df около 4-5? T-распределение с df->inf равно нормальному распределению. Когда я запускаю вашу модель и делаю: print trace['df'][10000:].mean(), я получаю 1.19876070951e+13, то есть что-то очень большое.

Одна из причин, по которой вы можете увидеть что-то другое, заключается в том, что сэмплер Metropolis, скорее всего, потерпит неудачу, если вы попытаетесь выполнить семплирование в совместном пространстве (которое раньше использовалось по умолчанию в pymc3). Если вы недавно не обновляли pymc3 из мастера, попробуйте обновить и снова запустить модель, так как Metropolis теперь по умолчанию неблокирует и производит выборку каждой переменной отдельно.

person twiecki    schedule 23.10.2014
comment
Теперь он отлично работает, даже если я заменю tst на tst = np.random.standard_t(3.0,size=250)*20000+1e4, как я изначально планировал. Спасибо! Пока не понимаю, почему предыдущая версия не работала. Вернуться к учебе. - person Alexey Goldin; 24.10.2014
comment
NUTS по-прежнему дает нелепые ответы (df ~ 0,22) - person Alexey Goldin; 24.10.2014
comment
Вы используете start = find_MAP(); sampler = NUTS(scaling=start)? Они почти всегда требуются. - person twiecki; 25.10.2014
comment
Прошу прощения, NUTS работает. Моя проблема заключалась в том, что мой априор был слишком ограничен для параметра масштаба в распределении Т Стьюдента, поскольку я не понял, что параметр масштаба составляет ~ 1/квадрат ширины квадрата. Слишком ограниченный априорный масштаб привел к очень низким значениям степеней свободы, поскольку это был единственный способ получить такие значения в хвостах. Теперь работает NUTS, мне пришлось заменить v = pm.T(pl,nu=df,mu = среднее, lam = 1,0/сигма, наблюдаемое = tst) на v = pm.T(pl,nu=df,mu = среднее , lam = 1,0/(сигма*сигма), наблюдаемое = tst) и pm.Flat с pm.Uniform или pm.Exponential Спасибо! - person Alexey Goldin; 26.10.2014