R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计

news/2024/3/29 8:46:53/文章来源:https://blog.csdn.net/tecdat/article/details/128090606

全文链接:http://tecdat.cn/?p=19664 

MCMC是从复杂概率模型中采样的通用技术。

  1. 蒙特卡洛

  2. 马尔可夫链

  3. Metropolis-Hastings算法点击文末“阅读原文”获取完整代码数据)。

问题

如果需要计算有复杂后验pdf p(θ| y)的随机变量θ的函数f(θ)的平均值或期望值。

d1c525d7ef41c4baad82d64e33a595e5.png

您可能需要计算后验概率分布p(θ)的最大值。

19c712bd5bc12cc73f55248e62bbf5e8.png

解决期望值的一种方法是从p(θ)绘制N个随机样本,当N足够大时,我们可以通过以下公式逼近期望值或最大值

59b8020c4621baf42005004875388421.png

将相同的策略应用于通过从p(θ| y)采样并取样本集中的最大值来找到argmaxp(θ| y)。

相关视频


解决方法

1.1直接模拟

1.2逆CDF

1.3拒绝/接受抽样

如果我们不知道精确/标准化的pdf或非常复杂,则MCMC会派上用场。


马尔可夫链

4ba56660c5b4e2c0072925dbb06f4e1a.png

为了模拟马尔可夫链,我们必须制定一个 过渡核T(xi,xj)。过渡核是从状态xi迁移到状态xj的概率。

 马尔可夫链的收敛性意味着它具有平稳分布π。马尔可夫链的统计分布是平稳的,那么它意味着分布不会随着时间的推移而改变。

Metropolis算法

 对于一个Markov链是平稳的。基本上表示

处于状态x并转换为状态x'的概率必须等于处于状态x'并转换为状态x的概率

9f1ffb3f268651a2c4f3efe69bf6afd7.png

或者

10cd0a2a22a87f44ee95f2eccd9e2470.png

方法是将转换分为两个子步骤;候选和接受拒绝。

令q(x'| x)表示 候选密度,我们可以使用概率 α(x'| x)来调整q  。

