1、 多模型問題的貝葉斯模型平均 組合預測法 預測方法單一模型預測方法：灰色預測模型 多元統計 BP神經網絡法等。多模型預測方法：加權平均法 最優模型修正 最小方差法等 貝葉斯模型平均組合預測 貝葉斯模型平均組合預測克服了其他預測模型的以下缺點：第一，未考慮主觀先驗信息。第二，沒有充分提取各預測方法正確的預測信息。第三，沒有考慮模型的不確定因素。貝葉斯模型平均組合預測的關鍵是計算后驗概率，計算后驗概率的關鍵是計算邊際似然，邊際似然是一個高維、復雜的積分。目前比較好的方法是馬爾科夫蒙特卡洛法。優點：1、當用不同的條件概率分布時，不需要改變運算法則。2、全面考慮了貝葉斯模型平均的權重和方差的后驗概率。

2、貝葉斯模型平均法（BMABMA）基本表述 其中 為BMA法對系統響應的組合預測值；最后一項是各單一模型 Mk 的預測值，t為變量； Pr(Mk|D) 為給定數據D下模型的后驗概率。BMA法的實質：以各模型的后驗概率為權重，對單一模型的預測值進行加權平均得到貝葉斯模型平均值 。BMAy 蒙特卡洛算法蒙特卡洛算法基本思想：當所求解問題是某種隨機事件出現的概率，或者是某個隨機變量的期望值時，通過某種“實驗”的方法，以這種事件出現的頻率估計這一隨機事件的概率，或者得到這個隨機變量的某些數字特征，并將其作為問題的解。步驟：（1）構造或描述概率過程 （2）實現從已知概率分布抽樣 （3）建立各種估計量 馬爾

3、科夫鏈簡介馬爾可夫鏈（Markov Chain），描述了一種狀態序列，其每個狀態值取決于前面有限個狀態。馬爾可夫鏈是具有馬爾可夫性質的隨機變量的一個數列。馬爾可夫性質： MCMC MCMC簡介MCMC方法是使用馬爾科夫鏈的蒙特卡洛積分，其基本思想是：構造一條Markov鏈，使其平穩分布為待估參數的后驗分布，通過這條馬爾科夫鏈產生后驗分布的樣本，并基于馬爾科夫鏈達到平穩分布時的樣本（有效樣本）進行蒙特卡洛積分。（1）構造Markov鏈，使其收斂到平穩分布（x）。（2）產生樣本：由空間中某一點出發，用（1）的Markov鏈進行抽樣模擬，產生點序列：x1,.,xn.（3）蒙特卡洛積分。馬爾科夫鏈及其

4、平穩分布我們先來看馬氏鏈的一個具體的例子。社會學家經常把人按其經濟狀況分成3類：下層(lower-class)、中層(middle-class)、上層(upper-class)，我們用1,2,3 分別代表這三個階層。社會學家們發現決定一個人的收入階層的最重要的因素就是其父母的收入階層。如果一個人的收入屬于下層類別，那么他的孩子屬于下層收入的概率是 0.65, 屬于中層收入的概率是 0.28, 屬于上層收入的概率是 0.07。事實上，從父代到子代，收入階層的變化的轉移概率如下 使用矩陣的表示方式，轉移概率矩陣記為：假設當前這一代人處在下層、中層、上層的人的比例是概率分布向量0=0(1),0(2)

5、,0(3)，那么他們的子女的分布比例將是 1=0*P,他們的孫子代的分布比例將是 2=1*P=0*P2, , 第n代子孫的收入分布比例將是 n=n-1*P=0*Pn。 根據這種算法，我們首先分別假設初始概率分布0=0.21,0.68,0.11和0=0.75,0.15,0.1 ，則我們可以迭代計算前n代人的分布狀況如下52. 036. 012. 018. 067. 015. 007. 028. 065. 0P計算結果馬氏鏈的收斂我們發現，兩次給定不同的初始概率分布，最終都收斂到概率分布 =0.286,0.489,0.225，也就是說收斂的行為和初始概率分布 0 無關。這說明這個收斂行為主要是由概

