关于正态分布均值 [降雪量数据] 贝叶斯推理的玩具 R 代码
Toy R code on Bayesian inference for mean of a normal distribution [data of snowfall amount]
我有一些降雪观测数据:
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
by
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 包 LearnBayes
.
中的 histprior
## 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")
哇,太美了!!
我有一些降雪观测数据:
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
by
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 包 LearnBayes
.
histprior
## 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")
哇,太美了!!