候选分布 Q(X'| X)是给定的候选X的状态X'的条件概率,

和 接受分布 α(x'| x)的条件概率接受候选的状态X'-X'。我们设计了接受概率函数,以满足详细的平衡。

该 转移概率 可以写成:

cd5224a25334f873a6028ddc7535fe61.png

插入上一个方程式,我们有

3d23a61e458ca1ee39669320379b7162.png

Metropolis-Hastings算法 

A的选择遵循以下逻辑。

1fc70e414ed4d9a4399a60c465f29447.png

在q下从x到x'的转移太频繁了。因此,我们应该选择α(x | x')=1。但是,为了满足 细致平稳,我们有

17d43c73676b68729dcdac808e92d85b.png

下一步是选择满足上述条件的接受。Metropolis-Hastings是一种常见的 选择:

cf4370ffd0e28259d333f45eb34a6ddb.png

即,当接受度大于1时,我们总是接受,而当接受度小于1时,我们将相应地拒绝。因此,Metropolis-Hastings算法包含以下内容:

  1. 初始化:随机选择一个初始状态x;

  2. 根据q(x'| x)随机选择一个新状态x';

3.接受根据α(x'| x)的状态。如果不接受,则不会进行转移,因此无需更新任何内容。否则,转移为x';

4.转移到2,直到生成T状态;

5.保存状态x,执行2。

原则上,我们从分布P(x)提取保存的状态,因为步骤4保证它们是不相关的。必须根据候选分布等不同因素来选择T的值。 重要的是,尚不清楚应该使用哪种分布q(x'| x);必须针对当前的特定问题进行调整。


属性

Metropolis-Hastings算法的一个有趣特性是它 仅取决于比率

bc1f05a0302056b006bb755b4cbe3644.png

是候选样本x'与先前样本xt之间的概率,

5a6c883526d189a33ff1f69e6242ef48.png

是两个方向(从xt到x',反之亦然)的候选密度之比。如果候选密度对称,则等于1。

马尔可夫链从任意初始值x0开始,并且算法运行多次迭代,直到“初始状态”被“忘记”为止。这些被丢弃的样本称为预烧(burn-in)。其余的x可接受值集代表分布P(x)中的样本


Metropolis采样

一个简单的Metropolis-Hastings采样

让我们看看从 伽玛分布 模拟任意形状和比例参数,使用具有Metropolis-Hastings采样算法。

下面给出了Metropolis-Hastings采样器的函数。该链初始化为零,并在每个阶段都建议使用N(a / b,a /(b * b))个候选对象。

cd2a535b9ac3580734b83b0dc0597515.png

基于正态分布且均值和方差相同gamma的Metropolis-Hastings独立采样

  1. 从某种状态开始xt。代码中的x。

  2. 在代码中提出一个新的状态x'候选

  3. 计算“接受概率”

    b989cf6df7c2a34067cf99a44bebe2c2.png

  4. 从[0,1] 得出一些均匀分布的随机数u;如果u <α接受该点,则设置xt + 1 = x'。否则,拒绝它并设置xt + 1 = xt。

MH可视化

set.seed(123)for (i in 2:n) {can <- rnorm(1, mu, sig)aprob <- min(1, (dgamma(can, a, b)/dgamma(x, a, b))/(dnorm(can, mu, sig)/dnorm(x, mu, sig)))u <- runif(1)if (u < aprob) x <- canvec[i] <- x

画图

设置参数。

nrep<- 54000
burnin<- 4000
shape<- 2.5
rate<-2.6

修改图,仅包含预烧期后的链

vec=vec[-(1:burnin)]
#vec=vec[burnin:length(vec)]
par(mfrow=c(2,1)) # 更改主框架,在一帧中有多少个图形
plot(ts(vec), xlab="Chain", ylab="Draws")
abline(h = mean(vec), lwd="2", col="red" )

59179db41bb06fdfbd7b5c3176855902.png


点击标题查阅往期内容

3940859f97c44a7a744e91a22889b82d.png

Python用MCMC马尔科夫链蒙特卡洛、拒绝抽样和Metropolis-Hastings采样算法

outside_default.png

左右滑动查看更多

outside_default.png

01

7a88fd3144168ccbfc7996747b360ecc.png

02

32ac5d7585ce6786f58c6f204df44955.png

03

626bc961ee95c57ca4c455ee98ae2bfe.png

04

af459a7045c1f0e1c0bcc9b9e80102e5.png

Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.007013 0.435600 0.724800 0.843300 1.133000 3.149000
var(vec[-(1:burnin)])
[1] 0.2976507

初始值

第一个样本 vec 是我们链的初始/起始值。我们可以更改它,以查看收敛是否发生了变化。

x <- 3*a/bvec[1] <- x

选择方案

如果候选密度与目标分布P(x)的形状匹配,即q(x'| xt)≈P(x')q(x'|),则该算法效果最佳。 xt)≈P(x')。如果使用正态候选密度q,则在预烧期间必须调整方差参数σ2。

通常,这是通过计算接受率来完成的,接受率是在最后N个样本的窗口中接受的候选样本的比例。

如果σ2太大,则接受率将非常低,因为候选可能落在概率密度低得多的区域中,因此a1将非常小,且链将收敛得非常慢。


示例2:回归的贝叶斯估计

Metropolis-Hastings采样用于贝叶斯估计回归模型。

201b848f72c4ffc755fff29490bcc1c6.png


设定参数


DGP和图

# 创建独立的x值,大约为零
x <- (-(Size-1)/2):((Size-1)/2)
# 根据ax + b + N(0,sd)创建相关值
y <-  trueA * x + trueB + rnorm(n=Size,mean=0,sd=trueSd)

4b1a564cf6fb81d3975ce556dfaf57df.png


正态分布拟然

pred = a*x + bsinglelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)sumll = sum(singlelikelihoods)

为什么使用对数

似然函数中概率的对数,这也是我求和所有数据点的概率(乘积的对数等于对数之和)的原因。

我们为什么要做这个?强烈建议这样做,因为许多小概率相乘的概率会变得很小。在某个阶段,计算机程序会陷入数值四舍五入或下溢问题。

因此, 当您编写概率时,请始终使用对数


示例:绘制斜率a的似然曲线

# 示例:绘制斜率a的似然曲线
plot (seq(3, 7, by=.05), slopelikelihoods , type="l")

c75597af5c37967619c9a54103ee0bba.png


先验分布

这三个参数的均匀分布和正态分布。

# 先验分布
# 更改优先级,log为True,因此这些均为log
density/likelihoodaprior = dunif(a, min=0, max=10, log = T)bprior = dnorm(b, sd = 2, log = T)sdprior = dunif(sd, min=0, max=30, log = T)

后验

先验和概率的乘积是MCMC将要处理的实际量。此函数称为后验函数。同样,这里我们使用和,因为我们使用对数。

posterior <- function(param){return (likelihood(param) + prior(param))
}

Metropolis算法

该算法是从 后验密度中采样最常见的贝叶斯统计应用之一 。

上面定义的后验。

  1. 从随机参数值开始

  2. 根据某个候选函数的概率密度,选择一个接近旧值的新参数值

  3. 以概率p(new)/ p(old)跳到这个新点,其中p是目标函数,并且p> 1也意味着跳跃

  4. 请注意,我们有一个 对称的跳跃/ 候选分布 q(x'| x)。

094608af4fe08984697c8228dc482ede.png

标准差σ是固定的。f11ea6737433084a149261930e1fb65b.png

所以接受概率等于

2f712d43d8cab239493fa76852a1dc43.png

######## Metropolis 算法 ################for (i in 1:iterations){probab = exp(posterior(proposal) - posterior(chain[i,]))if (runif(1) < probab){chain[i+1,] = proposal}else{chain[i+1,] = chain[i,]}

实施

(e)输出接受的值,并解释。

chain = metrMCMC(startvalue, 5500)
burnIn = 5000
accep = 1-mean(duplicated(chain[-(1:burnIn),]))

算法的第一步可能会因初始值而有偏差,因此通常会被丢弃来进行进一步分析(预烧期)。令人感兴趣的输出是接受率:候选多久被算法接受拒绝一次?候选函数会影响接受率:通常,候选越接近,接受率就越大。但是,非常高的接受率通常是无益的:这意味着算法在同一点上“停留”,这导致对参数空间(混合)的处理不够理想。

我们还可以更改初始值,以查看其是否更改结果/是否收敛。

startvalue = c(4,0,10)

小结

V1              V2                V3        Min.   :4.068   Min.   :-6.7072   Min.   : 6.787  1st Qu.:4.913   1st Qu.:-2.6973   1st Qu.: 9.323  Median :5.052   Median :-1.7551   Median :10.178  Mean   :5.052   Mean   :-1.7377   Mean   :10.385  3rd Qu.:5.193   3rd Qu.:-0.8134   3rd Qu.:11.166  Max.   :5.989   Max.   : 4.8425   Max.   :19.223
#比较:
summary(lm(y~x))
Call:
lm(formula = y ~ x)
Residuals:Min      1Q  Median      3Q     Max 
-22.259  -6.032  -1.718   6.955  19.892 
Coefficients:Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.1756     1.7566  -1.808    0.081 .  
x             5.0469     0.1964  25.697   <2e-16 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1
Residual standard error: 9.78 on 29 degrees of freedom
Multiple R-squared:  0.9579,    Adjusted R-squared:  0.9565 
F-statistic: 660.4 on 1 and 29 DF,  p-value: < 2.2e-16
summary(lm(y~x))$sigma
[1] 9.780494
coefficients(lm(y~x))[1]
(Intercept) -3.175555
coefficients(lm(y~x))[2]
x 
5.046873

总结:

### 总结: #######################
par(mfrow = c(2,3))
hist(chain[-(1:burnIn),1],prob=TRUE,nclass=30,col="109" 
abline(v = mean(chain[-(1:burnIn),1]), lwd="2")

e1b690459d262ad7fb592a40a517ce84.png


08c72e27db5f94a6051249eef8dbd2e9.jpeg

点击文末“阅读原文”

获取全文完整代码数据资料。

本文选自《R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计》。

点击标题查阅往期内容

Python用MCMC马尔科夫链蒙特卡洛、拒绝抽样和Metropolis-Hastings采样算法

R语言贝叶斯METROPOLIS-HASTINGS GIBBS 吉布斯采样器估计变点指数分布分析泊松过程车站等待时间

R语言马尔可夫MCMC中的METROPOLIS HASTINGS,MH算法抽样(采样)法可视化实例

python贝叶斯随机过程:马尔可夫链Markov-Chain,MC和Metropolis-Hastings,MH采样算法可视化

Python贝叶斯推断Metropolis-Hastings(M-H)MCMC采样算法的实现

Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

Matlab用BUGS马尔可夫区制转换Markov switching随机波动率模型、序列蒙特卡罗SMC、M H采样分析时间序列

R语言RSTAN MCMC:NUTS采样算法用LASSO 构建贝叶斯线性回归模型分析职业声望数据

R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机波动率SV模型、粒子滤波、Metropolis Hasting采样时间序列分析

R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型

R语言贝叶斯MCMC:用rstan建立线性回归模型分析汽车数据和可视化诊断

R语言贝叶斯MCMC:GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例

R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数

R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病

R语言中贝叶斯网络(BN)、动态贝叶斯网络、线性模型分析错颌畸形数据

R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

Python贝叶斯回归分析住房负担能力数据集

R语言实现贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析

Python用PyMC3实现贝叶斯线性回归模型

R语言用WinBUGS 软件对学术能力测验建立层次(分层)贝叶斯模型

R语言Gibbs抽样的贝叶斯简单线性回归仿真分析

R语言和STAN,JAGS:用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据

R语言基于copula的贝叶斯分层混合模型的诊断准确性研究

R语言贝叶斯线性回归和多元线性回归构建工资预测模型

R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例

R语言stan进行基于贝叶斯推断的回归模型

R语言中RStan贝叶斯层次模型分析示例

R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

R语言随机搜索变量选择SSVS估计贝叶斯向量自回归(BVAR)模型

WinBUGS对多元随机波动率模型:贝叶斯估计与模型比较

R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样

R语言贝叶斯推断与MCMC:实现Metropolis-Hastings 采样算法示例

R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化

视频:R语言中的Stan概率编程MCMC采样的贝叶斯模型

R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计

146fba30a63fa454295c44e89c6b9e4c.png

219863fb12c180bcd7639c836e8c8c89.jpeg

1e105b0f86b3b9350f1840e26185b0fe.png

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.luyixian.cn/news_show_38556.aspx

如若内容造成侵权/违法违规/事实不符,请联系dt猫网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

简单聊聊什么是react-redux,它能解决哪些问题

或许 在大多数人眼中 redux是一个相对复查很多的知识点 但确实如果你熟悉了流程 其实也比较简单的 redux是一个数据管理方案 我们先来举个例子 目前我们知道 react中有两种组件数据通信的方式 分别是 props 父传子 定义事件 子传父 通过事件将自己的数据传给父级 那如果是兄弟…

领悟《信号与系统》之 周期信号的傅里叶变换计算

周期信号的傅里叶变换计算一、周期信号的傅里叶变换存在的条件二、周期信号的傅里叶变换例题&#xff1a;一、周期信号的傅里叶变换存在的条件 典型非周期信号&#xff08;如指数信号&#xff0c;矩形信号等&#xff09;都是满足绝对可积&#xff08;或绝对可和&#xff09;条…

一种高选择性和灵敏的荧光生物标记物,可用于标记碱性磷酸酶 (ALP),5-FAM-Alkyne,510758-19-7,荧光生物标记物

【中文名称】5-羧基荧光素-炔烃【英文名称】 FAM alkyne,5-isomer&#xff0c;5-FAM alkyne【结 构 式】 【CAS号】510758-19-7 【分子式】C24H15NO6【分子量】413.39【基团】炔基基团【纯度】95%【规格标准】5mg&#xff0c;10mg&#xff0c;25mg&#xff0c;包装灵活&#x…

使用PyTorch实现简单的AlphaZero的算法(3):神经网络架构和自学习

神经网络架构和训练、自学习、棋盘对称性、Playout Cap Randomization&#xff0c;结果可视化 从我们之前的文章中&#xff0c;介绍了蒙特卡洛树搜索 (MCTS) 的工作原理以及如何使用它来获得给定棋盘状态的输出策略。我们也理解神经网络在 MCTS 中的两个主要作用&#xff1b;通…

xss-labs/level9

这一关界面感觉跟上一关很像 所以我们注入上一关的为编码的答案 javascript:alert(xss) 没能弹窗 查看源代码 他说我输入的链接不合法 我压根没有输入链接 我觉得后台应该是做了一个条件的判断 应该是要有链接才会在第二处输出点回显我们的输入 根据上面的猜测 我们构造如下…

Java.Integer.bitCount(int)源码解析

bitCount前言一、由易到难&#xff0c;头脑热身二、简单优化&#xff0c;一题多解三、分治优化四、bitCount(int)源码优化总结参考文献前言 如何求解一个二进制中1的个数&#xff1f;有常规的O(N)法&#xff0c;还有基于分治的O(logN)&#xff0c;即Java的bitCount(int)方法。…

我们为什么喜欢看疯狂科学家开飞艇?

很多人可能不是科幻迷&#xff0c;也在日常生活中接触过蒸汽朋克。为什么呢&#xff1f;很简单——蒸汽朋克几乎无处不在。相比其他科幻流派&#xff0c;蒸汽朋克可能算是最“出圈”的一种。简单地说&#xff0c;蒸汽朋克是一种科幻小说类型&#xff0c;由“蒸汽 ”(steam)和“…

《论文阅读》DeepSFM: Structure From Motion Via Deep Bundle Adjustment

留个笔记自用 DeepSFM: Structure From Motion Via Deep Bundle Adjustment 做什么 首先是最基础的&#xff0c;Structure-from-Motion&#xff08;SFM&#xff09;&#xff0c;SFM可以简单翻译成运动估计&#xff0c;是一种基于dui8序列图片进行三维重建的算法。简单来说就…

移动跨平台开发跨家选型参考建议

从 iPhone 诞生至今&#xff0c;智能手机风靡全球已将近20年&#xff0c;智能手机操作系统 iOS 和 Android 也成为当仁不让的顶流般的存在&#xff0c;而作为其背后的灵魂&#xff0c;移动应用也随着技术的发展已经越来越丰富。如果从技术层面来讲&#xff0c;移动 App 也从最开…

(1)点云库PCL学习——点云的格式、PCD文件的打开和显示

1、主要参考 (1)格式说明&#xff1a; 点云库PCL学习——点云的格式、PCD文件的打开和显示 ROS知识&#xff1a;点云文件.pcd格式_无水先生的博客-CSDN博客_pcd文件 &#xff08;2&#xff09;点云滤波&#xff0c;对nan的滤波 Python点云数据处理(三)滤波与RANSAC分割 - …

省 市 县 三级联动

大纲 一、导入省市县数据表(t_region) 二、引入jar包 三、导入所需util类&#xff08;整体框架&#xff09; 四、编写代码 1、配置数据库相关信息(数据库名、用户名、密码) config.propreties #oracle9i #driveroracle.jdbc.driver.OracleDriver #urljdbc:oracle:thin:loca…

jsp393学生宿舍管理系统mysql

两个权限 管理员和 学生 1. 学生信息管理 添加学生信息&#xff08;学生号&#xff0c;姓名 院系 班级入学日期 &#xff09;修改学生信息 学生退宿舍&#xff08;可以删除指定的学生也可以成批删除&#xff09; 2. 宿舍信息管理 宿舍的基本信息&#xff08;公寓数 宿舍…

第五届“强网”拟态防御国际精英挑战赛——特邀战队篇

第五届“强网”拟态防御国际精英挑战赛即将在南京隆重开赛&#xff01;本届大赛面向全球顶尖CTF战队&#xff0c;在创新应用场景与技术的基础上&#xff0c;拓展升级赛道&#xff0c;全面覆盖典型网络设备。大赛汇集国内外60支精英战队&#xff0c;参赛阵容、数量再创新高。 本…

科普下抖音的规则,为什么别人的内容很容易火,而我的很难?

今天给大家科普下抖音的规则&#xff0c;为什么别人的内容很容易火&#xff0c;而我的很难&#xff1f; 上一篇给大家讲了现在做抖音还来得及么&#xff1f;肯定的回答&#xff0c;一直都来得及。 既然来得及&#xff0c;那么我们怎么才能做好抖音呢&#xff1f; 在我看来&a…

5 - 2 单选题

1.下列线索二叉树中&#xff08;用虚线表示线索&#xff09;&#xff0c;符合后序线索树定义的是&#xff1a;B 后序线索二叉树的构建流程就是&#xff1a; 1.后序遍历二叉树&#xff1a;d b c a 2.第一个结点的前驱是NULL&#xff0c;即d的前驱&#xff0c;d的左孩子为NULL …

web表单(详解)

目录 1. 表单的概述 1.1 表单组成 2. 表单标记 2.1 input标记 2.2 select标记 2.3 textarea标记 3.HTML5新增标记 3.1 datalist标记 3.2 date 输入类型 3.3 color输入类型 3.4 button标记 3.5 details标记和summary标记 3.6 progress标记 3.7 meter标记 4 综合…

云原生微服务治理技术朝无代理架构的演进之路

摘要&#xff1a;本文基于对微服务治理技术从SOA, 微服务框架&#xff0c;到云原生架构的历史发展总结&#xff0c;提出了一种新的基于Javaagent技术的新一代无代理架构的服务治理技术&#xff0c;并介绍了其相关的代表性开源项目Sermant。本文分享自华为云社区《云原生微服务治…

Docker安装Redis集群失败经历汇总

在程序员的开发过程中&#xff0c;Redis可以说基本上是必不可少的缓存中间件。不管是二进制包还是docker安装Redis的文章在网上都是数不胜数。我之前自己玩Redis的时候基本不是二进制包安装就是docker安装&#xff0c;也没有尝试过集群方式。每次需要的时候&#xff0c;网上百度…

Cloud-computing 实验镜像 chinaskills_cloud_iaas.iso chinaskills_cloud_paas.iso

Cloud-computing 实验镜像 最近因新项目再次进行云计算环境的搭建&#xff0c; 找这两个镜像&#xff08; 找chinaskills_cloud_paas.iso chinaskills_cloud_iaas.iso&#xff09;颇为费劲&#xff0c;用尽九牛二虎之力总算找到了&#xff0c;该大侠还分享了诸多系统镜像和完…

Centos7搭建SVN代码控制服务器

Centos7搭建SVN代码控制服务器检查SVN是否安装创建SVN版本库配置代码库设置允许访问远程仓库的用户帐号密码设置权限控制设置SVN服务配置启动svn与停止启动SVN关闭SVN访问拉取远程仓库代码检查SVN是否安装 1、centos7系统自带SVN rpm -qa subversion2、如果没有则通过yum安装 …