多模型贝叶斯平均法及MCMC

上式两边扩大5倍，我们改写为
p(i)q(i,j)×0.5=p(j)q(j,i)×1
这样我们提高了接受率，而细致平稳条件并没有打破。这 告诉我们可以将两边的α同比例扩大并使一个较大的α扩大 到1。
Metropolis-Hastings算法的改进
illustration of the Metropolis sampler to sample from target density. (A)the current state of the chain is θ(t). (B) a proposal distribution around the current state is used to generate a proposal θ* (C) the proposal was accepted and the new state is set equal to the proposal, and the proposal distribution now centers on the new state.

以上算法已经能满足实际要求，但是有时候接受率α很低， 这样采样过程中马氏链容易原地踏步，收敛到平稳分布 p(x)的速度太慢，我们需要提升接受率来提高效率。 假设 α(i,j)=0.1,α(j,i)=0.2, 此时满足细致平稳条件，于是 p(i)q(i,j)×0.1=p(j)q(j,i)×0.2


Metropolis-Hastings算法

假设我们已经有一个转移矩阵为Q马氏链(q(i,j)表示从状 态 i转移到状态j的概率，也可以写为 q(j|i)或者q(i→j)), 显 然，通常情况下 p(i)q(i,j)≠p(j)q(j,i) 也就是细致平稳条件不成立，所以 p(x) 不太可能是这个 马氏链的平稳分布。我们可否对马氏链做一个改造，使 得细致平稳条件成立呢？譬如，我们引入一个 α(i,j), 我们 希望 p(i)q(i,j)α(i,j)=p(j)q(j,i)α(j,i)

x(k) = nextx ; x0 = nextx ;


k=k+1;
else pp = rand(); if pp < p x(k) = nextx ; x0 = nextx ; k=k+1; end end end
hist = histc(x(10000:k-1),0:0.1:100);
多模型问题的贝叶斯模型平均 组合预测法
预测方法

单一模型预测方法：灰色预测模型 多元统计 BP神经网络法等。


多模型预测方法：加权平均法
最优模型修正 最小方差法等
贝叶斯模型平均组合预测

贝叶斯模型平均组合预测克服了其他预测模型的以下缺 点：


第一，未考虑主观先验信息。
第二，没有充分提取各预测方法正确的预测信息。 第三，没有考虑模型的不确定因素。 贝叶斯模型平均组合预测的关键是计算后验概率，计算 后验概率的关键是计算边际似然，边际似然是一个高维、 复杂的积分。目前比较好的方法是马尔科夫蒙特卡洛法。 优点：1、当用不同的条件概率分布时，不需要改变运算 法则。2、全面考虑了贝叶斯模型平均的权重和方差的后 验概率。
plot(hist);
Metropolis-Hastings算法的简单例子
Metropolis-Hastings算法的简单例子

xx = 0:0.1:100; n=5; a=4; C=1; y = C .* xx.^(-n/2).*exp(a./(2*xx)); figure; plot(y);


马氏链定理

马氏链定理： 如果一个非周期马氏链具有转移概率矩阵P, 且它的任何两个状态是连通的，那么 lim P n ij 存在且与i无 n 关，记 lim P n ij =π(j),我们有
n

其中π称为马氏链的平稳分布。
Metropolis-Hastings算法

细致平稳条件：非周期马氏链的转移矩阵P和分布π(x) 满 足
使用矩阵的表示方式，转移概率矩阵记为：
0.65 0.28 0.07 P 0 . 15 0 . 67 0 . 18 0.12 0.36 0.52
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 π0=[π0(1),π0(2),π0(3)]，那么他们的子女的分布比例将是 π1=π0*P,他们的孙 子代的分布比例将是 π2=π1*P=π0*P^2, ……, 第n代子孙的收入分布比例将是 πn=πn-1*P=π0*P^n。 根据这种算法，我们首先分别假设初始概率分布π0=[0.21,0.68,0.11]和 π0=[0.75,0.15,0.1] ，则我们可以迭代计算前n代人的分布状况如下


