通信系统的计算机模拟第十二讲幻灯片.ppt

《通信系统的计算机模拟第十二讲幻灯片.ppt》由会员分享，可在线阅读，更多相关《通信系统的计算机模拟第十二讲幻灯片.ppt（59页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、通信系通信系统的的计算机模算机模拟第十二第十二讲1第1页，共59页，编辑于2022年，星期三关于期末报告关于期末报告l下下次课为下下次课为 Project 报告时间报告时间l今天请提交纸质报告一份（首页注明姓名、学号信息）今天请提交纸质报告一份（首页注明姓名、学号信息）l电子版和电子版和Matlab程序程序 email到到:l请将请将email标题标题 设为设为 Simulation2008姓名学号姓名学号l请准备好报告请准备好报告PPT,每人报告时间为（每人报告时间为（1022）分钟分钟l如有特殊时间需求请提前声明如有特殊时间需求请提前声明2第2页，共59页，编辑于2022年，星期三第十一讲

2、第十一讲 回顾回顾l将仿真产生的数据处理成将仿真产生的数据处理成有用有用有用有用的形式的形式l波形、眼图和散点图，书上的例子要学会，典型的标准模块波形、眼图和散点图，书上的例子要学会，典型的标准模块l估计估计l直方图直方图l功率密度估计功率密度估计l周期图周期图l带数据窗的周期图带数据窗的周期图l分段周期图分段周期图l增益，延迟和信噪比增益，延迟和信噪比3第3页，共59页，编辑于2022年，星期三蒙特卡罗方法导论蒙特卡罗方法导论l一个简单易懂的场院景中考察一些基本方法一个简单易懂的场院景中考察一些基本方法l对置信区间和收敛性等重要问题也作了简单的对置信区间和收敛性等重要问题也作了简单的讨论讨论

3、l假定设估计器所采有的观测值是相到独立的假定设估计器所采有的观测值是相到独立的4第4页，共59页，编辑于2022年，星期三9.1 蒙特卡罗方法的基本概念蒙特卡罗方法的基本概念 l蒙特卡罗仿真建立在几率游戏的基础上蒙特卡罗仿真建立在几率游戏的基础上l“蒙特卡罗蒙特卡罗”赌博著称的地中海城市。赌博著称的地中海城市。l蒙特卡罗仿真是那些利用蒙特卡罗方法估计系蒙特卡罗仿真是那些利用蒙特卡罗方法估计系统参数如误比特率（统参数如误比特率（BER）的仿真）的仿真l而蒙特卡罗估计同是指通过内在的随机试验来而蒙特卡罗估计同是指通过内在的随机试验来估计参数值的过程。估计参数值的过程。5第5页，共59页，编辑于20

4、22年，星期三9.1.1相对频率相对频率 l蒙特卡罗估计是基于概率的相对频率解释的蒙特卡罗估计是基于概率的相对频率解释的l为定义相对频率，第一步是确定随机试验和一个感兴趣的事件为定义相对频率，第一步是确定随机试验和一个感兴趣的事件A。l随机试验中试验结结果无法准确预测，用统计的方法加以描述。随机试验中试验结结果无法准确预测，用统计的方法加以描述。l一个最基本的随机试验就是掷硬币，正面、反面一个最基本的随机试验就是掷硬币，正面、反面l 随着机试验中的事件可以是一个结果或几个结果的集合随着机试验中的事件可以是一个结果或几个结果的集合l以数学字通信系统为例，随机试验可定度为发送一个二进制数以数学字通

5、信系统为例，随机试验可定度为发送一个二进制数1，接收机输出端结果为对所，接收机输出端结果为对所发送二进制符号的估计，它是二进制数发送二进制符号的估计，它是二进制数0或者二进制数或者二进制数1l我们所感兴趣的事件可能在发送我们所感兴趣的事件可能在发送1的过程中所产生的差错的过程中所产生的差错l要确定系统的要确定系统的BER，得估计在发送，得估计在发送1的条件下接收到的条件下接收到0这一事件的概率。大量的随机这一事件的概率。大量的随机试验，试次数为试验，试次数为N，以表示事件，以表示事件A发生的次数。将事件发生的次数。将事件A发生的概率近似为相对频率，发生的概率近似为相对频率，其定义为。其定义为。

