игрушки R код байесовского вывода для среднего нормального распределения [данные сумм снегопада]

Question

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

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") 

Answer 1

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

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 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") 

Ничего себе, это красиво !!