取什么样的 α(i,j) 以上等式能成立呢？最简单的，按照对 称性，我们可以取
α(i,j)=p(j)q(j,i)，α(j,i)=p(i)q(i,j)

Metropolis-Hastings算法

于是我们把原来具有转移矩阵Q的一个很普通的马氏链， 改造为了具有转移矩阵Q′的马氏链，而 Q′恰好满足细 致平稳条件，由此马氏链Q′的平稳分布就是p(x)。 在改造 Q 的过程中引入的 α(i,j)称为接受率，物理意义可 以理解为在原来的马氏链上，从状态 i 以q(i,j) 的概率转跳 转到状态j 的时候，我们以α(i,j)的概率接受这个转移，于 是得到新的马氏链Q′的转移概率为q(i,j)α(i,j)。 假设我们已经有一个转移矩阵Q(对应元素为q(i,j)), 把以上 的过程整理一下，我们就得到了如下的用于采样概率分 布p(x)的算法。
计算结果
马氏链的收敛

我们发现，两次给定不同的初始概率分布，最终都收敛到概 率分布 π=*0.286,0.489,0.225+，也就是说收敛的行为和初始 概率分布 π0 无关。这说明这个收敛行为主要是由概率转移 矩阵P决定的。我们计算一下 P^n P^20=P^21=⋯=P^n= 我们发现，当 n 足够大的时候，这个P^n矩阵的每一行都是 稳定地收敛到π=*0.286,0.489,0.225+ 这个概率分布。自然的， 这个收敛现象并非是我们这个马氏链独有的，而是绝大多数 马氏链的共同行为。

蒙特卡洛算法基本思想：当所求解问题是某种随机事件 出现的概率，或者是某个随机变量的期望值时，通过某 种“实验”的方法，以这种事件出现的频率估计这一随 机事件的概率，或者得到这个随机变量的某些数字特征， 并将其作为问题的解。 步骤：（1）构造或描述概率过程 （2）实现从已知概率分布抽样 （3）建立各种估计量
%% 利用Metropolis-Hastings 算法产生样本 % 这里的样本分布为：p(x) = C* x.^(-n/2)*exp(-a/(2x)); % 对于参数n=5，参数a=4.
xlen = 200000;
x = zeros(1,xlen); x0 = 1 ; len = length(x); k=1; while k <= len nextx = 100*rand();
Metropolis-Hastings算法的简单例子

我们要用 MH 算法对标准高斯分布进行采样，转移函数(对称)是方差 为 0.05 的高斯矩阵
matlab 代码如下：
n = 250000; x = zeros(n, 1); x(1) = 0.5;


for i = 1: n-1


马尔科夫链及其平稳分布

我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状 况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)， 我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收 入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属 于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收 入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子 代，收入阶层的变化的转移概率如下


贝叶斯模型平均法（BMA）

基本表述

其中yBMA为BMA法对系统响应的组合预测值；
最后一项是各单一模型 Mk 的预测值，t为变量； Pr(Mk|D) 为给定数据D下模型的后验概率。 BMA法的实质：以各模型的后验概率为权重，对单一模型的预测值进行 加权平均得到贝叶斯模型平均值 。



蒙特卡洛算法
x_c = normrnd(x(1), 0.05); if rand < min(1, normpdf(x_c)/normpdf(x(i)))
x(i+1) = x_c;
else x(i+1) = x(i); end end
Metropolis-Hastings算法的简单例子



%compute the p from x0
pn = (nextx^(-2.5))*(exp(-2/(nextx))); p0 = (x0^(-2.5))*(exp(-2/(x0))); p = min(pn/p0,1); if p >= 1
Metropolis-Hastings算法的简单例子


π(i)*Pij=π(j)*Pji (for all i,j)
细致平稳条件的物理含义就是对于任何两个状态i,j, 从 i转 移出去到j 而丢失的概率质量，恰好会被从 j 转移回i 的概 率质量补充回来，所以状态i上的概率质量π(i)是稳定的， 从而π(x)是马氏链的平稳分布。 简单验证：m=π(1)*P12-π(2)*P21 =0.286*0.280.489*0.15=0.0067