Как построить байесовскую имитационную модель для трехкратного подбрасывания монеты - PullRequest
Новая
       122

Как построить байесовскую имитационную модель для трехкратного подбрасывания монеты

0 голосов
/ 13 февраля 2019

Представьте, что мы подбросили смещенную монету 8 раз (мы не знаем, насколько она смещена), и мы до сих пор записали от 5 голов (H) до 3 хвостов (T).Какова вероятность того, что следующие 3 броска будут хвостами?Другими словами, нас интересует ожидаемая вероятность наличия 5Hs и 6Ts после 11-го броска.

Я хочу построить имитационную модель MCMC, используя pyMC3, чтобы найти байесовское решение.Существует также аналитическое решение в рамках байесовского подхода к этой проблеме.Таким образом, я смогу сравнить результаты, полученные в результате симуляции, аналитического метода, а также классического метода измерения частоты.Позвольте мне кратко объяснить, что я мог сделать до сих пор:

  1. Частотное решение:

Если мы рассмотрим вероятность одного броска: E (T) = p = (3/8) = 0,375 Тогда окончательный ответ E ({T, T, T}) = p ^ 3 = (3/8) ^ 3 = 0,052.

2.1.Байесовское решение с аналитическим способом:

Пожалуйста, предположите, что неизвестный параметр «p» представляет для смещения монеты.Если мы рассмотрим вероятность одного броска: E (T) = Integral0-1 [p * P (p | H = 5, T = 3) dp] = 0,400 (я вычислил результат после некоторой алгебраической манипуляции) Точно так жеокончательный ответ: E ({T, T, T}) = Integral0-1 [p ^ 3 * P (p | H = 5, T = 3) dp] = 10/11 = 0,909.

2.2.Байесовское решение с симуляцией MCMC: Когда мы рассматриваем вероятность для одного броска, я построил модель в pyMC3, как показано ниже: 

Head: 0 
Tail: 1
data = [0, 0, 0, 0, 0, 1, 1, 1]
import pymc3 as pm

with pm.Model() as coin_flipping:
    p = pm.Uniform(‘p’, lower=0, upper=1)
    y = pm.Bernoulli(‘y’, p=p, observed=data)
    trace = pm.sample(1000)
    pm.traceplot(trace)

После того, как я запустил этот код, я получил, что последним средним является E (T) = 0,398, что очень близко к результату аналитического решения (0,400).Я счастлив до сих пор, но это не окончательный ответ.Мне нужна модель, которая имитирует вероятность E ({T, T, T}).Я ценю, если кто-нибудь поможет мне на этом этапе.

1 Ответ

0 голосов
/ 15 февраля 2019

Одним из способов является сделать это эмпирически с помощью апостериорной прогнозной выборки PyMC3.То есть, если у вас есть задняя выборка, вы можете генерировать выборки из случайных параметризаций модели.Метод pymc3.sample_posterior_predictive() создаст новые выборки размером с ваши исходные данные наблюдений.Поскольку вас интересуют только три сальто, мы можем просто игнорировать дополнительные сальто, которые он генерирует.Например, если вы хотите 10000 случайных наборов предсказанных сальто, вы должны сделать 

with pm.Model() as coin_flipping:
    # this is still uniform, but I always prefer Beta for proportions
    p = pm.Beta(‘p’, alpha=1, beta=1)

    pm.Bernoulli(‘y’, p=p, observed=data)

    # chains looked a bit waggly at 1K; 10K looks smoother
    trace = pm.sample(10000, random_seed=2019, chains=4)

    # note this generates (10000, 8) observations
    post_pred = pm.sample_posterior_predictive(trace, samples=10000, random_seed=2019)

Чтобы затем увидеть, как часто следующие три сальто (1,1,1), мы можем сделать 

np.mean(np.sum(post_pred['y'][:,:3], axis=1) == 3)
# 0.0919

Аналитическое решение

В этом примере, поскольку у нас есть аналитическое апостериорное предиктивное распределение (Beta-Binomial[k | n, a=4, b=6] - см. таблицу сопряженных распределений в Википедии), мы можем точноРассчитайте вероятность наблюдения трех хвостов в следующие три сальто следующим образом: 

from scipy.special import comb, beta as beta_fn

n, k = 3, 3  # flips, tails
a, b = 4, 6  # 1 + observed tails, 1 + observed heads

comb(n, k) * beta_fn(n + a, n - k + b) / beta_fn(a, b)
# 0.09090909090909091

Обратите внимание, что beta_fn является бета-функцией Эйлера , в отличие от бета-распределения.

...