ch蒙特卡罗方法导论概要实用

会计学1ch蒙特卡罗方法导论概要实用2引言本章的目的是简单地介绍参数估计的蒙特卡罗方法基础。本章将涵盖蒙特卡罗估计方法的几个重要方面。本章的目标是给出蒙特卡罗方法的定义，并在一个简单易懂的背景中考察一些基本方法。对置信区间和收敛性等重要问题也作了简要的讨论。在本章的所有论述中，都假设估计器所采用的观测值是相互独立的，而在下一章我们将放宽这一假设。第1页/共59页3蒙特卡罗仿真建立在几率游戏的基础上，它之所以命名为“蒙特卡罗”，因为它以一个以赌博著称的地中海城市得以命名。蒙特卡罗仿真是指那些利用蒙特卡罗方法估计系统参数如误比特率（BER）的仿真，而蒙特卡罗估计则是指通过内在的随机试验来估计参数值的过程（一）相对频率蒙特卡罗估计是基于概率的相对频率。为定义相对频率，第一步是确定随机试验和一个感兴趣的事件A。由概率论可知，随机试验中试验结果无法准确预测，而只能用统计方法加以描述。一、基本概念第2页/共59页4一个最基本的随机试验就是掷硬币，感兴趣的试验结果有两个，可用集合｛正面，反面｝描述。在掷硬币之前，并不知道会出现哪一个结果，但是，如果已知硬币是“公平”的，或者说是无偏的，则集合｛正面，反面｝中每一种结果出现概率会相等，且相互独立。随机试验的性能决定了试验结果。

随机试验中的事件可以是一个结果或几个结果的集合。以数字通信系统为例，随机试验可定义为发送一个二进制数1,接收机输出端的结果为对所发送二进制符号的估计，它是二进制数0或者二进制数1，我们所感兴趣的事件可能是在发送1的过程中所产生的差错。要确定系统的BER，得估计在发送1的条件下接收到0这一事件的概率。第3页/共59页5在定义了随机试验和感兴趣的事件A以后，我们来考虑蒙特卡罗方法的下一个步骤，即进行大量的随机试验，试验次数为N，以NA表示事件A发生的次数。将事件A发生的概率近似为相对频率，其定义为NA/N。在相对频率的意义下，事件A的概率可以通过重复无限多次随机试验来求得，即在估计数字传输系统里的差错概率这个例子中，N是总的发送比特数或符号数（系统实际发送的或仿真的），NA是差错发生的次数（测量出来的或仿真得到的）。第4页/共59页6在蒙特卡罗仿真中，显然必定有N＜∞，数值NA/N就是Pr(A)的估计器，该估计器记作Pr(A)。必须注意到，由于试验的随机性，在试验次数N有限的情况下，NA是随机变量，因此Pr(A)也是随机变量。该随机变量的统计性质决定了估计器的精度，因而也决定了仿真的质量。（二）无偏和一致估计器

蒙特卡罗估计器必须满足几个重要的性质,其在实际中才有用。首先，我们希望蒙特卡罗估计器是无偏的。也就是说，如果A是参数A的估计值，我们希望第5页/共59页7换句话说，在平均意义下可以得到正确的结果。假设进行了多次蒙特卡罗仿真，得到感兴趣的随机变量的一组估计值，显然我们希望这些估计值有比较小的方差。如果估计值是无偏的并具有小的方差，则估计器所产生的估计值会聚集在待估计参数真值的周围，且有较小的散布范围。除非所研究的事件是统计独立的，用解析的方法确定蒙特卡罗估计器的方差通常是困难的，但是几乎可以肯定，估计值的方差会随着仿真运行时间（相应随机试验重复的次数）的增加而减少。我们把满足这一性质的估计器称为一致的。对于一致估计器，当N→∞时，A2→0，其中N代表随机试验重复的次数。而对于无偏和一致估计器，误差第6页/共59页8具有零均值，而其方差e2在N→∞收敛到零。但是，这一收敛过程通常非常缓慢。（三）蒙特卡罗估计考虑确定一个形状不规则的区域的面积作为一个蒙特卡罗估计器的简单例子。假设待估计面积的区域完全包含在一个面积已知的方框中，随机试验定义为在包围方框中随机抽取采样，事件A定义为采样点落在面积待确定的区域内。要得到一个未知面积的无偏估计器，随机采样点必须均匀分布在面积已知的包围区域中，利用具有两个均匀随机数发生器的计算机程序就可以很容易地实现这一点。产生N＝500个均匀分布采样点的结果如图9-1所示，这500个点是利用如下MATLAB代码产生的：第7页/共59页9x=rand（1，500)；y=rand（1，500)：plot(x，Y，‘k+’）axissquare下一步是定义感兴趣的事件A。我们要估计如图9-2所示的爆炸形（sunburst)区域的面积，分别定义Nbox和Nsunburst为落在包围方框中和爆炸形区域的采样点数。因为采样点在包围方框中是均匀分布的，Asunburst/Abox近似为采样点落在爆炸形区域内的数目与落在包围方框中的数目之比，即Nsunburst/Nbox。也就是从而有第8页/共59页10在采样点为均匀分布这一条件下，随着采样点数的增加，近似精度会不断提高。为了以简单直接的方式说明蒙特卡罗方法，我们考虑值的蒙特卡罗估计器。注意，估计器是随机仿真，因为它是随机试验的仿真，因而这个例子可作为本书以后章节的导论。第9页/共59页11图9-1均匀分布的随机样点第10页/共59页12第11页/共59页13（四）的估计估计数值的方法之一是用一个具有单位面积的正方形包围一个馅饼状（pie-shaped）的区域，即单位圆的第一象限。如图9-3所示为这种情况，以及总的采样点数Nbox。如果正方形在x轴上所占的区间是(0,1)，在y轴上所占的区间也是(0,1)，显然Abox＝1，馅饼状区域（四分之一圆）的面积为由此可得第12页/共59页14第13页/共59页15假设采样点是均匀分布的，则Npie_slice和Nbox之比构成了Apie_slice/Abox的无偏估计，所以的估计器，记作,可写成所以对值的估计，只要在包围方框中撒上均匀分布的采样点，再对落在内切圆中的采样点计数，再利用式（9-9）即可完成。第14页/共59页16

例9-1上述过程可以方便地用MATLAB程序实现。如图9-4所示的结果中，有5个的估计，每一个都是基于500次重复随机试验，所得的的五个估计值用以下向量表示如果对5个结果进行平均，则＝3.0736，这样的结果等价于2500次的试验结果。用来产生该结果的MATLAB程序如下：c9_estimatepi.m■

第15页/共59页17图9-4的蒙特卡罗估计，5轮实验，每轮500次重复第16页/共59页18图9-4的蒙特卡罗估计，10轮实验，每轮10000次重复第17页/共59页19尽管这个例子非常简单，它却具有所有蒙特卡罗仿真的许多重要特点，其中都有一个测试条件以及两个计数器。随机试验每进行一次，第一个计数器就增加1，如果测试条件每满足一次，第二个计数器就增加1。在数字通信系统的仿真中，仿真目的是估计误比特率，测试条件就是在发送给定的比特或数据符号时是否发生了差错。在仿真中，每当发送了一个比特或数据符号，第一计数器就增加1。每当观测到一次差错，第二个计数器就增加1。第18页/共59页20二、在通信系统中的应用--AWGN信道为了用蒙特卡罗方法估计通信系统的性能，让N个符号通过系统（实际上是系统的计算机仿真模型），并计算发送差错的个数Ne。如果在N次的符号发送中有Ne次差错，则符号差错概率为问题是这一估计是有偏还是无偏呢？是否为一致估计？为了在一个尽可能简单的背景下研究这些重要的问题，我们假设信道是AWGN（加性高斯白噪声）信道。在AWGN条件下，由信道噪声所产生的差错相互独立，而且在N次符号发送中出现的差错次数Ne可以用二项分布来描述，所以我们转而先对二项分布作稍微详细一点的了解。在讨论完二项分布之后，我们以式（9-11)作为两个高度理想化的通信系统的符号差错概率估计器。第19页/共59页21（一）二项式分布我们下面的任务是确定PE的统计特性。首先是确定Ne的均值和方差。对于相互独立的差错事件，N次符号发送中有Ne次差错的概率服从以下的二项式分布其中是二项式分布系数，PE是单次发送时的差错概率。第20页/共59页22服从二项式分布的随机变量的均值和方差很容易求得。Ne的均值为Ne的方差为把这些结果代入式（9-11），则差错概率的蒙特卡罗估计器的均值为第21页/共59页23将式（9-14）代入式（9-16）可得可见差错概率的蒙特卡罗估计器是无偏的。而差错概率蒙特卡罗估计器的方差为将式（9-15）代入式（9-18）有第22页/共59页24可知估计器是一致的，因为当N->∞

