Вычислить апостериорную вероятность с учетом вероятности, распределенной по Бернулли.

При подбрасывании монеты мы хотели бы вычислить p (theta | Data), где theta - основной параметр.

  • Априор следует бета-распределению с параметрами a и b.
  • Вероятность соответствует распределению Бернулли, которое дает нам вероятность выпадения орла.

Вот реализация кода:

    a = 1  # a and b are the beta distribution's parameters
    b= 1
    num = 1e5 #Number of candidate theta values
    z= 17220 #Number of heads
    N= 143293 #Total number of flips

    Theta = seq(0.07,0.12, length.out= num)
    prior = dbeta(Theta, a,b) #Compute the prior at each value 

    likelihood = Theta^z *(1-Theta)^(N-z)

    pData = likelihood * prior /sum(likelihood * prior) #Compute evidence
    posterior = likelihood*prior / pData

Я хотел бы убедиться, что апостериорное решение равно бета аналитических решений (a + z, N-z + b). Однако, поскольку вероятность равна 0, потому что значения тета малы, вероятность доказательства равна Nan, как и апостериорная.

Я попытался вычислить логарифмическую вероятность, но это дает мне большое отрицательное число, которое равно 0 при использовании экспоненты.

 Theta = seq(0.07,0.12, by= num_steps)
 lprior = log(dbeta(Theta, a,b)) #Compute the log prior at each value 

 llikelihood = log(Theta)*z + log(1-Theta)*(N-z) #log likelihood

 lpData = llikelihood + lprior - sum(llikelihood + lprior) #compute evidence

  lposterior = llikelihood+lprior - lpData
  posterior = exp(lposterior)
  plot(Theta, posterior, type="l")
  lines(Theta, exp(llikelihood), type="l")
  lines(Theta, exp(lprior), type="l")

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

это

Как мне вычислить каждое значение?


