Блог

У меня есть ряд наблюдений за снегопадом:

 x <- c(98.044, 107.696, 146.050, 102.870, 131.318, 170.434, 84.836, 154.686,
       162.814, 101.854, 103.378, 16.256)
  

и мне сказали, что это соответствует нормальному распределению с известным стандартным отклонением в 25,4, но неизвестным средним mu . Я должен сделать вывод о mu использовании байесовской формулы.

Это информация о предыдущем mu

 mean of snow |  50.8  | 76.2  | 101.6 | 127.0 |  152.4 | 177.8  
---------------------------------------------------------------
probability  |   0.1  | 0.15  | 0.25  |0.25   |  0.15  |  0.1 
---------------------------------------------------------------
  

Следующее — это то, что я пробовал до сих пор, но последняя строка о post работает некорректно. Результирующий график просто представляет собой горизонтальную линию.

 library(LearnBayes)
midpts <- c(seq(50.8, 177.8, 30))
prob <- c(0.1, 0.15, 0.25, 0.25, 0.15, 0.1)
p <- seq(50, 180, length = 40000)
histp <- histprior(p, midpts, prob)
plot(p, histp, type = "l")

# posterior density
post <- round(histp * dnorm(x, 115, 42) / sum(histp * dnorm(x, 115, 42)), 3)
plot(p, post, type = "l")
  

Мое первое предложение: убедитесь, что вы понимаете статистику, стоящую за этим. Когда я увидел ваш

 post <- round(histp * dnorm(x, 115, 42) / sum(histp * dnorm(x, 115, 42)), 3)
  

Я полагал, что вы перепутали несколько концепций. Похоже, что это формула Байеса, но у вас неверный код для определения вероятности. Правильная функция правдоподобия равна

 ## likelihood function: `L(obs | mu)`
## standard error is known (to make problem easy) at 25.4
Lik <- function (obs, mu) prod(dnorm(obs, mu, 25.4))
  

Обратите внимание, mu является неизвестным, поэтому оно должно быть переменной этой функции; кроме того, вероятность является произведением всей индивидуальной плотности вероятности при наблюдениях. Теперь мы можем оценить вероятность, например, при mu = 100 с помощью

 Lik(x, 100)
# [1] 6.884842e-30
  

Для успешной реализации R нам нужна векторизованная версия функции Lik . То есть функция, которая может вычислять на векторном входе для mu , а не просто на скалярном входе. Я буду просто использовать sapply для векторизации:

 vecLik <- function (obs, mu) sapply(mu, Lik, obs = obs)
  

Давайте попробуем

 vecLik(x, c(80, 90, 100))
# [1] 6.248416e-34 1.662366e-31 6.884842e-30
  

Теперь пришло время получить предварительное распределение для mu . В принципе, это непрерывная функция, но, похоже, мы хотим ее дискретное приближение, используя histprior из R package LearnBayes .

 ## prior distribution for `mu`: `prior(mu)`
midpts <- c(seq(50.8, 177.8, 30))
prob <- c(0.1, 0.15, 0.25, 0.25, 0.15, 0.1)
mu_grid <- seq(50, 180, length = 40000)  ## a grid of `mu` for discretization
library(LearnBayes)
prior_mu_grid <- histprior(mu_grid, midpts, prob)  ## discrete prior density
plot(mu_grid, prior_mu_grid, type = "l")
  

Прежде чем применять формулу Байе, мы сначала вычисляем нормализующую константу NC в знаменателе. Это было бы интегралом от Lik(obs | mu) * prior(mu) . Но поскольку у нас есть дискретное приближение для prior(mu) , мы используем сумму Римана для аппроксимации этого интеграла.

 delta <- mu_grid[2] - mu_grid[1]    ## division size
NC <- sum(vecLik(x, mu_grid) * prior_mu_grid * delta)    ## Riemann sum
# [1] 2.573673e-28
  

Отлично, все готово, мы можем использовать формулу Байеса:

 posterior(mu | obs) = Lik(obs | mu) * prior(mu) / NC
  

Опять же, поскольку prior(mu) дискретизировано, posterior(mu) тоже дискретизируется.

 post_mu <- vecLik(x, mu_grid) * prior_mu_grid / NC
  

Ха-ха, давайте сделаем набросок после mu , чтобы увидеть результат вывода:

 plot(mu_grid, post_mu, type = "l")
  

Вау, это красиво!!