时，方差会不断下降。注意式（9-17）和式（9-19）都假设了二项式分布，这只有在差错事件相互独立这一条件下才成立。在利用蒙特卡罗仿真来估计通信系统某一性能参数时，例如符号差错概率，我们都希望估计是无偏和一致的，如果估计器是无偏的，则蒙特卡罗仿真可以在平均意义下提供正确的结果。除此以外，为使估计器实用，它还必须具有小方差，这样估计值就会以较高的概率落在参数真值的附近。

如果估计器具有无偏性和一致性，我们可以仿真更多的发送符号，以便记录到更多的差错，进而减小估计器的方差。式（9-19）表明，为了达到一个给定的方差，必须记录到一定数量的差错，这反过来又要求一定的仿真执行时间。但是对于给定方差，式（9-19)存在的一个现实问题是无法事先确定所需要的试验次数N，因为在进行仿真之前无法得知PE。不过，在实际过程中，我们也许可以通过运用界或者其他的分析方法，第23页/共59页25将PE的值确定在一个数量级之内，所以式（9-19)可能还是有用的。对于有偏但是一致的估计器，估计值收敛到一个不正确的值，这显然是一种非常不理想的情况，除非我们知道如何消除这个偏差。尽管知道估计器的性质非常重要，但在许多情况下很难证明估计器是无偏的和一致的。必须强调，尽管本节所得到的所有结果很理想，但它们成立的条件是信道噪声引起的差错相互独立，从而对应的差错概率是二项分布的。如果差错事件是相关的，例如在带限信道的情况下，这里所给出的结果不再正确，从而面临的问题就更加困难。不过，即使差错事件不是相互独立的，式（9-11)仍然是一个有效的差错概率估计器。

第24页/共59页26