r statistics bayesian distribution coin-flipping
person Salma Bouzid    schedule 16.01.2019    source источник
comment
Вы определенно захотите поработать с вероятностями журнала. Вы могли бы проверить это log_posterior = log(analytic_solution). Если возвращение журналов к исходному масштабу вызывает проблемы, вы можете просто сделать все в масштабе журнала.   -  person Joseph Clark McIntyre    schedule 16.01.2019
comment
Ваш априор имеет одинаковое значение (1) для каждой точки тэты. Это намеренно?   -  person Jonny Phelps    schedule 16.01.2019
comment
Я предполагаю, что информативный априор является распространяемым бета-версией, поэтому dbeta всегда будет давать значения, равные 1. Как вы думаете, мне следует вычислять его по-другому?   -  person Salma Bouzid    schedule 16.01.2019
comment
@JosephClarkMcIntyre Я отредактировал сообщение, чтобы включить вычисление логарифмической шкалы. Я все еще нахожу нули в конце. Я уверен, что в моих вычислениях есть ошибка.   -  person Salma Bouzid    schedule 16.01.2019
comment
Но вы возводите в степень в конце с posterior = exp(lposterior), и вы сказали, что именно здесь возникают проблемы. Так что просто работайте с lposterior и сравните его с журналом аналитического решения.   -  person Joseph Clark McIntyre    schedule 16.01.2019
comment
@JosephClarkMcIntyre: dbeta (Theta, z + a, N-z + b) дает мне 1, следовательно, журнал дает мне 0. Считаете ли вы, что использование бета-априорной вероятности / вероятности Бернулли не может работать для малых значений тета? Если моя конечная цель - построить апостериорную, вероятностную и априорную. Как мне вычислить каждое из этих значений, чтобы везде не было нулей?   -  person Salma Bouzid    schedule 16.01.2019
comment
Вы пробовали dbeta(Theta, z+a, N-z+b, log = TRUE)? Иногда не хватает точности в способе хранения чисел, чтобы сначала получить PDF-файл, а затем зарегистрировать его, но R вернет журнал напрямую (я чувствую, что вы должны получить dbeta(Theta, z+a, N-z+b) = 0, но 'dbeta (Theta, z + a, N -z + b, log = TRUE) = - (какое-то большое число) '; когда я использовал Theta = 0,07, N = 143293, z = 17220, a = 1, b = 1, я получил -2308)   -  person Joseph Clark McIntyre    schedule 16.01.2019
comment
Вы правы, log = True решает проблему аналитической формулы, но не формулы Бернулли: llikelihood = log (Theta) * z + log (1-Theta) * (N-z) #log likelihood. Нули везде .... И как мне нанести на график априорную и апостериорную вероятность?   -  person Salma Bouzid    schedule 16.01.2019
comment
@JosephClarkMcIntyre, может быть, что-то не так с моим расчетом априорной вероятности? Почему я не могу получить график Бернулли, который напоминал бы нормальные распределения, поскольку у меня много образцов?   -  person Salma Bouzid    schedule 16.01.2019
comment
Возможно, попробуйте log1p(-Theta) вместо log(1-Theta), если это проблема числовой стабильности.   -  person merv    schedule 17.01.2019
comment
Я их совсем не понимаю. Когда я пробую те же тета, N и z, как я упоминал выше, я получаю -54941,66. Учитывая используемые вами значения тета, журнал всегда будет отрицательным, и вы никогда не должны получать число, близкое к 0 (если вы не возведете в степень, что я считаю ошибкой).   -  person Joseph Clark McIntyre    schedule 17.01.2019
comment
Моя конечная цель - фактически воспроизвести результаты библиотеки (bayesAB) для A / B-тестирования. Чтобы вычислить апостериорную обработку и контроль, я хотел сделать то же самое, что и в сценарии подбрасывания монеты. Вы хоть представляете, как вычисляется апостериор в этом пакете? library(bayesAB) bayes.test <- bayesTest(control$converted,treatment$converted, priors = c('alpha' = 1, 'beta' = 1), n_samples = 1e5, distribution = 'bernoulli') print(bayes.test)   -  person Salma Bouzid    schedule 17.01.2019
comment
Думаю, я понял. Я скоро отправлю ответ. Спасибо @JosephClarkMcIntyre за ваши подробные разъяснения и merv за ваши правки и трюк с числовой стабильностью :)   -  person Salma Bouzid    schedule 17.01.2019

Ответы (1)


arrow_upward
0
arrow_downward

Этот ответ был предоставлен в разделе комментариев @JosephClarkMcIntyre.

Вот краткое изложение:

  • В испытании Бернулли, когда N - общее количество попыток - и z - общее количество успехов велико, а базовый параметр тета мал, лучше работать только в пространстве журнала и никогда не брать экспоненту.
  • Более того, поскольку логарифмическая функция увеличивается, сравнение логарифмических апостериорных значений двух распределений аналогично сравнению апостериорных значений.
  • Вышеупомянутая реализация была неправильной, потому что формула для вычисления доказательства неверна. p (свидетельство) = сумма (вероятность * априор), p (log_evidence) = сумма (log_likelihood + log_prior)

Это последний код, в котором в логе находятся априор, вероятность и свидетельство:

  a = 1  # a and b are the beta distribution's paramteres
  b= 1
  num_steps = 1e5
  z= 17220 #Number of heads
  N= 143293 #Total number of flips

  Theta = seq(from=0.07,to=0.12, length.out= num_steps)
  lprior = dbeta(Theta, a,b,log=TRUE) #Compute the log prior at each value 

  llikelihood = log(Theta)*z + log1p(-Theta)*(N-z) #log likelihood

  lpData = sum(llikelihood + lprior) #compute log of the evidence

  lposterior = llikelihood+lprior - lpData
  plot(Theta,log(dbeta(Theta,a+z,N-z+b)))
  plot(Theta, lposterior, type="l")

Однако аналитический и вычисленный апостериорный лог не совпадают, как показано на графике.

Не стесняйтесь комментировать, если вы считаете, что в этом ответе есть недостаток, или объясните, почему аналитический и вычисленный апостериорный лог не совпадают. ^^

person Salma Bouzid    schedule 17.01.2019