r - 正态分布均值贝叶斯推理的玩具 R 代码 [降雪量数据]

我有一些降雪观察:

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 的先验分布的时候了。原则上这是一个连续函数，但看起来我们想要一个离散的近似值，使用来自 R 包 histprior 的 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")

哇，这个好美!!

关于r - 正态分布均值贝叶斯推理的玩具 R 代码 [降雪量数据]，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/40189329/

Mu

r - 正态分布均值贝叶斯推理的玩具 R 代码 [降雪量数据]