R语言MCMC:Metropolis |
您所在的位置:网站首页 › 混沌经txt › R语言MCMC:Metropolis |
全文链接:http://tecdat.cn/?p=19664 最近我们被客户要求撰写关于Metropolis-Hastings采样的研究报告,包括一些图形和统计输出。 MCMC是从复杂概率模型中采样的通用技术。 蒙特卡洛 马尔可夫链 Metropolis-Hastings算法 问题如果需要计算有复杂后验pdf p(θ| y)的随机变量θ的函数f(θ)的平均值或期望值。 您可能需要计算后验概率分布p(θ)的最大值。 解决期望值的一种方法是从p(θ)绘制N个随机样本,当N足够大时,我们可以通过以下公式逼近期望值或最大值 将相同的策略应用于通过从p(θ| y)采样并取样本集中的最大值来找到argmaxp(θ| y)。 解决方法1.1直接模拟 1.2逆CDF 1.3拒绝/接受抽样 如果我们不知道精确/标准化的pdf或非常复杂,则MCMC会派上用场。 马尔可夫链为了模拟马尔可夫链,我们必须制定一个 过渡核T(xi,xj)。过渡核是从状态xi迁移到状态xj的概率。 马尔可夫链的收敛性意味着它具有平稳分布π。马尔可夫链的统计分布是平稳的,那么它意味着分布不会随着时间的推移而改变。 Metropolis算法对于一个Markov链是平稳的。基本上表示 处于状态x并转换为状态x'的概率必须等于处于状态x'并转换为状态x的概率 或者 方法是将转换分为两个子步骤;候选和接受拒绝。 令q(x'| x)表示 候选密度,我们可以使用概率 α(x'| x)来调整q 。 候选分布 Q(X'| X)是给定的候选X的状态X'的条件概率, 和 接受分布 α(x'| x)的条件概率接受候选的状态X'-X'。我们设计了接受概率函数,以满足详细的平衡。 该 转移概率 可以写成: 插入上一个方程式,我们有 Metropolis-Hastings算法A的选择遵循以下逻辑。 在q下从x到x'的转移太频繁了。因此,我们应该选择α(x | x')=1。但是,为了满足 细致平稳,我们有 下一步是选择满足上述条件的接受。Metropolis-Hastings是一种常见的 选择: 即,当接受度大于1时,我们总是接受,而当接受度小于1时,我们将相应地拒绝。因此,Metropolis-Hastings算法包含以下内容: 初始化:随机选择一个初始状态x; 根据q(x'| x)随机选择一个新状态x'; 3.接受根据α(x'| x)的状态。如果不接受,则不会进行转移,因此无需更新任何内容。否则,转移为x'; 4.转移到2,直到生成T状态; 5.保存状态x,执行2。 原则上,我们从分布P(x)提取保存的状态,因为步骤4保证它们是不相关的。必须根据候选分布等不同因素来选择T的值。 重要的是,尚不清楚应该使用哪种分布q(x'| x);必须针对当前的特定问题进行调整。 属性Metropolis-Hastings算法的一个有趣特性是它 仅取决于比率 是候选样本x'与先前样本xt之间的概率, 是两个方向(从xt到x',反之亦然)的候选密度之比。如果候选密度对称,则等于1。 马尔可夫链从任意初始值x0开始,并且算法运行多次迭代,直到“初始状态”被“忘记”为止。这些被丢弃的样本称为预烧(burn-in)。其余的x可接受值集代表分布P(x)中的样本 Metropolis采样一个简单的Metropolis-Hastings采样让我们看看从 伽玛分布 模拟任意形状和比例参数,使用具有Metropolis-Hastings采样算法。 下面给出了Metropolis-Hastings采样器的函数。该链初始化为零,并在每个阶段都建议使用N(a / b,a /(b * b))个候选对象。 基于正态分布且均值和方差相同gamma的Metropolis-Hastings独立采样 从某种状态开始xt。代码中的x。 在代码中提出一个新的状态x'候选 计算“接受概率” 从[0,1] 得出一些均匀分布的随机数u;如果u |t|) (Intercept) -3.1756 1.7566 -1.808 0.081 . x 5.0469 0.1964 25.697 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |