统计计算 课件 4.2 Metropolis-Hasting算法

4.2Metropolis-Hasting算法(梅特罗波利斯－黑斯廷斯)4.2.1Metropolis-Hasting算法4.2.2模拟退火算法

马尔可夫链蒙特卡洛法蒙特卡洛法（MonteCarlomethod)，是通过从概率模型的随机抽样进行近似数值计算的方法。马尔可夫链蒙特卡洛法是以马尔可夫链（Markovchain)为概率模型的蒙特卡罗法。马尔可夫链蒙特卡洛法构建一个马尔可夫链，使其平稳分布就是要进行抽样的分布，首先基于该马尔可夫链进行随机游走，产生样本的序列，之后使用该平稳分布的样本进行近似数值计算Metropolis-Hastings算法是最基本的马尔可夫链蒙特卡罗法吉布斯抽样（Gibbssampling）是更简单、使用更广泛的马尔可夫链蒙特卡罗法，马尔可夫链蒙特卡罗法被应用于概率分布的估计、定积分的近似计算、最优化问题的近似求解等问题，特别是被应用于统计学习中概率模型的学习与推理，是重要的统计学习计算方法。M-H采样算法是1953年由美国物理学家尼古拉斯·梅特罗波利斯与加拿大统计学家W·K·黑斯廷斯提出的，这是第一个基于马尔科夫链的蒙特卡洛方法，也被评选为20世纪最重要的算法之一。该算法的思想是基于接受-拒绝抽样方法（舍选法）。在接受-拒绝抽样方法中，有建议概率密度函数g、目标概率密度函数f和一个概率函数（接受函数）h。Metropolis算法假设在样本空间上存在目标(采样的分布)概率密度函数f，需要在Ω上产生样本的马尔科夫链，使得它的平稳分布恰好是该采样分布。4.2.1Metropolis—Hasting采样模拟退火方法来源于统计力学。它包含两个部分：退火部分和Metropolis，分别对应于外循环和内循环。内循环指的是在每次温度下迭代数次，寻找在该温度下能量的最小值，即使内部能量最小。外循环指的是先将固体达到一个较高的温度（初始温度），然后使温度按照一定的比例下降，一直到温度变为终止温度停止。退火过程可以使用蒙特卡洛方法进行模拟。它以某一高温为起点x0（初始温度），根据x0使用某种算法求出候选解x*，利用预先设计的算法得到一个新的状态。将新状态的能量和旧状态的能量对比，如果新状态的能量小于旧状态的能量（能量减少），则接受新状态，这就是下一步的状态。如果新状态的能量大于旧状态的能量（能量变大，偏离全局最优位置，与目标不符），进行概率判断，生成一个[0,1]上的随机数u，如果u小于接受概率α，则接受新状态，作为本步的状态，否则拒绝状态转移，使用前面的状态。当状态稳定之后，达到了目前状态的最优解，即平稳分布。可定义系统的状态转移的接受概率为平稳分布充分条件

detailedbalanceij构造转移矩阵q(i,j)ijAcceptionRate，生成一个新

采样分布为

，

，利用一个建议分布候选状态X*，根据接受概率选择接受还是拒绝x*。密度函数值变大，此时接受概率α=1，接受X*，Xt=X*。若此时接受概率先生成一个服从[0,1]上的均匀分布的随机数u，如果u<α，则接受X*，Xt=X*，否则，Xt=Xt-1

，使用前面的状态。这样的过程一直持续到收敛即可，得到的样本就是采样的分布的样本。如果建议分布为一个条件概率，把它作为转移概率与采样分布的乘积处理方法是左边乘以右边，右边乘以左边，可以保证相等，即不满足细致平稳条件，即称为接受概率，则称为转移概率，构造出了满足细致平稳条件的马尔科夫链的转移概率。建议分布具有对称性时，即只需比较

大小建议分布具有对称性时，即p(x)因此取接受概率为,此时为Metropolis算法如果分布不具有对称性，可以要求即状态之间转换是可逆的，则接受概率为此时的算法称为Metropolis－Hastings算法对称分布指的是分布的密度函数具有对称性，如正态分布、柯西分布、t分布、均匀分布等，这些均可以作为建议分布。例1

从密度函数为中采样，建议分布采用U(0,1)，即使用Metropolis算法进行采样第一步。初始值随机选取U(0,1)上的随机数X*为0.3，此时接受概率随机选取U(0,1)上的随机数为0.8，因为u>α，则X1=X0=0.5。

第二步。随机选取U(0,1)上的随机数X*为0.7，此时接受概率随机选取U(0,1)上的随机数为0.3，因为u<α，则X3=X*=0.4。

依次进行下去，可抽取出n个服从密度函数为

的样本随机选取U(0,1)上的随机数为0.8，因为u<α，则X2=X*=0.7。

第三步。随机选取U(0,1)上的随机数X*为0.4，此时接受概率importrandomfromscipy.statsimportuniformimportmatplotlib.pyplotaspltN=5000x=[0foriinrange(N)]x[0]=0.5t=0#统计拒绝候选点的个数kk=0y=[]#执行遍历过程whilet<N-1:t=t+1#建议分布为均匀分布x_star=uniform.rvs(loc=0,scale=1,size=1,random_state=None)q=x_star/x[t-1]#计算接受概率alpha=min(1,q)#从均匀分布中随机抽取一个数uu=random.uniform(0,1)ifu<alpha:#接受候选点x[t]=x_starelse:#拒绝候选点x[t]=x[t-1]k+=1#计算y=2xforiinrange(N):k=2*x[i]y.append(k)#采样分布的密度函数的绘制plt.scatter(x,y,color='red')num_bins=10#采样得到的样本值的直方图plt.hist(x,num_bins,density=True,color='blue',alpha=0.7)plt.show()print('拒绝的候选点的个数为',k,'拒绝的比例为',k/N)拒绝的候选点的个数为[1.13680721]拒绝的比例为[0.00022736]例2目标平稳分布（采样分布）是一个期望为5，标准差为3的正态分布，而建议分布仍然为正态分布，标准差为3，期望等于上一个状态的分布的期望，使用Metropolis算法进行采样。首先从建议分布N(x，3)中生成候选值,接受概率为再随机选取U(0,1)上的随机数，最后判断是否接受该候选状态。importrandomfromscipy.statsimportnormimportmatplotlib.pyplotasplt#定义平稳分布正态分布的密度函数defnorm_p(x):y=norm.pdf(x,loc=5,scale=3)returny#采样的次数N=5000x=[0foriinrange(N)]sigma=3#设置初始值t=0#统计拒绝候选点的个数kk=0#执行遍历过程whilet<N-1:t=t+1#状态转移，进行随机抽样，建议分布为正态分布x_star=norm.rvs(loc=x[t-1],scale=sigma,size=1,random_state=None)alpha=min(1,(norm_p(x_star[0])/norm_p(x[t-1])))#计算接受概率#随机从均匀分布里抽取一个数uu=random.uniform(0,1)ifu<alpha:x[t]=x_star[0]else:x[t]=x[t-1]k+=1plt.scatter(x,norm.pdf(x,loc=5,scale=3),color='red')num_bins=100plt.hist(x,num_bins,density=True,color='blue',alpha=0.7)plt.show()print('拒绝的候选点的个数为',k,'拒绝的比例为',k/N)拒绝的候选点的个数为1430，拒绝的比例为0.286例3目标平稳分布（采样分布）是Rayleigh分布，密度函数为σ为4。建议分布采用自由度为xt的卡方分布，因为卡方分布不具有对称性，故使用Metropolis-Hastings算法进行采样，生成样本。首先从建议分布

中生成候选值X*计算接受概率再随机选取U(0,1)上的随机数，最后判断是否接受该候选状态。importnumpyasnpfromscipy.statsimportchi2fromscipy.statsimportrayleighimportrandomimportmatplotlib.pyplotaspltN=5000sigma=4x=[0foriinrange(N)]x[0]=1t=0#定义平稳分布为loc=0的Rayleigh分布defRayleigh_p(x,sigma):y=(x/(sigma**2))*np.exp(-x**2/(2*(sigma**2)))returnyk=0#执行遍历过程whilet<N-1:t=t+1#状态转移，进行随机抽样，建议分布为卡方分布x_star=chi2.rvs(df=x[t-1],size=1)alpha=min(1,Rayleigh_p(x_star,sigma)*chi2.pdf(x[t-1],df=x_star)/(Rayleigh_p(x[t-1],sigma)*chi2.pdf(x_star,df=x[t-1])))#这里进行接受概率的计算u=random.uniform(0,1)ifu<alpha:x[t]=x_starelse:x[t]=x[t-1]#k+=1plt.scatter(x,rayleigh.pdf(x,loc=0,scale=sigma),color='red')num_bins=50plt.hist(x,num_bins,density=True,alpha=0.7)plt.show()print('拒绝的候选点的个数为',k,'拒绝的比例为',k/N)拒绝的候选点的个数为2054，拒绝的比例为0.4108例4目标平稳分布（采样分布）是t分布，密度函数为自由度为n。建议分布采用正态分布，因为正态分布具有对称性，故使用Metropolis算法进行采样，生成样本。首先从建议分布

中生成候选值X*计算接受概率再随机选取U(0,1)上的随机数，最后判断是否接受该候选状态。importrandomfromscipy.statsimportnormfromscipy.statsimporttimportmatplotlib.pyplotasplt#定义平稳分布t分布deft_p(x):y=t.pdf(x,n)returnyN=2000x=[0foriinrange(N)]n=4sigma=0.01s=0k=0whiles<N-1:s=s+1#状态转移，进行随机抽样,建议分布为正态分布x_star=norm.rvs(loc=x[s-1],scale=sigma,size=1,random_state=None)alpha=min(1,(t_p(x_star[0])/t_p(x[s-1])))#计算接受概率u=random.uniform(0,1)ifu<alpha:x[s]=x_star[0]else:x[s]=x[s-1]k+=1plt.scatter(x,t.pdf(x,n),color='red')#采样分布的密度函数num_bins=10plt.hist(x,num_bins,density=True,color='blue',alpha=0.7)#采样数据的直方图plt.show()print('拒绝的候选点的个数为',k,'拒绝的比例为',k/N)标准差为0.01时，拒绝的候选点的个数为3，拒绝的比例为0.0015。（2）标准差为0.5时拒绝的候选点的个数为295拒绝的比例为0.1475（3）标准差为2时，拒绝的候选点的个数为1145拒绝的比例为0.5725。（4）标准差为10时拒绝的候选点的个数为1691拒绝的比例为0.8455。例5从密度函数为的混合正态分布中抽取容量为30的样本，使用样本数据估计比例p，根据同等无知的原则，尤其是p的范围为[0，1]，可以使用beta分布作为先验分布，我们这里使用U(0,1)。使用p的先验分布U(0,1)作为建议分布，后验分布作为采样分布。首先从建议分布

中生成候选值X*计算接受概率再随机选取U(0,1)上的随机数，最后判断是否接受该候选状态。先验分布为U(0,1)，先写出X和p的联合分布律：然后求X的边缘分布律：最后求出p的后验分布：使用p的先验分布U(0,1)作为建议分布，因此接受概率为importmatplotlib.pyplotaspltimportnumpyasnpm=3000r=[0.2,0.8]c=[]mu=[0,5]sigma=[1,1]foriinrange(m):k=np.random.multinomial(1,r)#先从0.2,0,8中取一个数，返回的是它的位置，然后取最大值z_i=np.argmax(k)#生成期望和方差均为指定数据的正态分布的随机数c_i=np.random.normal(mu[z_i],sigma[z_i])c.append(c_i)ig,ax=plt.subplots(figsize=(8,4))ax.hist(c,bins=100,density=True)#直方图importnumpyasnpimportmatplotlib.pyplotaspltmu=[0,5]#期望sigma=[1,1]#方差r=[0.2,0.8]#比例#正态分布的密度函数defmix_normal(x,mean,variance):return((1./np.sqrt(2*np.pi*variance))*np.exp(-(x-mean)**2/(2*variance)))a=np.arange(-5,18,0.01)#混合分布的密度函数y=r[0]*mix_normal(a,mean=mu[0],variance=sigma[0]**2)+r[1]*mix_normal(a,mean=mu[1],variance=sigma[1]**2)ig,ax=plt.subplots(figsize=(8,4))ax.plot(a,y)importrandomfromscipy.statsimportnormfromscipy.statsimportbetaimportnumpyasnpimportmatplotlib.pyplotaspltN=5000x=[0foriinrange(N)]mu=[0,5]sigma=[1,1]a=1b=1t=0k=0z=[]#比例r=[0.2,0.8]m=300#存放m个混合正态分布的随机数c=[]q=d(y1/y2)*(y3/y4)alpha=min(1,q)#随机从均匀分布里抽取一个数uu=random.uniform(0,1)ifu<alpha:x[t]=x_star[0]else:x[t]=x[t-1]k+=1print('p的估计值为',l)x[0]=0.1foriinrange(m):#先从0.2,0,8中取一个数，返回的是它的位置，然后取最大值z_i=np.argmax(np.random.multinomial(1,r))#生成期望和方差均为指定数据的正态分布的随机数c_i=np.random.normal(mu[z_i],sigma[z_i])c.append(c_i)#执行遍历过程whilet<N-1:t=t+1#状态转移，进行随机抽样，建议分布为beta分布，即先验分布x_star=beta.rvs(a,b,size=1)y1=x_star*norm.pdf(c,mu[0],sigma[0])+(1-x_star)*norm.pdf(c,mu[1],sigma[1])y2=x[t-1]*norm.pdf(c,mu[0],sigma[0])+(1-x[t-1])*norm.pdf(c,mu[1],sigma[1])y3=(x[t-1]**(a-1))*((1-x[t-1])**(b-1))y4=(x_star**(a-1))*((1-x_star)**(b-1))

ifu<alpha:x[t]=x_star[0]else:x[t]=x[t-1]k+=1#forjinrange(0,N):#plt.plot(x)num_bins=10#采样数据的直方图plt.hist(x,num_bins,density=True,color='blue',alpha=0.7)plt.show()print('拒绝的候选点的个数为',k,'拒绝的比例为',k/N)l=np.mean(x)print('p的估计值为',l)拒绝的候选点的个数为4627拒绝的比例为0.9254p的估计值为0.19821454111333214拒绝的候选点的个数为4627拒绝的比例为0.9254p的估计值为0.19821454111333214（2）已知采样分布为混合正态分布，建议分布为正态分布，因为正态分布具有对称性，故使用Metropolis算法进行采样，生成样本。首先从建议分布N(3，10)中生成候选值,计算接受概率再随机选取U(0,1)上的随机数，最后判断是否接受该候选状态importrandomfromscipy.statsimportnormimportmatplotlib.pyplotaspltN=5000x=[0foriinrange(N)]mu=[0,5,3]sigma=[1,1,10]#比例pp=0.2#设置初始值t=0k=0x[0]=0defmix_norm(x):y=p*norm.pdf(x,loc=mu[0],scale=sigma[0])+(1-p)*norm.pdf(x,loc=mu[1],scale=sigma[1])#混合正态分布的密度函数

returny#执行遍历过程whilet<N-1:t=t+1#状态转移，进行随机抽样，建议分布为正态分布x_star=norm.rvs(loc=mu[2],scale=sigma[2],size=1,random_state=None)#计算接受概率alpha=min(1,(mix_norm(x_star[0])/mix_norm(x[t-1])))#随机从均匀分布里抽取一个数uu=random.uniform(0,1)ifu<alpha:x[t]=x_star[0]else:x[t]=x[t-1]k+=1#forjinrange(0,N):#plt.plot(x)num_bins=10plt.hist(x,num_bins,density=True,alpha=0.7)plt.show()print('拒绝的候选点的个数为',k,'拒绝的比例为',k/N)l=np.mean(x)print('p的估计值为',l)拒绝的候选点的个数为3999拒绝的比例为0.7998模拟退火算法的思想最早是由Metropolis等提出的。其出发点是基于物理中固体物质的退火过程与一般的组合优化问题之间的相似性。物理退火过程由以下三部分组成：1、加温过程。其目的是增强粒子的热运动，使其偏离平衡位置。当温度足够高时，固体熔为液体，从而消除系统原先可能存在的非均匀状态，使随后进行的冷却过程以某一平衡态为起点。溶解过程与系统的熵增过程相联系，系统能量随温度升高而增大。2、等温过程。对于与周围环境交换热量而温度不变的封闭系统，系统状态的自发变化总是朝着自由能减少的方向进行的，当自由能达到最小时，系统达到平衡状态。3、冷却过程。使粒子热运动减弱并渐趋有序，系统能量下降，得到低能的晶体结构。模拟退火算法（SimulatedAnnealing，SA）其中，加温过程对应算法的设定初温，等温过程对应算法的Metropolis抽样过程，冷却过程对应控制参数的下降。这里能量的变化就是目标函数，要得到的最优解就是能量最低态。Metropolis准则是SA算法收敛于全局最优解的关键所在，Metropolis准则以一定的概率接受恶化解，这样就使算法跳离局部最优的陷阱。在模拟退火算法中，一般采用Metropolis准则，具体如下当T趋向于0时，由概率分布可知温度的最小点。在模拟退火方法中，最主要的部分是要缓慢降低算法的虚拟温度。这就引入一个方案：冷却进度表。它是模拟退火算法能否得到全局最优解的关键方法。冷却进度表是指从某一高温状态T向低温状态冷却时的降温函数。有两种衰减函数的选择，一是几何冷却方式的衰减函数：α用来控制降温的速率，为温度衰减系数，α越大，衰减速率越小。α是接近于1的常数，一般为0.5-0.99。该衰减函数对控制参数的衰减量是随算法进程递减的。二是逆对数冷却方式。假设初始温度为T0，当前的温度为T(t),设计冷却进度表要遵循三个原则：（1）有足够大的初始温度(2）选择合适的衰减函数。（3）合适的马尔可夫链长度k。马尔可夫链的长度选取原则是应使温度ti上产生的序列达到准平衡态。（4）停止温度应该足够小。例6

求一元函数在区间[0,1]上的最小值和最大值。可以大体看出，在x=0.8991时取得最小值，在x=0.21时取得最小值。接下来使用模拟退火算法进行计算，（1）设置初始温度t0=1000，温度的最小值tmin为0.000001，达到最小值就不需要再降了；（2）用x表示自变量，在0到1之间，low=0，high=1，y表示函数值；（3）当温度t0大于tmin时，每个温度迭代k=50次，执行如下步骤：1）计算函数值y。2）定义x的扰动函数，任取一个[0,1]上的随机数u，如果大于0.5，则xnew=

x+(high-x)*u，否则xnew=

x-(x-low)*u。如果xnew在low到high之间，计算函数值ynew；否则舍去。3）如果ynew<y，则x=xnew，否则，采用Metropolis准则：4）对每个温度，判断条件ynew>ymax，如果成立，则ymax=ynew，

xbestM=xNew。求出最大值ymax及对应的x。判断条件ynew<ymin，如果成立，则ymin=ynew，xbestm=xNew。求出最大值ymin及对应的x。5）选择降温函数α为温度衰减系数importnumpyasnpimportmathimportmatplotlib.pyplotaspltdefFunction(x):y=x+10*math.sin(5*x)+7*math.cos(4*x)returny#定义自变量x的扰动defDisturb(low,high,x):ifnp.random.random()>0.5:xnew=x+(high-x)*np.random.random()else:xnew=x-(x-low)*np.random.random()returnxnewt0=1000#初始温度tmin=0.000001#温度的最小值#初始值，求最值的区间low=0high=1x=np.random.uniform(low,high)alpha=0.98#迭代次数k=50y=0t=0#最大值和最小值的初值ymax=0ymin=0#模拟退火算法whilet0>=tmin:#每一个温度循环k轮foriinrange(k):y=Function(x)#利用变换函数生成一个xnewxNew=Disturb(low,high,x)if(low<=xNewandxNew<=high):yNew=Function(xNew)ifyNew<y:x=xNewelse:#metropolis算法p=np.exp(-(yNew-y)/t0)r=np.random.uniform(low=0,high=1)ifr<p:x=xNew#此处>求最大值,<求最小值ifyNew>ymax:ymax=yNewxbestM=xNewifyNew<ymin:ymin=yNewxbestm=xNewt+=1t0=alpha*t0#降温函数print('自变量x=',xbestM,'函数值y=',ymax,'迭代次数为',t)print('自变量x=',xbestm,'函数值y=',ymin,'迭代次数为',t)变量x=0.22573930378085352函数值y=13.599327598073529迭代次数为228自变量x=0.895162543645095函数值y=-15.162367482544628

若在[0,2]范围内求最小和最大值，则把N，high都改为2即可。自变量x=1.5734304429268087函数值y=18.572174578287733迭代次数为1026自变量x=0.8936258923609317函数值y=-15.16377975376529迭代次数为1026若在[0,6]范围内求最小和最大值，则把N，high都改为6即可。自变量x=1.581805