6、率轉移矩陣P決定的。我們計算一下 PnP20=P21= =Pn=我們發現，當 n 足夠大的時候，這個Pn矩陣的每一行都是穩定地收斂到=0.286,0.489,0.225 這個概率分布。自然的，這個收斂現象并非是我們這個馬氏鏈獨有的，而是絕大多數馬氏鏈的共同行為。馬氏鏈定理： 如果一個非周期馬氏鏈具有轉移概率矩陣P,且它的任何兩個狀態是連通的，那么 存在且與i無關，記 =(j),我們有其中稱為馬氏鏈的平穩分布。馬氏鏈定理ijnnPlimijnnPlim Metropolis-Hastings Metropolis-Hastings算法算法 細致平穩條件：非周期馬氏鏈的轉移矩陣P和分布(x) 滿足

7、 (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.28-0.489*0.15=0.0067Metropolis-HastingsMetropolis-Hastings算法算法假設我們已經有一個轉移矩陣為Q馬氏鏈(q(i,j)表示從狀態 i轉移到狀態j的概率，也可以寫為 q(j|i)或者q(ij), 顯然，通常情況下 p(

8、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)取什么樣的 (i,j) 以上等式能成立呢？最簡單的，按照對稱性，我們可以取 (i,j)=p(j)q(j,i)，(j,i)=p(i)q(i,j)Metropolis-HastingsMetropolis-Hastings算法算法于是我們把原來具有轉移矩陣Q的一個很普通的馬氏鏈，改造為了具有轉移矩陣Q的馬氏鏈，而 Q恰好滿足細致平穩條

9、件，由此馬氏鏈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)的算法。Metropolis-HastingsMetropolis-Hastings算法算法上述過程中 p(x),q(x|y) 說的都是離散的情形，事實上即便這兩個分布是連續的，以上算法仍然是有效，于是就

10、得到更一般的連續概率分布 p(x)的采樣算法，而 q(x|y) 就是任意一個連續二元概率分布對應的條件分布。Metropolis-HastingsMetropolis-Hastings算法的改進算法的改進以上算法已經能滿足實際要求，但是有時候接受率很低，這樣采樣過程中馬氏鏈容易原地踏步，收斂到平穩分布p(x)的速度太慢，我們需要提升接受率來提高效率。假設 (i,j)=0.1,(j,i)=0.2, 此時滿足細致平穩條件，于是 p(i)q(i,j)0.1=p(j)q(j,i)0.2上式兩邊擴大5倍，我們改寫為 p(i)q(i,j)0.5=p(j)q(j,i)1這樣我們提高了接受率，而細致平穩條件并

11、沒有打破。這告訴我們可以將兩邊的同比例擴大并使一個較大的擴大到1。Metropolis-HastingsMetropolis-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 ac

12、cepted and the new state is set equal to the proposal, and the proposal distribution now centers on the new state.Metropolis-HastingsMetropolis-Hastings算法的簡單例子算法的簡單例子我們要用 MH 算法對標準高斯分布進行采樣，轉移函數(對稱)是方差為 0.05 的高斯矩陣matlab 代碼如下：n = 250000;x = zeros(n, 1);x(1) = 0.5;for i = 1: n-1x_c = normrnd(x(1), 0.05)

13、;if rand min(1, normpdf(x_c)/normpdf(x(i)x(i+1) = x_c;elsex(i+1) = x(i);endendMetropolis-HastingsMetropolis-Hastings算法的簡單例子算法的簡單例子% 利用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 = 1Metropolis-HastingsMetropolis-Hastings算法的簡單例子算法的簡單例子 x(k) = nextx ; x0 = nextx ; k = k + 1 ; else pp = rand(); if pp p x(k) = nextx ; x0 = nextx ; k = k + 1 ; end endendhist = histc(x(10000:k-1