例9-2在差错事件相互独立的情况下，可用掷硬币事件来模拟二进制发送，N个符号的发送可建模为对一枚不均匀的硬币进行N次投掷。可以假设第i次投掷结果为“反面”相当于对第i次发送作出了正确的判决，而第i次投掷结果为“正面”相当于对第i次发送作出了错误的判决。在这个例子中，掷币过程所具有的统计特性由仿真来决定。因为投掷是相互独立的，所以这个试验模拟了AWGN信道中的二进制数据传输。假设投掷结果为“反面”（没有差错）的概率为1-p，结果为“正面”（有差错）的概率为p。我们希望通过投掷N次硬币来估计p值。p的蒙特卡罗估计器为第25页/共59页27其中NHeads表示在连续投掷N次硬币中出现“正面”的次数。当然，在N次投掷中，NHeads可以是0到N中的任何数，但是出现k次“正面”的概率为因此，为了估计NHeads的统计分布并确定p的估计器p，我们必须重复进行这个试验多次，比如M次。以下是仿真投掷硬币试验的MATLAB程序：c9_cointoss.m执行该程序的结果如图9-5所示。上部分显示的直方图，每次试验的结果以及理论结果如图下部所示。注意，正如拉普拉斯近似所预测的那样，理论结果近似高斯分布。二项式系数是利用函数nkchoose计算出来的。第26页/共59页28图9-5投掷硬币试验的结果(M=2000,N=500)第27页/共59页29图9-5投掷硬币试验的结果(M=20000,N=500)第28页/共59页30图9-5投掷硬币试验的结果(M=200000,N=500)第29页/共59页31用这样的方法来计算二项式的系数，是为了阐明在n值很大时比较有用的一个算法。尽管MATLAB具有很大的动态范围，这个例子也不需要前面程序所示的方法，但在采用其他程序语言时，此方法通常是有用的。■（二）两个简单的蒙特卡罗仿真在这一节里，我们首次考虑通信系统的蒙特卡罗仿真。先作以下假设：.在发射机中没有进行脉冲成形。.假设信道是AWGN信道。.信源输出端的数据符号是相互独立和等概率的。.在系统中没有滤波处理，因此不存在码间干扰。有了这些假设，本节考虑的系统和相应的仿真都极其简单。系统是易于解析处理的，任何初学数字通信原理的读者都可以根据目测写出这个差错概率。第30页/共59页32尽管下面的例子很简单，但是还是很重要的，因为从中得出的几个重要观察对后面的研究都很有用。除此以外，还确立了仿真程序的基本结构。我们还将看到，当蒙特卡罗仿真应用于更贴近我们要集中研究的主题（数字通信系统）方面的问题时，其性能究竟如何。如图9-6所示是基本的仿真模型。由于没有滤波，系统的时延为零，因此，在比较发送符号d[n]和接收符d[n]之前，没有必要用时延模块来将相应符号进行对齐或同步。但是在图9-6中我们仍然用点画线表示出时延模块，以提醒读者几乎所有的仿真都需要这个重要部分。由于前面的假设，唯一的差错源是信道噪声。因此我们采用定义同相和正交信号分量xd(t)和xq(t)的方法，以便使它们给出的是信号的信号空间分量，而不是时域波形的采样。这么做的好处是信号空间分量只要用每个发送符号的一个采样就可以给出了。基于每个符号一个采样的仿真执行起来非常快。第31页/共59页33利用这一方法，调制器输出端假设的带通信号可表示为其中Ac表示载波幅值，km是跟调制有关的常数，d[n]是第n个数据符号（d[n]＝0或1），是参考相位。直观上可知，x(t,n)的复包络只是符号序数n的函数，表示为在本例中，我们假设＝0，因此，在图9-6中有第32页/共59页34第33页/共59页35对如图9-6所示的系统，为了确定和画出以Eb/N0为函数的BER，保持Eb恒定，并让噪声功率在感兴趣的范围内递增，这就要求对图9-6中的噪声发生器输出端的噪声功率进行校准。噪声方差和噪声功率谱密度（PSD）的关系为信噪比SNR定义为Eb/N0，其中fs为采样频率。因此有如果能量Eb和采样频率都归一化为1，则有第34页/共59页36该表达式用于确定后面仿真中的噪声标准偏差。

例9-3（二进制相移键控，BPSK）为了产生二进制PSK信号的同相和正交信号空间分量，在式（9-24）和式（9-25)中令Ac=1和km＝，由此得到第35页/共59页37于是，在d[n]=0和d[n]=1这两种情况下都有xq(t)=0,这就给出了二进制PSK的信号空间表示，如图9-7所示，图中1是信号空间的基函数。由于信号空间是一维的（有一个基函数产生），在仿真中只需要产生信号和噪声的直接分量，而图9-6中所示的正交通道可以丢弃掉。图9-7也同时示出了判决区域0和1，如果接收信号点落在0区域（在1=0的右边区域），接收机做出判决d[n]=0;如果接收信号点落在1区域（在1=1的左边区域），接收机做出判决d[n]=1。对于AWGN条件下，等概率、等能量的信号，接收机的阈值总是为零。所以判决规则为以上这些讨论给出如下MATLAB仿真程序：c9-MCBPSK.m第36页/共59页38图9-8二进制相移键控第37页/共59页39对于每一个SNR值，发送Nsymbols＝10000个符号，运行结果如图9-8所示。注意，当SNR增加时，BER估计器的可靠度会变差，这是由于差错发生次数减少的缘故。这一现象表明，仿真所发送的符号数可能与SNR有关，或者是连续运行仿真程序，直到对每一个SNR值都记录到相同数目的差错为止。■

