проблемы с началом работы с простым примером pymc3

Question

Я новичок в использовании пакета PyMC3 и просто пытаюсь реализовать пример из курса по неопределенности измерений, который я беру.(Обратите внимание, что это необязательный курс обучения сотрудников через работу, а не класс, где я не должен находить ответы в Интернете).В курсе используется R, но я считаю, что Python предпочтительнее.

(простая) задача ставится следующим образом:

Скажем, у вас есть конечный измеритель фактической (неизвестной) длины в комнате.температура length и измеренная длина m.Соотношение между ними:

length = m / (1 + alpha*dT)

, где alpha - коэффициент расширения, а dT - отклонение от комнатной температуры, а m - измеренная величина.Цель состоит в том, чтобы найти апостериорное распределение на length, чтобы определить его ожидаемое значение и стандартное отклонение (т. Е. Неопределенность измерения)

Задача задает предыдущие распределения по альфа и dT (гауссианы с небольшим стандартным отклонением)и свободный предшествующий length (гауссовский с большим стандартным отклонением).Проблема указывает, что m был измерен 25 раз со средним значением 50.000215 и стандартным отклонением 5.8e-6.Мы предполагаем, что измерения m обычно распределяются со средним значением истинного значения m.

Одна из проблем, с которыми я столкнулся, заключается в том, что вероятность, по-видимому, не может быть определена только на основеэти статистические данные в PyMC3, поэтому я сгенерировал некоторые фиктивные данные измерений (в итоге я сделал 1000 измерений вместо 25).Опять же, вопрос заключается в том, чтобы получить апостериорное распределение на length (и в процессе, хотя и менее интересном, обновленные постеры на alpha и dT).

Вот мой код, который неработа и проблемы сходимости:

from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as stats
import pymc3 as pm
import theano.tensor as tt

basic_model = pm.Model()

xdata = np.random.normal(50.000215,5.8e-6*np.sqrt(1000),1000)

with basic_model:
    #prior distributions
    theta = pm.Normal('theta',mu=-.1,sd=.04)
    alpha = pm.Normal('alpha',mu=.0000115,sd=.0000012)
    length = pm.Normal('length',mu=50,sd=1)

mumeas = length*(1+alpha*theta)


with basic_model:
    obs = pm.Normal('obs',mu=mumeas,sd=5.8e-6,observed=xdata)
    #yobs = Normal('yobs',)
    start = pm.find_MAP()
    #trace = pm.sample(2000, step=pm.Metropolis, start=start)
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=200000,step=step,start=start,njobs=1)


length_samples = trace['length']

fig,ax=plt.subplots()
plt.hist(length_samples, histtype='stepfilled', bins=30, alpha=0.85,
     label="posterior of $\lambda_1$", color="#A60628", normed=True)

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

library(rjags)

#Data
jags_data <- list(xbar=50.000215) 

jags_code <- jags.model(file = "calibration.txt",
                    data = jags_data, 
                    n.chains = 1, 
                    n.adapt = 30000)


post_samples <- coda.samples(model = jags_code, 
                         variable.names = 
c("l","mu","alpha","theta"),#,"ypred"),
                         n.iter = 30000)


summary(post_samples)

mean(post_samples[[1]][,"l"])
sd(post_samples[[1]][,"l"])

plot(post_samples)

и модель calib.txt:

model{
  l~dnorm(50,1.0)
  alpha~dnorm(0.0000115,694444444444)
  theta~dnorm(-0.1,625)
  mu<-l*(1+alpha*theta)
  xbar~dnorm(mu,29726516052)
}

(заметьте, я думаю, что распределение dnorm занимает 1 / сигма ^ 2,отсюда и странные отклонения)

Любая помощь или понимание того, почему выборка PyMC3 не сходится и что я должен делать по-другому, была бы чрезвычайно признательна.Спасибо!

Answer 1

У меня также были проблемы с получением чего-либо полезного из сгенерированных данных и модели в коде. Мне кажется, что уровень шума в фальшивых данных в равной степени можно объяснить различными источниками дисперсии в модели. Это может привести к ситуации с сильно коррелированными задними параметрами. Добавьте к этому крайние масштабы дисбаланса, тогда имеет смысл, что это будет иметь проблемы с выборкой.

Однако, глядя на модель JAGS, кажется, что они действительно используют только это одно входное наблюдение . Я никогда не видел этот метод (?) Раньше, то есть ввод сводной статистики данных вместо самих необработанных данных. Я полагаю, что это сработало для них в JAGS, поэтому я решил попробовать запустить точно такой же MCMC, в том числе с использованием точной (tau) параметризации гауссиана.

Оригинальная модель с метрополисом

with pm.Model() as m0:
    # tau === precision parameterization
    dT = pm.Normal('dT', mu=-0.1, tau=625)
    alpha = pm.Normal('alpha', mu=0.0000115, tau=694444444444)
    length = pm.Normal('length', mu=50.0, tau=1.0)

    mu = pm.Deterministic('mu', length*(1+alpha*dT))

    # only one input observation; tau indicates the 5.8 nm sd
    obs = pm.Normal('obs', mu=mu, tau=29726516052, observed=[50.000215])

    trace = pm.sample(30000, tune=30000, chains=4, cores=4, step=pm.Metropolis())

Несмотря на то, что выборка length и dT все еще не так хороша, она, по крайней мере, в целом выглядит сходящейся:

Я думаю, что здесь следует отметить, что, несмотря на относительно слабый априор length (sd=1), сильные априоры по всем другим параметрам, по-видимому, распространяют жесткую неопределенность, ограничивающую апостериор length. В конечном счете, это является последним интересом, так что, похоже, это согласуется с целью упражнения. Также обратите внимание на то, что mu выходит сзади, точно так же, как описано распределение, а именно N(50.000215, 5.8e-6).

График трасс

Лесной участок

Парный участок

Здесь, однако, вы можете видеть, что основная проблема все еще существует. Между стандартными ошибками существует сильная корреляция между length и dT, а также разница в 4 или 5 порядков шкалы. Я определенно проделал бы долгий путь, прежде чем я действительно доверился результату.

Альтернативная модель с гайками

Чтобы запустить это с NUTS, вам нужно решить проблему с масштабированием. То есть, каким-то образом нам нужно провести повторную параметризацию, чтобы приблизить все значения tau к 1. Затем вы запустили бы сэмплер и преобразовали его обратно в интересующие вас единицы. К сожалению, у меня нет времени играть с этим прямо сейчас (я бы тоже должен был это выяснить), но, возможно, это то, что вы можете начать изучать самостоятельно.