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

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

0 голосов
/ 16 января 2019

При подбрасывании монеты мы бы хотели вычислить 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")

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

this

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

1 Ответ

0 голосов
/ 17 января 2019

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

Вот резюме:

  • В испытании Бернулли, когда N - общее количество испытаний - и z - общее количество успешных действий велико, а базовый параметр тета - маленький, лучше работать только в лог-пространстве и никогда не принимать экспоненциальную ,
  • Более того, поскольку функция log увеличивается, сравнение постеров журнала двух распределений аналогично сравнению апостериорных.
  • Вышеприведенная реализация была неправильной, потому что формула для вычисления доказательств не верна. 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")

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

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

...