6、l在相对频率的意义下，事件在相对频率的意义下，事件A的概率可以通过重复无限多次随机实验获得的概率可以通过重复无限多次随机实验获得N是总的发送比特数或符号数，是总的发送比特数或符号数，NA错误数错误数6第6页，共59页，编辑于2022年，星期三讨论讨论 由于试验的随机性，在试验次数由于试验的随机性，在试验次数N有有限的情况下，限的情况下，NA是随机变量，是随机变量，Pr(A)也是随机也是随机变量。变量。l估计器估计器 l（仿真）随机变量的统计性质决定了估计器精（仿真）随机变量的统计性质决定了估计器精度，因而也决定了仿真的质量。度，因而也决定了仿真的质量。7第7页，共59页，编辑于2022年，星期

7、三9.1.2无偏和一致估计器无偏和一致估计器 l蒙特卡罗估计器必须满足几个重要性质在实际中才有意义蒙特卡罗估计器必须满足几个重要性质在实际中才有意义l首先，我们希望蒙特卡罗估计器是无偏的。如果首先，我们希望蒙特卡罗估计器是无偏的。如果A的估计值的估计值 ，我们希望，我们希望换句话说，在平均意义下可以得到正确的结果。换句话说，在平均意义下可以得到正确的结果。l估值有比较小的方差估值有比较小的方差l如果估计值是无偏的并具有小的方差，则估计器所产生的估计值会聚集在待估计如果估计值是无偏的并具有小的方差，则估计器所产生的估计值会聚集在待估计参数真值的周围，且有较小的散布范围。参数真值的周围，且有较小的

8、散布范围。l除非所研究的事件是统计独立的，用解析的方法确定蒙特卡罗估计器的方差通常是困除非所研究的事件是统计独立的，用解析的方法确定蒙特卡罗估计器的方差通常是困难的，但是几乎可以肯定，随着估计值数量的增加而减少，称为一致的。难的，但是几乎可以肯定，随着估计值数量的增加而减少，称为一致的。l对于一致估计器，当对于一致估计器，当 时，时，l而对于无偏和和一致估计器，误差而对于无偏和和一致估计器，误差 具有零均值，而其方差具有零均值，而其方差 在在 收敛为收敛为0。收敛过程通常非常缓慢。收敛过程通常非常缓慢 8第8页，共59页，编辑于2022年，星期三9.1.3蒙特卡罗估计蒙特卡罗估计 l考虑确定一

9、个形状不规则的区域的面积作为一个蒙特卡罗估考虑确定一个形状不规则的区域的面积作为一个蒙特卡罗估计器的简单例子。计器的简单例子。l侍估计面积的区域完全包含在一个面积已知的方框中，随侍估计面积的区域完全包含在一个面积已知的方框中，随机试验定义在包围方框中随机抽取采样，事件机试验定义在包围方框中随机抽取采样，事件A定义为采样定义为采样点落在面积待确定区域内。点落在面积待确定区域内。l要得到一个未知面积的无偏估计器，随机采用样点必须均要得到一个未知面积的无偏估计器，随机采用样点必须均匀分布在面积已知范围区域中，利用具有两个均匀随机数匀分布在面积已知范围区域中，利用具有两个均匀随机数发生器的计算机程序就

10、可以很容易地实现这一点。发生器的计算机程序就可以很容易地实现这一点。9第9页，共59页，编辑于2022年，星期三Step 110第10页，共59页，编辑于2022年，星期三l下一步是定义感兴趣的事件下一步是定义感兴趣的事件A。我们要估计如图。我们要估计如图9-2所示的爆炸形（所示的爆炸形（sunburst）区域的面）区域的面积，分别定义积，分别定义 和和 落在包围中和爆炸形区域的采样点数落在包围中和爆炸形区域的采样点数 l因为采样点在包围方框中是均匀分布的，因为采样点在包围方框中是均匀分布的，近似为采样点落在钻石形区域内的数目与落近似为采样点落在钻石形区域内的数目与落在包围方框中的数目之比在包