例9-4（二进制频移键控，BFSK）为了产生二进制FSK信号空间的同相和正交分量，令式（9-24）和式（9-25）中的km=/2，从而得到第38页/共59页40类似地有，这样所得到的二进制FSK信号空间表示如图9-9所示，其中1和2是信号空间的基函数。（根据第4章的内容，对于二维空间，可以认为基函数定义了低通复包络信号的同相和正交分量。）因为二进制FSK信号空间是二维的，仿真必须产生信号和噪声的同相和正交分量。判决区域也在图9-9中表示出来了，如果接收信号点落在0区域内（判决边界的右下部分），接收机判决d[n]=0;如果接收号点落在1区域内（判决边界的左上部分），接收机判d[n]=1。注意，对于一个表示接收信号（图9-6中的yd[n]和yq[n]）给定信号空间点，判决规则为第39页/共59页41仿真由以下的MATLAB程序实现：c9-MCBFSK.m第40页/共59页42对于每一个SNR值，设定Nsymbols＝10000个符号，程序执行的结果如图9-10所示。同样要注意到，随着SNR的增加，由于发生较少的差错数，估计器的可靠性会下降。在前一个例子中已经介绍了适当的校正方法。■图9-10二进制频移键控第41页/共59页43三、蒙特卡罗积分蒙特卡罗积分这一主题很自然地出现在我们的通信仿真研究中。我们已经知道对于AWGN信道，通过采样积分--清除（integrate-and-dump）检测器的输出所得到的充分统计量V，是一个高斯随机变量，其均值由数据符号决定，方差由信道噪声决定。在d[n]=0和d[n]=1的条件下，条件概率密度函数如图9-11所示，其中kT是接收机阈值。在d[n]=1条件下的条件差错概率为对于Pr(E|d[n]=0)，也可推出类似的表达式。从而，系统的差错概率估计第42页/共59页44对蒙特卡罗积分的简要研究有助于我们进一步加深对蒙特卡罗仿真方法的理解。例如，蒙特卡罗积分的研究提供了一个简单的环境，可用来说明蒙特卡罗估计器的收敛性。牵涉到积分值的估计问题。第43页/共59页45（一）基本概念假设我们想求以下积分其中g(x)是一个在积分区域内有界的函数。由基本的概率论，可知函数g(x)的总体均值为其中fx(x)是随机变量X的概率密度函数。如果随机变量X的概率密度函数在区间（0,1)上满fx(x)=1，而在其他地方为零，则E{g(x)}=I。因此，如果U是在区间（0，1)上均匀分布的随机变量，那么就有第44页/共59页46利用相对频率的观点可以得到因此，对被积函数进行仿真，再在（0,1）区间上对它进行N次采样，采样点的平均值可用来估计积分值。系统的蒙特卡罗仿真基本上也是这个思路，因为通常无法获得充分统计量在差错区域内的解析表达式，所以采用系统仿真来产生统计量采样。如果无法求式（9-41)中的极限（实际应用中往往如此），则得到积分的一个近似。将此近似记作I,得蒙特卡罗估计器为第45页/共59页47总之，通过对函数g(x)在N个均匀分布的采样点上求值再取平均，就可实现积分的估计器。这种方法适用于任何常义积分。通过简单的变量代换，可用蒙特卡罗方法来估算任意区间上的常义积分。例如，积分可以通过变量代换y=(x-a)/(b-a)变成标准形式并求得

例9-5为了用蒙特卡罗方法估计值，只需找到一个定积分值为的函数。很容易就可想到如下积分第46页/共59页48因此，我们就可以利用式（9-42）定义的算法，求解积分值I，并将所得结果乘以4。显然我们所能做到的顶多只是选用一个很大但有限的N，据此获得的不是的准确值，而是其近似值。因此，这个积分值，也就是的估计值，是一个随机变量。如图9-12所示为五次估计值的结果，每一次都基于500个试验。这五个估计值为对五次结果作平均得利用蒙特卡罗积分估计值的MATLAB程序如下：c9_example5．m第47页/共59页49图9-12用蒙特卡罗积分估计5轮Moncar实验，每轮500次实验第48页/共59页50图9-12用MC积分估计10轮Moncar实验，每轮10000次实验第49页/共59页51（二）收敛性假设要估计积分值I,同时取得了N个随机观测或采样值Xi。I的估计器为假设N个观测Xi为独立同分布。样本的算术平均为因此估计值是无偏的。由于已经假设观测的独立性，样本方差为第50页/共59页52该式表明积分估计器是一致的。假设Xi=g(Ui)，样本的方差，记作x2，表达为因此，给定被积函数g(u)，可以确定要达到某一给定误差方差所需要的N值。因为估计器是一致的，只要N充分的大，就可以得到精确的积分估计。同时可知，如果样本g(ui)具有小方差，也可以得到对I的精确