网站备案照相怎么照,百度广告投放公司,wordpress 标签别名 id,工装公司是做什么的文章目录 贝叶斯简介Why R理论基础一、三种先验分布和对应后验的计算1. 离散先验2.Beta先验#xff08;共轭先验#xff09;3. 直方图先验 二. 后验抽样1. 网格点采样法2. 其他方法 三、贝叶斯推断1. 参数估计(1) 后验均值(2) 后验方差(3) 后验区间 2. 假设检验3. 预测(1) 先… 文章目录 贝叶斯简介Why R理论基础一、三种先验分布和对应后验的计算1. 离散先验2.Beta先验共轭先验3. 直方图先验 二. 后验抽样1. 网格点采样法2. 其他方法 三、贝叶斯推断1. 参数估计(1) 后验均值(2) 后验方差(3) 后验区间 2. 假设检验3. 预测(1) 先验预测(2) 后验预测 贝叶斯简介
贝叶斯学派和频率学派的核心区别用信念度而不是长期频率来衡量某件事发生的概率这个信念度可能会随着新信息的加入而发生变化。贝叶斯的三个核心概念
先验概率prior在没有进行任何尝试时一个人对某件事的信念例如人群中X病毒携带率为3%似然likelihood观测到的数据例如A检测结果为阳性后验posterior通过观测到的数据对某件事的信念进行修正例如A携带X病毒的概率为75%。
贝叶斯法则根据数据更新你关于模型的概率的公式 ( M O D E L ∣ D A T A ) ∝ ( M O D E L ) × ( D A T A ∣ M O D E L ) (MODEL ∣ DATA) ∝ (MODEL) × (DATA ∣ MODEL) P(MODEL∣DATA)∝P(MODEL)×P(DATA∣MODEL) Why R
众多贝叶斯计算引擎都可以在R中直接调用。
贝叶斯计算引擎 贝叶斯引擎是用于贝叶斯统计建模和推断的工具这些引擎能帮助用户构建复杂的概率模型并利用马尔可夫链蒙特卡罗MCMC等方法进行推断。
BUGSWinBUGS/OpenBUGSJAGS后续会更新相关文章专门介绍JAGSStan
用于贝叶斯分析/计算的R包
MCMCpack: 包含针对常用模型的 MCMC 算法通常在社会和行为科学中使用coda: 包含一套用于收敛诊断和输出分析的函数可用于总结、绘图和诊断 MCMC 样本的收敛性rjags: 实现在 R 中完全交互式地使用JAGS的 R 包r2jags: 基于 rjags但提供了更为简便的接口简化了 JAGS 的调用流程runjags: 一个包含许多 JAGS 增强功能的 R 包包括自定义 JAGS 模块rstan: Stan 的 R 语言接口。它允许用户方便地从 R 中拟合 Stan 模型并访问输出结果Nimble: 一个提供通用 MCMC 系统的 R 包允许对用 BUGS/JAGS 模型语言编写的模型进行自定义 MCMC。用户可以选择采样器并编写新的采样器。模型和采样器通过生成的C 代码自动编译以提高速度LaplacesDemon: 一个旨在提供完整贝叶斯环境的 R 包包括众多 MCMC算法、具有多个优化算法的拉普拉斯近似、大量示例、数十种额外的概率分布以及众多 MCMC 诊断工具等 理论基础
单参数贝叶斯模型仅包含一个待估计参数的概率模型的R语言建模实现一个直观的例子总体比例p的贝叶斯推断希望从样本数据中得到p的点估计和区间估计不确定性。本文将从总体比例p估计问题出发探讨单参数贝叶斯建模过程。
贝叶斯分析的三步走
指定先验分布写出似然函数计算参数的后验分布
上面三个步骤中我们需要研究的问题
先验分布如何合理表达先验信念
选择什么先验分布方便推导无信息先验、层级先验、共轭先验等
似然函数如何用分布描述数据生成过程
使用什么分布建模正态、二项、泊松等
后验分布如何计算后验分布和后验预测
是否能够解析计算有共轭形式时可以解析否则需数值方法如 MCMC
对于总体比例估计问题我们的似然函数就是二项分布即独立重复试验中成功次数。若我们观察到y 次成功在 n 次试验中 ∼ B i n o m i a l ( , ) ∼Binomial(,) y∼Binomial(n,p)其似然函数为 ( ) ( ∣ ) ( , ) ( 1 − ) − ()(∣)(,)^ (1−)^{−} L(p)P(y∣p)(n,y)py(1−p)n−y
贝叶斯模型 ( , ) ( ) ( ∣ ) (, ) ()(|) p(,y)()p(y∣) ( ) () ()是先验分布 ( ∣ ) (|) p(y∣)是似然函数贝叶斯模型即和y的联合分布 ( , ) (, ) p(,y)
我们关心的是参数的后验分布 ( ∣ ) ( ∣ ) (∣y)它可以由贝叶斯公式计算出 ( ∣ ) ( , ) / ( ) ( ) ( ∣ ) / ( ) ( ∣ ) (, )/() ()( ∣ )/() (∣y)p(,y)/p(y)()p(y∣)/p(y) 其中() 不依赖于当给定观测数据时可视作常数从而可以看作是一个归一化常数用于将( ∣ )化成有效的概率密度函数。 一、三种先验分布和对应后验的计算
构建先验分布的常见方法是假设先验分布属于特定的参数密度函数族。然后选择先验分布的参数称为超参数 hyper parameters使先验分布尽可能地表示先验信念。
在研究总体概率p时有三种典型先验1. 离散先验2. Beta先验共轭先验3. 直方图先验。本文下面将详细讲解这三种先验以及对应后验的计算。
1. 离散先验
假设只能取有限个值如 0.1, 0.2, …, 0.9我们给这些值一个先验概率分布 P ( p p i ) w i P(pp_i)w_i P(ppi)wi w i w_i wi求和为1。
离散先验分布
p seq(0.05, 0.95, by 0.1)
prior c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior prior/sum(prior)
plot(p, prior, type h, ylabPrior Probability)p为概率值向量prior为p对应的先验概率向量归一化处理type h为直方图选项用于绘制概率质量函数pmf 输出
已知观测数据后计算后验分布
library(LearnBayes)
data c(11, 16)
post pdisc(p, prior, data)
round(cbind(p, prior, post),2)data为观测数据11次成功16次失败pdisc用于计算离散先验下的后验概率三个参数分别是概率点向量、先验概率、观测数据 输出
对比先验和后验
library(lattice)
PRIORdata.frame(prior,p,prior)
POSTdata.frame(posterior,p,post)
names(PRIOR)c(Type,P,Probability)
names(POST)c(Type,P,Probability)
datarbind(PRIOR,POST)
xyplot(Probability~P|Type,datadata, layoutc(1,2), typeh, lwd3,colblack)data.frame用于创建一个三列的数据框names用于给数据框命名 rbind表示将两个数据框按行合并 xyplot是lattice包的核心函数用于绘制散点图或线图这里Probability~P|Type表示纵轴为Probability,横轴为P 按照Type分面板绘制 输出
2.Beta先验共轭先验 从上图介绍我们可以得到为什么要选择Beta先验它和似然函数二项分布有相同形式的密度函数这样似然*先验得到的后验还是这个形式。
问题如何确定Beta分布的两个参数
策略1: 画出来好多条密度函数找一条满足期望的策略2: 根据想要的均值/样本量确定 - B e t a ( α , β ) Beta(\alpha,\beta) Beta(α,β)分布的均值为 α / ( α β ) \alpha/(\alpha\beta) α/(αβ) - B e t a ( α , β ) Beta(\alpha,\beta) Beta(α,β)分布的等效先验样本量为 α β − 2 \alpha\beta-2 αβ−2其中 α − 1 \alpha-1 α−1次成功 β − 1 \beta-1 β−1次失败
在R中可以用LearnBayes包中的beta.select()函数找到与先验知识对分位数的认知匹配的两个超参数。
quantile1list(p.5,x.3) # 中位数
quantile2list(p.9,x.5) # 90%分位数
beta.select(quantile1,quantile2)输出 [1] 3.26 7.19 先验信息研究者认为比例 的中位数为 0.390% 分位数为 0.5选出的 α 3.26 \alpha3.26 α3.26 β 7.19 \beta7.19 β7.19 假设我们选出了一个合适的先验 B e t a ( α , β ) Beta(\alpha,\beta) Beta(α,β)观测 n n n次中有 y y y次成功就可以用贝叶斯公式求出后验 π ∝ p y ( 1 − p ) n − y p α − 1 ( 1 − p ) β − 1 p α y − 1 ( 1 − p ) β n − y − 1 \pi ∝ p^y(1-p)^{n-y}p^{\alpha-1}(1-p)^{\beta-1}p^{\alphay-1}(1-p)^{\betan-y-1} π∝py(1−p)n−ypα−1(1−p)β−1pαy−1(1−p)βn−y−1
这是 B e t a ( α y , β n − y ) Beta(\alphay,\betan-y) Beta(αy,βn−y)的核这就是共轭性的含义先验分布和后验分布来自同一种分布类型。
3. 直方图先验
就是每个区间内指定一个先验概率
midpt seq(0.05, 0.95, by 0.1) # 按照0.1的间隔划分0.05到0.95的区间
prior c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0) # 每个区间的先验概率
prior prior/sum(prior) # 归一化
curve(histprior(x,midpt,prior),from0, to1,ylabPrior density,ylimc(0,.3))输出 计算后验分布一样乘上似然函数的等价beta分布表示
curve(histprior(x,midpt,prior) * dbeta(x,s1,f1),from0, to1,ylabPosterior density)输出 二. 后验抽样
后验抽样是指在获取参数后验分布后从后验分布中采样参数进而用于后续贝叶斯推断。
1. 网格点采样法
网格点采样法也称为网格近似法 (Grid Approximation)是一种通过在参数空间中定义一个离散的网格然后计算每个网格点上的未归一化后验概率即似然函数值乘以先验概率密度值来近似后验分布的方法。一旦计算出这些概率就可以将它们归一化使其总和为 1从而得到一个离散化的后验概率分布从这个离散分布中抽取样本。
p seq(0, 1, length500) # 网格点500个0到1之间均匀的数值点
post histprior(p, midpt, prior)*dbeta(p, s1, f1) # 以直方图先验为例计算每个网格点的后验概率
post post/sum(post) # 后验概率归一化
ps sample(p, 5000, replace TRUE,prob post) # 从post分布冲抽取5000个样本
hist(ps, breaks 30, xlabp,main)输出
2. 其他方法
后续会专门写一篇文章讲各种MCMC方法。 三、贝叶斯推断
通过确定先验分布后得到的后验分布包含了关于未知参数的所有当前信息所有贝叶斯推断都基于后验分布包括
参数估计 点估计后验均值 (posterior mean)、众数 (posterior mode / MAP)、中位数 (posterior median)离散程度后验方差可信区间(credible interval)后验分位数 假设检验后验概率预测估计潜在可观测但当前未观测的量模型选择/变量选择Bayes factor,DIC
1. 参数估计
(1) 后验均值
后验分布的均值常用作 p p p点估计 B e t a ( α y , β n − y ) Beta(\alphay,\betan-y) Beta(αy,βn−y)的均值为 α y / α β n {\alphay}/{\alpha\betan} αy/αβn
(2) 后验方差
后验方差是后验分布离散程度的一种度量后验方差越大我们对参数的不确定性就越大。 B e t a ( α , β ) Beta(\alpha,\beta) Beta(α,β)的方差为 α β / ( α β ) 2 ( α β 1 ) \alpha\beta/(\alpha\beta)^2(\alpha\beta1) αβ/(αβ)2(αβ1) 对于 B e t a ( α , β ) Beta(\alpha,\beta) Beta(α,β)先验和二项似然n 次试验中有 y 次成功后验分布为 B e t a ( α y , β n − y ) Beta(\alphay,\betan-y) Beta(αy,βn−y)从而后验方差为 ( α y ) ( β n − y ) / ( α β n ) 2 ( α β n 1 ) (\alphay)(\betan-y)/(\alpha\betan)^2(\alpha\betan1) (αy)(βn−y)/(αβn)2(αβn1) 注意对于均匀先验和二项似然后验方差几乎总是小于先验方差 (3) 后验区间
两种类型的可信区间
等尾可信区间(equal-tail credible interval):
100(1 − α)% 的等尾可信区间是从后验分布的 α/2 分位数到 (1 − α/2) 分位数的区间若要构建 95% 的等尾可信区间我们需要计算后验分布的 0.025 和 0.975 分位数。
例如求后验均值的95%置信区间有两种方法一是直接根据beta分布来计算
a 3.26; b 7.19 # 选出的先验参数
s 11; f 16 # 观测数据似然
qbeta(c(0.05, 0.95), a s, b f)输出[1] 0.2555267 0.5133608
二是用模拟的方法
ps rbeta(1000, a s, b f)
quantile(ps, c(0.05, 0.95))输出5% 95% 0.2490099 0.5152768 beta表示从beta分布中生成指定个随机数这里生成了1000个样本然后计算样本的0.05和0.95分位数 对于贝叶斯学派95% 可信区间的解释是观察到观测数据后真实参数 p 落在该区间内的概率为 0.95
最高后验密度区间(highest posterior density, HPD):
区间内任意点的密度大于区间外任意点的密度在包含指定概率的所有区间中长度最短有时称为最小长度可信区间。 当后验分布高度偏斜或多峰时优于等尾可信区间 计算可以使用 TeachingDemos 包中的 R 函数hpd()或emp.hpd()
2. 假设检验
假设我们要检验关于参数 p 的以下假设 H 0 : p ≤ 0.10 v . s H 1 : p 0.10 H0 : p ≤ 0.10 \ \ \ \ v.s \ \ \ \ H1 : p 0.10 H0:p≤0.10 v.s H1:p0.10 我们只需计算后验分布在这两个区间上的概率使用R函数 pbeta即可。
3. 预测
(1) 先验预测 以离散先验为例
library(LearnBayes)
p seq(0.05, 0.95, by.1)
prior c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.8, 0.3, 0.1, 0.05)
prior prior/sum(prior)
ns 20
ys 0:20
pred pdiscp(p, prior, ns, ys)
round(cbind(0:20, pred), 3)前面都和第一部分的离散先验一样ns表示未来样本量ys表示未来可能的成功个数从020pdiscp用于计算离散先验下的后验预测分布四个参数分别为概率点向量、先验概率、未来样本量、未来成功次数 输出
(2) 后验预测
以Beta先验为例有解析计算和模拟两种方法
解析计算本质就是更新beta分布的两个参数
library(LearnBayes)
ab c(3.26, 7.19)
ns 20; ys 0:20
pred pbetap(ab, ns, ys)
round(cbind(0:20, pred), 3)输出 pbetap函数计算 Beta 先验下的后验预测分布Beta-Binomial 分布三个参数分别是Beta先验的参数向量、未来试验次数、需要预测的成功次数 2.基于模拟模拟的方法适合任何先验分布下例展示了beta先验的后验模拟
m 1000
a 3.26; b 7.19 # 先验Beta分布的参数
s 11; f 16
ns 20
p rbeta(m, as, bf)
y rbinom(m, ns, p)
freq table(y) # 统计各成功次数y出现的频次
ys as.integer(names(freq)) # 提取成功次数去重
predprob freq/sum(freq) # y的概率分布
# 成功概率分布的直方图绘制
plot(ys, predprob, typeh, xlaby,
ylabPredictive Probability,
cex.lab1, cex.axis1.0)
# 解析解验证
library(LearnBayes)
TruePred1 pbetap(c(as,bf), ns, ys) # beta-二项后验分布的概率
points(ys, TruePred1, typep, colred)rbeta从后验分布中生成1000个p的随机样本rbinom对每个采样的 p生成 ns20 次试验的成功次数 y 输出 如果后验是混合分布那么后验预测分布是一个连续混合分布或复合分布下例中展示了用模拟的方法从混合正态分布0.3N(0, 1) 0.7N(3, 1/2)中随机抽样
n - 10000; alpha - 0.3
k - sample(1:2, size n, replace T, probc(alpha, 1-alpha)) # 根据概率分到1组或2组
mu - numeric(length(k))
sd - numeric(length(k))
mu - c(0,3)[k] # 根据k选择均值0或3
sd - c(1, 1/2)[k]# 根据k选择方差1或1/2
x - rnorm(n, mu, sd) # 生成混合样本
hist(x, breaksFD, probT, xlimc(-4,6), ylimc(0,0.8), mainHistogram of a mixture normal)
# 混合密度函数
p - function(x,alpha){
alpha*dnorm(x, 0, 1)(1-alpha)*dnorm(x, 3, 1/2)
}
curve(dnorm(x,0,1), col4, lty2, lwd2, addT) # N01
curve(dnorm(x,3,1/2), col4, lty4, lwd2, addT) # N31/2
curve(p(x,0.3), col2, lty1, lwd 2, addT) # 混合密度
legend(topleft, legend c(N(0,1),N(1,1/2), 0.3N(0,1)0.7N(1,1/2)), cex0.6, lty c(2,4,1), colc(4, 4, 2), lwd c(2,2,2))输出