11、围方框中的数目之比 ，即也就是，即也就是在在采样点为均匀分布这一条件下，随着采样采样点为均匀分布这一条件下，随着采样点数的增加，近似精度会不断提高。点数的增加，近似精度会不断提高。11第11页，共59页，编辑于2022年，星期三 9.1.4 的估计的估计 l估计数值估计数值的方法之一是用一个具有单位面积的正方形包围一个馅饼状（的方法之一是用一个具有单位面积的正方形包围一个馅饼状（pie-shaped）的区域，即单位圆的第一象限。如图的区域，即单位圆的第一象限。如图9-3所示为这种情况，以及总的采样点数所示为这种情况，以及总的采样点数Nbox。如果。如果正方形在正方形在x轴上所点的区间是（轴上所

12、点的区间是（0，1），在），在y轴上所占的区间也是（轴上所占的区间也是（0，1），显然馅饼状区），显然馅饼状区域（四分之一圆）的面积为域（四分之一圆）的面积为山巅一寺一壶酒，尔乐苦煞吾，把山巅一寺一壶酒，尔乐苦煞吾，把酒吃，酒杀尔，杀不死，乐尔乐。酒吃，酒杀尔，杀不死，乐尔乐。31415926535897932384626 3141592653589793238 12第12页，共59页，编辑于2022年，星期三l假设采样点是均匀分布的，则假设采样点是均匀分布的，则 和和 Nbox 的比构成的比构成 的无偏估计，所以的无偏估计，所以13第13页，共59页，编辑于2022年，星期三例例9-1lm=

13、input(Enter M,the number of experiments );ln=input(Enter N,the number of trials per experiment );lz=zeros(1,m);ldata=zeros(n,m);lfor j=1:ml x=rand(1,n);l y=rand(1,n);l k=0;l for i=1:nl if x(i)2+y(i)2=1%Fall inside pie slice?l k=k+1;l endl data(i,j)=4*(k/i);%jth estimate of pil endl z(j)=data(n,j);%S

14、tore datalendlplot(data,k)%Plot curveslxlabel(Number of Trials)lylabel(Estimate of pi)14第14页，共59页，编辑于2022年，星期三Result如果对5个结果进行平均，则如果对如果对5个结果进行平均，则个结果进行平均，则 ，这样的结果等价于，这样的结果等价于2500次的试验结果次的试验结果 15第15页，共59页，编辑于2022年，星期三讨论讨论l 这个例子具有所有蒙特卡罗仿真的许多重要特点这个例子具有所有蒙特卡罗仿真的许多重要特点l其中都有一个测试条件以及两个计数器其中都有一个测试条件以及两个计数器l随机

15、试验每进行一次，第一个计数器就增加随机试验每进行一次，第一个计数器就增加1，如果测试条件每满足一次，第二个计数，如果测试条件每满足一次，第二个计数器就增加器就增加1。l在数字通信系统的仿真中，仿真目的是估计误比特率。测试条件就是在数字通信系统的仿真中，仿真目的是估计误比特率。测试条件就是发送给定的比特或数据符号时是否发生了差错。发送给定的比特或数据符号时是否发生了差错。l在仿真中，每当发送了一个比特或数据符号，第一计数器就增加在仿真中，每当发送了一个比特或数据符号，第一计数器就增加1。每当观。每当观测到一次差错，第二个计数器就增加测到一次差错，第二个计数器就增加1。l下一节中的两个例子会演示这

16、个方法。但我们先考察一下下一节中的两个例子会演示这个方法。但我们先考察一下AWGN信道信道的特点。的特点。16第16页，共59页，编辑于2022年，星期三9.2在通信系统中的应用在通信系统中的应用AWGN信道信道 l为了用蒙特卡罗方法估计通信系统的性能，让为了用蒙特卡罗方法估计通信系统的性能，让N个符号通过系统（实际上是个符号通过系统（实际上是系统的计算机仿真模型），并计算发送差错的个数系统的计算机仿真模型），并计算发送差错的个数Ne。如果在。如果在N次的符号次的符号发送中有发送中有Ne次差错，则符号差错概率为次差错，则符号差错概率为l这一估计是有偏的还是无偏的呢？它是否为一致估计？这一估计是

17、有偏的还是无偏的呢？它是否为一致估计？l为了在这一个尽可能简单的场景下研究这些重要的问题，我们假设信道是为了在这一个尽可能简单的场景下研究这些重要的问题，我们假设信道是AWGN（加性高斯（加性高斯白噪声）信道。白噪声）信道。l在在AWGN条件下，由信道噪声所产生的差错相互独立，而且在次符号发送中出现的条件下，由信道噪声所产生的差错相互独立，而且在次符号发送中出现的差错次数可以用二项分布来描述差错次数可以用二项分布来描述17第17页，共59页，编辑于2022年，星期三9.2.1 二项式分布二项式分布 l我们下面的任务是确定我们下面的任务是确定 的统计特性。的统计特性。l首先是确定的均值和方差。首

18、先是确定的均值和方差。l对于相互独立的差错事件，对于相互独立的差错事件，N次符号发送中有次符号发送中有Ne次差错的概率服从以下的次差错的概率服从以下的二项式分布二项式分布是二项式分布系数，是二项式分布系数，PE是单次发送的差错概率。是单次发送的差错概率。18第18页，共59页，编辑于2022年，星期三二项式分布的性质二项式分布的性质及及蒙特卡罗估计器的性质蒙特卡罗估计器的性质代入代入则差错概率的蒙特卡罗估计器的均值为则差错概率的蒙特卡罗估计器的均值为差错概率的蒙特卡罗估计器是无偏的。差错概率的蒙特卡罗估计器是无偏的。差错概率蒙特卡罗估计器的方差为差错概率蒙特卡罗估计器的方差为估计器是一致的估计

19、器是一致的 假设了二项假设了二项式分布，这式分布，这只有在差错只有在差错事件相互独事件相互独立这一条件立这一条件下才成立下才成立 19第19页，共59页，编辑于2022年，星期三讨论讨论l希望蒙特卡罗估计器无偏、一致希望蒙特卡罗估计器无偏、一致l无偏无偏平均意义下正确平均意义下正确l小的方差小的方差实用的实用的l一致一致随着随着N的增大，方差减小的增大，方差减小lPE未定，给定要求的方差，未定，给定要求的方差，N不能确定不能确定l有偏一致的呢？有偏一致的呢？不正确，除非知道如何消不正确，除非知道如何消除偏差除偏差20第20页，共59页，编辑于2022年，星期三l尽管知道估计器的性质非常重要，但

20、在许多情况下很难证明估尽管知道估计器的性质非常重要，但在许多情况下很难证明估计器是无偏的和一致的计器是无偏的和一致的l本节所得到的所有结果很理想本节所得到的所有结果很理想成立的条件是信噪声引起的差成立的条件是信噪声引起的差错相互独立，从而对应的差错概率是二项分布的错相互独立，从而对应的差错概率是二项分布的l如果差错事件是相关的，例如在带限信道的情况下，这里所给如果差错事件是相关的，例如在带限信道的情况下，这里所给出的结果不再正确，从而面临的问题就更加困难出的结果不再正确，从而面临的问题就更加困难l不过，即使差错事件不是相互独立的，式（不过，即使差错事件不是相互独立的，式（9-11）仍然是一个有

21、）仍然是一个有效的差错概率估计器。效的差错概率估计器。21第21页，共59页，编辑于2022年，星期三例例9-2l在差错事件相互独立的情况下，可用掷硬币事件来模拟二进在差错事件相互独立的情况下，可用掷硬币事件来模拟二进制发送，制发送，N个符号的发送可建模为对一枚不均匀的硬币进行个符号的发送可建模为对一枚不均匀的硬币进行N次投掷。可以假设第次投掷。可以假设第i次投掷结果为次投掷结果为“反面反面”相当于对第相当于对第i次次发送作出了正确的判决，而第发送作出了正确的判决，而第i次投掷结果为次投掷结果为“正面正面”相当相当于对第于对第i次发送作出了错误的判决。在这个例子中，掷币过程次发送作出了错误的判

22、决。在这个例子中，掷币过程所具有的统计特性由仿真来决定。因为投掷是相互独立的，所所具有的统计特性由仿真来决定。因为投掷是相互独立的，所以这个试验模拟了以这个试验模拟了AWGN信道中的二进制数据传输。信道中的二进制数据传输。22第22页，共59页，编辑于2022年，星期三Solutionl假设投掷结果为假设投掷结果为“反面反面”（没有差错）的概率为（没有差错）的概率为1-p，结果为，结果为“正面正面”（有差错）的概率为（有差错）的概率为p。我们希望通过投掷次硬币来估计。我们希望通过投掷次硬币来估计p值。值。p的蒙的蒙特卡罗估计器为特卡罗估计器为其中其中Nheads表示在连续投掷次硬币中出现表示在

23、连续投掷次硬币中出现“正面正面”的次数。当然，在的次数。当然，在N次投掷次投掷中，可以是中，可以是0到到N中的任何数，但是出现中的任何数，但是出现k次次“正面正面”的概率为的概率为因此，为了估计因此，为了估计Nheads的统计分布并确定的统计分布并确定p的估计器的估计器 ，我们必须重复，我们必须重复进行这个试验多次，比如是进行这个试验多次，比如是M次。次。23第23页，共59页，编辑于2022年，星期三Matlab 程序程序lM=2000;%number of experimentslN=500;%Number of tosses/experimentlH=zeros(1,M);%Initia

24、lize arraylH_theor=zeros(1,M);%Initialize arraylfor j=1:Ml A=rand(1,N);l heads=0;l for k=1:Nl if A(k);lsnr=10.(snrdB/10);%convert from dBlh=waitbar(0,SNR Iteration);llen_snr=length(snrdB);34第34页，共59页，编辑于2022年，星期三lfor j=1:len_snr%increment SNRl waitbar(j/len_snr)l sigma=sqrt(1/(2*snr(j);%noise standa

25、rd deviationl error_count=0;l lfor k=1:Nsymbols%simulation loop beginsld=round(rand(1);%datal if d=0l x_d=1;%direct transmitter outputl x_q=0;%quadrature transmitter outputl elsel x_d=0;%direct transmitter outputl x_q=1;%quadrature transmitter outputl end l n_d=sigma*randn(1);%direct noise component

26、l n_q=sigma*randn(1);%quadrature noise componentl y_d=x_d+n_d;%direct receiver inputl y_q=x_q+n_q;%quadrature receiver inputl if y_d y_q%test condition35第35页，共59页，编辑于2022年，星期三l d_est=0;%conditional data estimatel elseld_est=1;%conditional data estimatel endl if(d_est=d)l error_count=error_count+1;%e

27、rror counterl endl end%simulation loop endsl errors(j)=error_count;%store error count for plotlendlclose(h)lber_sim=errors/Nsymbols;%BER estimatelber_theor=q(sqrt(snr);%theoretical BERlsemilogy(snrdB,ber_theor,snrdB,ber_sim,o)laxis(snrdB_min snrdB_max 0.0001 1)lxlabel(SNR in dB)lylabel(BER)llegend(T

28、heoretical,Simulation)36第36页，共59页，编辑于2022年，星期三对于每一个对于每一个SNR值，发送值，发送Nsymbols=1000个符号，运行结果个符号，运行结果如图如图9-8所示。所示。注意，当注意，当SNR增加时，增加时，BER估计器估计器的可靠度会变差，这是由于差错发的可靠度会变差，这是由于差错发生次数减少的缘故。生次数减少的缘故。这一现象表明，仿真所发送的符这一现象表明，仿真所发送的符号数可能与号数可能与SNR有关，或者是连有关，或者是连续运行仿真程序，直到对每一个续运行仿真程序，直到对每一个SNR值都记录到相同数目的差错为值都记录到相同数目的差错为止。止

29、。37第37页，共59页，编辑于2022年，星期三例例9-4（二进制频移键控，（二进制频移键控，BFSK）l为了产生二进制为了产生二进制FSK信号空间的同相和正交分量，信号空间的同相和正交分量，令式（令式（9-24）和式（）和式（9-25）中的）中的km=/2/2，从而得，从而得到到l （9-33）l类似地有，类似地有，l （9-34）38第38页，共59页，编辑于2022年，星期三l信号空间表示如图信号空间表示如图9-9所示，其中所示，其中 1和和 2是信号是信号空间的基函数。空间的基函数。l（根据第（根据第4章的内容，对于二维空间，可以认章的内容，对于二维空间，可以认为基函数定义了低通复包

30、络信号的同相和正为基函数定义了低通复包络信号的同相和正交分量。）交分量。）l因为二进制因为二进制FSK信号空间是二维的，仿真信号空间是二维的，仿真必须产生信号和噪声的同相和正交分量。必须产生信号和噪声的同相和正交分量。l判决区域也在图判决区域也在图9-9中表示出来了中表示出来了39第39页，共59页，编辑于2022年，星期三lsnrdB_min=-3;snrdB_max=8;%SNR(in dB)limits lsnrdB=snrdB_min:1:snrdB_max;lNsymbols=input(Enter number of symbols );lsnr=10.(snrdB/10);%co

31、nvert from dBlh=waitbar(0,SNR Iteration);llen_snr=length(snrdB);Matlab program40第40页，共59页，编辑于2022年，星期三lfor j=1:len_snr%increment SNRl waitbar(j/len_snr)l sigma=sqrt(1/(2*snr(j);%noise standard deviationl lerror_count=0;l for k=1:Nsymbols%simulation loop beginsl d=round(rand(1);%datal x_d=2*d-1;%tran

32、smitter outputl n_d=sigma*randn(1);%noisel y_d=x_d+n_d;%receiver inputl if y_d 0%test conditionl d_est=1;%conditional data estimatel elsel d_est=0;%conditional data estimatel endl if(d_est=d)l error_count=error_count+1;%error counterl endl end41第41页，共59页，编辑于2022年，星期三l%simulation loop endsl errors(j)

33、=error_count;%store error count for plotlendlclose(h)lber_sim=errors/Nsymbols;%BER estimatelber_theor=q(sqrt(2*snr);%theoretical BERlsemilogy(snrdB,ber_theor,snrdB,ber_sim,o)laxis(snrdB_min snrdB_max 0.0001 1)lxlabel(SNR in dB)lylabel(BER)llegend(Theoretical,Simulation)42第42页，共59页，编辑于2022年，星期三9.3 蒙特

34、卡罗积分蒙特卡罗积分l AWGN信道，通过采样积分清除（信道，通过采样积分清除（integrate-and-dump）检测器的输出所得到的充分）检测器的输出所得到的充分统计量统计量V，是一个高斯随机变量，其均值由数据符号决定，方差由信道噪声决定。，是一个高斯随机变量，其均值由数据符号决定，方差由信道噪声决定。l在在dn=0和和dn=1的条件下，条件概率密度函数如图的条件下，条件概率密度函数如图9-11所示，其中所示，其中KT是接收机阈值。在是接收机阈值。在dn=1的条件下的条件差错概率为的条件下的条件差错概率为系统的差错概率估计系统的差错概率估计43第43页，共59页，编辑于2022年，星期三

35、蒙特卡罗积分基本概念蒙特卡罗积分基本概念 l假设要求假设要求 （9-38）l其中其中g（x）是一个在积分区域内有界的函数。由基本的概率论，可知函数的总体均）是一个在积分区域内有界的函数。由基本的概率论，可知函数的总体均值为值为l （9-39）l其中是其中是fx(x)随机变量的概率密度函数。在区间（随机变量的概率密度函数。在区间（0，1）上满足）上满足fx(x)=1，而在其他地方为零，而在其他地方为零，则则 。因此，如果。因此，如果U是在区间上均匀分布的随机变量，那么就有是在区间上均匀分布的随机变量，那么就有l （9-40）l利用相对频率的观点可以得到利用相对频率的观点可以得到l （9-41）l

36、因此，我们对被积函数进行仿真，再在区间上对它进行次采样，采样点的平均值因此，我们对被积函数进行仿真，再在区间上对它进行次采样，采样点的平均值可用来估计积分值。系统的蒙特卡罗仿真基本上也是这个思路，因为通常无法获可用来估计积分值。系统的蒙特卡罗仿真基本上也是这个思路，因为通常无法获得充分统计量在差错区域内的解析表达式，所以采用系统仿真来产生统计量的采得充分统计量在差错区域内的解析表达式，所以采用系统仿真来产生统计量的采样。样。44第44页，共59页，编辑于2022年，星期三l如果无法求式（如果无法求式（9-41）中的极限（实际应用中往往如此），）中的极限（实际应用中往往如此），则得到积分的一个近

37、似。将此近似记作得蒙特卡罗估计器为则得到积分的一个近似。将此近似记作得蒙特卡罗估计器为l （9-42）l总之，通过对函数在总之，通过对函数在N个均匀分布的采样点上求值再取平均，个均匀分布的采样点上求值再取平均，就可实现积分的估计器。这种方法适用于任何常义积分。通就可实现积分的估计器。这种方法适用于任何常义积分。通过简单的变量代换，可用蒙特卡罗方法来估算任意区间上的过简单的变量代换，可用蒙特卡罗方法来估算任意区间上的常义积分。例如，积分常义积分。例如，积分 l （9-43）l可以通过变量代换可以通过变量代换 变成标准形式并求得变成标准形式并求得l （9-44）45第45页，共59页，编辑于202

38、2年，星期三例例9-5l为了用蒙特卡罗方法估计为了用蒙特卡罗方法估计值，只需找到一个定值，只需找到一个定积分值为积分值为的函数。很容易就可想到如下积分的函数。很容易就可想到如下积分因此，我们就可以利用式（因此，我们就可以利用式（9-42）定义的算法，求解积分值，并将所得结）定义的算法，求解积分值，并将所得结果乘以果乘以4。显然我们所能做到的顶多只是选用一个很大但有限的，据此获。显然我们所能做到的顶多只是选用一个很大但有限的，据此获得的不是的准确值，而是其近似值。得的不是的准确值，而是其近似值。因此，这个积分值，也就是的估计值，是一个随机变量。如图因此，这个积分值，也就是的估计值，是一个随机变量

39、。如图9-12所示为所示为五次估计值的结果，每一次都基于五次估计值的结果，每一次都基于500个试验。个试验。(9-42)46第46页，共59页，编辑于2022年，星期三Matlab 程序程序lM=5;%Number of experimentslN=500;%Trials per experimentlu=rand(N,M);%Generate random numbersluu=1./(1+u.*u);%Define functionldata=zeros(N,M);%Initialize arrayl%The following four lines of code determinel%

40、M estimates as a function of j,0j=N.ldata(1,:)=4*uu(1,:);lfor j=2:Nl data(j,:)=4*sum(uu(1:j,:)/j;lendlest=data(N,:)%M estimates of pilest1=sum(est)/M%Average estimatelplot(data,k)%Plot resultslxlabel(Number of Trials)lylabel(Estimate of pi)47第47页，共59页，编辑于2022年，星期三48第48页，共59页，编辑于2022年，星期三收敛性收敛性 l假设要估

41、计积分值假设要估计积分值I，同时取得了，同时取得了N个随机观测或采样值个随机观测或采样值X。I的估计的估计器为器为l （9-48）l假设假设N个观测为独立同分布（个观测为独立同分布（IID）。样本的算术平均为）。样本的算术平均为l （9-49）l因此估计值因此估计值 是无偏的。由于已经假设观测的独立性，样本方差是无偏的。由于已经假设观测的独立性，样本方差为为l （9-50）l该式表明积分估计器是一致的。该式表明积分估计器是一致的。49第49页，共59页，编辑于2022年，星期三l假设假设 ，样本的方差，记作，样本的方差，记作 ，表达为，表达为 （9-51）l因此，给定被积函数因此，给定被积函数

42、 g(u)，可以确定要达到某一给定误差方，可以确定要达到某一给定误差方差所需要的差所需要的N值。值。l因为估计器是一致的，只要因为估计器是一致的，只要N充分得大，就可以得到精确的充分得大，就可以得到精确的积分估计。同时可知如果样本积分估计。同时可知如果样本g（ui）具有小方差，也可以）具有小方差，也可以得到对得到对I的精确估计。的精确估计。l对已给定的对已给定的N，如果，如果g(u)在积分区域内近似恒定（平滑），在积分区域内近似恒定（平滑），积分的蒙特卡罗估计就会非常精确。即使，只要积分的蒙特卡罗估计就会非常精确。即使，只要g(u)在积分区在积分区域内恒定，积分的估计也是精确的。域内恒定，积分

43、的估计也是精确的。50第50页，共59页，编辑于2022年，星期三置信区间置信区间 l估计器的性能常常用置信区间来表示。估计器的性能常常用置信区间来表示。l置信区间表明了估计值以概率置信区间表明了估计值以概率(1-)落在数值的给落在数值的给定范围内定范围内(i i).).l因此，其中有两个关键参数：由因此，其中有两个关键参数：由确定的概率和确定的概率和由由确定的范围。确定的范围。l置信区间定义为如下表达式置信区间定义为如下表达式51第51页，共59页，编辑于2022年，星期三置信区间。我们把区间我们把区间 称为称为 置信区间。置信区间。52第52页，共59页，编辑于2022年，星期三l现在来确

44、定参数现在来确定参数 。式（。式（9-52）可写成误差）可写成误差 的形式的形式l （9-53）l其中其中 ，这里由，这里由 式（式（9-51）确定。假设误差）确定。假设误差 是高斯随机变是高斯随机变量，则量，则 的概率密度函数近似为的概率密度函数近似为l l因为积分估计的无偏性，因为积分估计的无偏性，是零均值随机变量，所以有是零均值随机变量，所以有l （9-54）l通过变量代换通过变量代换 可得可得l （9-55）l其中其中 Q()代表高斯代表高斯Q-函数。函数。53第53页，共59页，编辑于2022年，星期三l从图从图9-13可得可得l （9-56）l从而从而l （9-57）lI的估计落在

45、区间的估计落在区间 的概率是的概率是1-，其中，其中x x由式（由式（9-51）决定。表达式）决定。表达式 的的数量大小决定了上下置信界。数量大小决定了上下置信界。l必须确定必须确定x x。54第54页，共59页，编辑于2022年，星期三 例例9-6 l为说明前面的概念，考虑以下积分为说明前面的概念，考虑以下积分l （9-58）l利用数值积分的方法，（如利用数值积分的方法，（如MATLAB函数函数quad），求得），求得I值为值为l （9-59）l用同样的方法得用同样的方法得 l （9-60）l将式（将式（9-59）和式（）和式（9-60）代入式（）代入式（9-51）可得）可得l （9-61）

46、l所以积分估计的标准偏差为所以积分估计的标准偏差为l （9-62）l上下置信界为上下置信界为l （9-63）l当然，在这个例子中存在一个严重的问题。置信区间的上下界依赖于真值，这在本例中恰好是已知的。当然，在这个例子中存在一个严重的问题。置信区间的上下界依赖于真值，这在本例中恰好是已知的。一般情况下是无法得知这个信息的，因而必须采用其他方法。常用的方法是用一个估计值来近似积分的一般情况下是无法得知这个信息的，因而必须采用其他方法。常用的方法是用一个估计值来近似积分的精确值，这个估计值通过长时间的仿真获得以保证一定的精度。精确值，这个估计值通过长时间的仿真获得以保证一定的精度。55第55页，共5

47、9页，编辑于2022年，星期三lmean=sqrt(pi)*(0.5-q(sqrt(2);%resultlint2=sqrt(pi/2)*(0.5-q(2);%2nd integrallvarx=int2-mean*mean;%variance of estimatelstdx=sqrt(varx);%standard deviationlalpha=0.1;%for 90%confidece intervall%lnsum=0;%initialize nsumlnppseg=1000;%samples per segmentlnseg=100;%number of segmentslest=

48、zeros(1,nseg);%initialize vectorl%Matlab Program56第56页，共59页，编辑于2022年，星期三lfor j=1:nseg%increment segmentl ui=rand(1,nppseg);%uniform samplesl gui=sum(exp(-ui.*ui);%integrand samplesl nsum=nsum+gui;%sum samplesl est(j)=nsum/(j*nppseg);%normalizelend%end loopl%lnn=nppseg*(1:nseg);%sample indexlub=mean+

49、stdx*qinv(alpha/2)./sqrt(nn);%upper boundllb=mean-stdx*qinv(alpha/2)./sqrt(nn);%lower boundlmeanv=mean*ones(1,nseg);%exact resultlsi=1:nseg;%seg.index for plotlplot(si,est,k-,si,meanv,k-,si,ub,k:,si,lb,k:)lxlabel(Number of Segments)%x axis labellylabel(Estimate of Integral)%y axis labelllegend(estim

50、ated value,true value,upper bound,lower bound);57第57页，共59页，编辑于2022年，星期三58第58页，共59页，编辑于2022年，星期三小结小结 l蒙特卡罗估计和仿真，建立在进行随机试验的基础上蒙特卡罗估计和仿真，建立在进行随机试验的基础上l首先确定感兴趣的事件，然后大量地重复进行相应的随机试验。感兴趣事首先确定感兴趣的事件，然后大量地重复进行相应的随机试验。感兴趣事件的相对频率为该事件发生的次数和总的随机试验次数之比，这个相对频件的相对频率为该事件发生的次数和总的随机试验次数之比，这个相对频率（为随机变量）就是感兴趣事件发生概率的估计器。