matlab随机过程.pdf-资料库

随机模拟讲义1 刘继成 jcliu hust@sina.com 2008 年 3 月 13 日 12008年上半年为华中科技大学数学系本科生讲授

前言随机模拟若想检验数学模型是否反映客观现实，最自然的方法是比较由模型计算的理论概率和由客观试验得到的经验频率。不幸的是，这两件事都往往是费时的、昂贵的、困难的，甚至是不可能的。此时，计算机模拟在这两方面都可以派上用场：提供理论概率的数值估计与接近现实试验的模拟。模拟的第一步自然是在计算机程序的算法中如何产生随机性。程序语言，甚至计算器，都提供了“随机”生成[0，1]区间内连续数的方法。这些数被称为伪随机数，因为每次运行程序常常生成相同的“随机数”。尽管如此，对于多数的具体问题这样的随机数已经够用。我们将假定计算机已经能够生成[0，1]上的均匀随机数。也假定这些数是独立同分布的，尽管它们常常是周期的、相关的、……。。。。。。。。。。。。本讲义的安排如下，第一章是Matlab简介，从实践动手角度了解并熟悉Matlab环境、命令，这将方便于Matlab的初学者。第二章是简单随机变量的模拟。只是给了常用的Matlab模拟语句，没有堆砌同一种变量的多种模拟方法。对于没有列举的随机变量的模拟，以及有特殊需求的读者应该由这些方法得到启发，或者参考更详细的其他文献资料。第三章是简单随机过程的模拟。主要是简单独立增量过程的模拟，多维的推广是直接的。第四章是Markov过程的模拟。包括服务系统，生灭过程、简单分支过程等。第五章包括这些模拟的应用，计算概率、估计积分、模拟现实、误差估计，以及减小方差技术。第六章给读者提供了一些经典的模拟问题，通过这些问题的时机模拟将会更加牢固地掌握实际模拟的步骤。平稳过程的模拟、以及利用平稳过程来预测的内容并没有包含在本讲义之内，但这丝毫不影响该内容的重要性，这也是将会增补进来的主要内容之一。希望读者碰到类似的问题时能够查阅相关资料解决之。各章的内容包括了模拟的基本思路和Matlab代码。源程序包展示了对各种随机过程与随机机制的有效模拟和可视化的基本技术。试图强调matlab自然处理矩阵和向量的方法。目标是为涉及应用随机模拟的读者在准备自己的程序代码时找到灵感和想法。建议读者在了解了模拟的基本方法之后就着手解决自己感兴趣的实际问题。对实际具体问题的解决有助于更深刻理解模拟的思想、也会在具体应用中拓展现有的模拟方法。

第一章 Matlab 简介若你想在计算机上运行Matlab，点击：开始/程序/MATLAB 6.5，这样将会出现有三个窗口的交互界面。如果你是初学者，可以先浏览一下Matlab的指导材料，点击：Help/ MATLAB Help,打开窗口左边的“MATLAB”一节即可看到相关的内容。就自己而言，我学习Matlab更喜欢的方式是：输入并运行一些命令、观察出现的结果，然后查阅想了解的帮助文件。这也正是本节的方法。在“command window”窗口中显示有提示符“>>”, 在提示符后输入下面的命令，按回车键即可运行并显示相应的结果。当然，不要输入行号、也不必输入后面的注释。在这个部分讨论的Matlab 文件有: rando.m，vrando.m，show.m。一、Matlab 初步 1：2*9 Matlab当作计算器用。 2：sin(1) Matlab仅显示四位小数，但保存的更多！ 3：format long 显示更多位小数。 4：sin(1) 5：2^999 6：format short 7：x=sin(1) 将计算结果存在变量中。 8：x 显示x的值。 9：x=rand(10,1) x是包含有10个[0，1] 上均匀分布随机数的集合，它是一个列向量或者是10×1的 10：x+5 x的每个分量都加5。 11：1000*x x的每个分量都乘以1000。 12：x=rand(10,7) 10×7的随机数矩阵。若想重复此命令或其他命令，按住向上的光标键直至看到矩阵。想重复的命令。 13：x=rand(1000,1) 将1000换成更大的数试试。 14：x=rand(1000,1); 用分号“；”可以不显示结果。 15：help 显示标准的帮助列表。 16：help elmat 显示关于初等矩阵的帮助，包括命令“rand”。 17：help rand 直接显示“rand”的帮助。 18：x(1:20,1) 取出x的第一列中的1-20个数。 19：help punct Matlab中关于标点符号的用法。 20：max(x) 21：mean(x) 22：sum(x) 23：median(x) x的中位数。 24：cumsum(x) x的分量累计和向量。 25：y=sort(rand(10,1)) 由小到大排序后的向量。 26：hist(x) 作出x 的直方图。 27：hist(x,30) 用30个方柱代替缺省的10个。 28：y=-log(x) 对x的分量取自然对数。 29：hist(y,30) 多数的y的分量只是接近0的，但有些是和6差不多大的，y中的数被称为指数分布随机数。 30：z=randn(1000,1); 生成1000个标准正态分布随机数。 31：hist(z,30) 直方图是钟形的。对大于 1000的数试试结果。二、获取更多帮助 32：如果你想查找不会使用的命令，可以点击：：Help/ MATLAB Help，打开左边的“MATLAB” 节，选择“Functions – Categorical List”即可。据我所知，这是寻求帮助的最好方法。

三、画出数据点 33：plot(x(1:10),’*’); 用“*”描出x的前10个点。注意两个单引号为英文的单引号，下同。 34：plot(x-0.5); 向下平移 0.5，描出述据点，且将其连成线。 35：hold on 将下面的图形加到上面的图形中。 36：plot(cumsum(x-0.5),’r’); 将这个结果图画到上面的图形中。“ ’r’ ”表示用红色的线绘出，而缺省的颜色为蓝色。 37：zoom on 用鼠标点击可放大图形，双击回到原始的尺寸。 38：clf 清除当前的图形。 39：z=randn(1000,1); 生成1000个标准正态分布随机数。 40：w=z+randn(1000,1); 生成依赖z的随机数。 41：plot(z,w,’*’); 作出（z,w ）的图形。 42：axis([-3 3 -4 4]); 显示x 在 [-3，3]与y在[-4 ，4]范围的图形。. 四、作图函数 43：clf 44：ezplot(’sin(x)’,[0 3*pi]); 画出正弦函数的图形。 45：hold on 46：t=0:0.01:3*pi; 定义一个时间点向量，间隔为0.01 。 47：t t为一行向量。 48：t=t’ 现在t为一列向量。 49：plot(t,sin(5*t),’r’); 用红色画sin(t)关于t的函数。显然，函数ezplot不能设置图形的颜色。 50：title(’sin(x) and sin(5x)’) 给图形加上更恰当的标题。五、运行现有的Matlab程序 51：上网下载或者拷贝一些编辑好的Matlab程序到自己的电脑中。 52：如果在你电脑的某个文件夹中有现成的Matlab程序（*.m）,可以设置“Current Directory” (Command Window窗口的上面)为该文件夹即可运行这些程序。 53：如果在你电脑里的几个文件夹里都有Matlab程序，点击菜单中：File/ Set Path/Add Folder, 加入所有这些文件夹，最后选择“Save”。当你在Command Window窗口键入一命令后，Matlab 会在所有的这些文件夹中查找这个命令名。六、抛硬币 54：3<5 不等式满足结果为 1。 55：5<3 结果是0。 56：x=rand(20,1) 前面已输入过类似的命令。输入“x=”，然后用向上的光标键往回翻看，找到后将1000改为20。 57：x>0.5 对x的所有分量检查该不等式。 58：z=1+(x>0.5) z的值为1或者2。这有点像抛硬币，1为正面，2为反面。 59：show(z,’正反’) 这是一个名字为show的程序，有两个变量，一个是自然数向量，一个是用来与每个数相对应显示的字符串。它是自己编制的程序，保存在：d:\MATLAB6p5\work\show.m。 60：show(1+(rand(1500,1)>0.5),’正反’) 生成1500个抛硬币的结果。现在按下向上的光标键/ 回车，就会得到很多抛硬币的结果。你找到连续出现正面最多的个数了吗？ 61：show(1+(rand(1500,1)>0.5),’O-’) 可以通过改变显示的字符来简化刚才的问题。用向上的光标键很容易更改前面的命令来实现它。这些语句对抛硬币的问题当然是足够了，因为它只有两个结果。但对其他，像掷色子，的随机试验，“rando.m”将更加有用，这也是自己编制的程序，保存在：d:\MATLAB6p5\work\show.m。 62：d=[1 1 1 1 1 1]/6 掷色子的结果概率是一个行向量（或者1 ×6矩阵）。 63：sum(d) 确认它们的和为1！ 64：rando(d) 用这些概率去模拟掷色子的每个结果。用向上的光标键重复这个命令几次。

模拟掷色子的另一个简单的方法是放大均匀分布随机数后取整，floor(1+6*rand(1))。 65：vrando([1 1 1 1]/4,20) 程序rando的向量版本。每个数是等概率出现的。 66：show(vrando([1 1 1 1]/4,100),’BGSU’) 随机地生成字符B、G、S和U。出现BUGS之前， BGSU出现了吗？七、写一个Matlab程序 86：end 87：show(x,’BGSU’); 88：hist(x,4); 你将创建一个新的Matlab程序，名字为mywalk.m，用它来模拟100步的随机游动。在“file”菜单下有一个空白的按钮，按下它即打开一个新的编辑窗口。在那个窗口里，分行输入下面的命令，然后保存该程序为mywalk.m。如果你保存在新的文件夹里，请确认这个文件夹是否已加入到Path中或者改变为Current Directory。 67：n = 100; 选取步数。 68：x = rand(n,1); 生成均匀分布随机数。 69：y = 2*(x > 0.5) - 1; 转换这些数到为-1和+1。 70：z = cumsum(y); 计算y 的累积和。 71：clf 72：plot(z) 画出z的第1, 2, 3, ...等的值。在command window窗口中输mywalk，运行（按回车）该程序，然后用光标键多次重复它。如果有错误提示，检查你的输入是否是正确的。 73：运行几次后，你或许想一次就生成一个更长的字符串。到此目的的一个好的方法是将mywalk.m 改为带参数的Matlab function，这样就可以调用它。 74：在你的程序中，将行“n = 100;”替换为 function [z] = mywalk(n) 这样，mywalk是一个带参数n的函数(生成序列的长度)，返回变量z。函数里面的变量（比如y）是内部变量，它的值不被带到函数外面。就像sin和rand一样，函数 mywalk 返回一个值（向量z）。回到command window窗口输入： 75：mywalk(1000); 运行参数为1000的程序mywalk。八、矩阵 76：M=rand(6,6) 6×6的随机数矩阵。 77：M(2,:) 取出矩阵M的第2行。 78：M(:,4) 取出矩阵M的第4列。 79：diag(M) 取出矩阵M的对角线元素。 80：sum(M) 矩阵列求和。 81：sum(M’)’ 对矩阵M的行求和。“’”表示转置九、Markov链在第66行中，序列中字母的出现是相互独立的。我们将建立下面的一种情形，B通常跟随在U之后，但决不跟在G之后。出现B后，依概率向量[0.2 0.6 0.2 0]选择下一个字母。G出现后，又以另一不同的概率向量出现下一个字母，以此类推。为此，我们将创建名字为BGSU_markov.m的新程序。打开一个新的编辑窗口，输入下面的命令，然后再命令窗口输入BGSU_markov运行之。 82：P=[[0.2 0.6 0.2 0]; [0 0.2 0.6 0.2]; [0.2 0 0.2 0.6]; [0.6 0.2 0 0.2]]; P是一个4×4矩阵。每一行表明将以多大的概率选择下一个字母。第一行即是数字1之后（对应字母B）的概率，第二行是数字2之后（G）的概率等等。 83：x(1) = rando([1 1 1 1]/4); 随机地选择第一个状态。 84：for i=1:399, 85：x(i+1) = rando(P(x(i),:)); 这是非常明智的：无论在哪个时刻，x(i)的值是多少，P(x(i),:)总是矩阵P的第x(i)行。该行的概率作为rando的参数来生成下一个状态。

第二章简单分布的模拟 Matlab里生成[0，1]上的均匀随机数的语句是：rand(1,1); rand(n,m)。一旦有了[0，1] 上均匀随机数，则我们就能够做下面的事情。在这个部分讨论的Matlab 文件有: simexp.m, simpareto.m，simparetonrm.m, simdiscr.m, simbinom.m, simgeom.m, distrmu.m, distrstat.m。一、一般连续分布（逆变换法、拒绝法、Hazard率方法）生成有连续分布函数随机数的一般方法是用反函数法。设 G(y)=F^{-1}(y)，如果 u(1)..., u(n) 是服从（0，1））上均匀分布的随机数，那么 G(u(1)), ..., G(u(n))就是分布函数为 F(x)的随机数。比如，指数分布，Pareto 分布等。 1、指数分布 simexp.m 事件以强度lambda的时间随机地发生，即事件在[t，t＋h]时间内发生的可能性是lambda×h，令t为事件发生前的等待时间。 t=-log(rand)/lambda; % 服从参数为lambda的指数分布Exp(lambda)的随机数。 t=-log(rand(1,m))./lambda; % 服从Exp(lambda)的m维行向量。 2、Pareto分布 simpareto.m 概率密度函数: f(x)=alpha/(1+x)^(1+alpha), x>0 累积分布函数: F(x) = 1-(1+x)^(-alpha)。这是带有所谓重尾分布中最简单的分布列子。产生一个均匀分布的样本，并用分布函数的反函数： sample = (1-rand(1, npoints)).^(-1/alpha)-1; 3、标准Pareto分布 simparetonrm.m 概率密度函数: f(x)=gamma*alpha/(1+gamma*x)^(1+alpha), x>0 累积分布函数: F(x) = 1-(1+gamma*x)^(-alpha) 其中，参数gamma是用来控制期望值的。在分布有重尾的情况下，若1cumprob(j)) & (uni<=cumprob(j+1))); sample(ind)=j; end 1、0-1 分布 (rand(1,m)<=p); % 生成 m 个以概率 p 为 1，概率 1-p 为 0 的随机数（m 维行向量）。三、特殊分布 1、二项分布 simbinom.m 将每次成功的概率为p的试验独立做n次，设x是成功的个数 x=sum(rand(n,m)<=p); % x是服从Bin(n,p)的m维随机数向量。

2、几何分布 simgeom.m 实验每次成功的概率为p，设x为第一次成功前失败的次数。 x=floor(-log(rand(1,m))./(-log(1-p))); %服从参数为p的几何分布Ge(p)的m维随机数行向量。floor 是取小于它的最小整数的函数。 3、泊松分布 Matlab的统计工具箱含有产生泊松分布随机数的命令，为poissrnd。 poissrnd(lambda); poissrnd(lambda, n, m); % 产生参数为lambda的泊松分布Po(lambda)随机数的n×m矩阵。如果没有上面的命令，也可以用如下的命令替代之。 arrival=cumsum(-log(rand(1,5)./lambda)); n=length(find(arrival<=lambda)); %find是找出非0值所在的位置，length是它的维数。 4、正态分布高斯分布，或正态分布的随机数用Matlab生成的命令是randn。 randn(1,m); % 服从标准正态分布N(0,1)的m维随机数行向量。 randn(n,m); % 每个分量是服从N(0,1)的n×m矩阵。 mu+sigma.*randn(1,m); % m个服从N(mu,sigma^2)分布的随机数四、离散试验的模拟 1、从{1,…,n}中任取一个。 int(n*rand(1,m)+1); 从{1,…,n}中任取不可重复两个。 a=int(n*rand(1)+1); b=int(n*rand(1)+1); while(a=b) b=int(n*rand(1)+1); end 2、随机子集模拟集合{1,…,n}的随机子集，我们是定义序列S(j)={0,1}，S(j)=1即表示将j在S中。每个S(j)， j=1,…,n,以1/2的概率独立选择0或1。 for j=1:n s(j)=int(rand(1)+1); j=j+1; end 3、随机排列假如我们向随机地排列a(1),…,a(n)，一个快速的方法是每一次互换两个数的位置，共n-1次。 for j=n:2 N=int(j*rand(1)+1); y=a(N); a(N)=a(j); a(j)=y; end 五、外部参数的随机数产生器 distrmu.m 1、在一些模拟程序中(比如更新过程)，把概率分布作为一个外部参数来传递是很方便的。这是通过创建一个MATLAB函数来实现。例如，@rand, @simpareto。分布的参数是作为数组来传递的(在rand

中为空数组{}，simpareto中为参量alpha{1.4}。 distrmu.m 是一个表-查找函数，从它的参数列表中提取期望参数的外部随机数发生器，如： mu = distrmu(@simparetonrm, {1.4, 2.5}); 2、平稳分布 distrstat.m 假设有一个分布函数为F(x)、期望值为mu的分布，则它的平稳或均衡分布的分布函数是G(x) = 1/mu * int_0^x (1-F(y))dy。例如，密度函数为f(x)=2-2*x, 0<=x<=1的线性分布是（0，1）上均匀分布上的平稳分布。参数为(alpha-1) Pareto分布是参数为alpha的Pareto分布，参数为lambda 的指数分布的平稳分布就是自己。这将出现在平稳更新过程或者排队系统的平稳版本的例子中。 distrstat.m 是一个表-查找函数，给定一个外部随机数生成器，返回它的平稳分布随机数生成器。两个参数都以数组的形式给出。至于应用，可参见平稳更新计数过程。例如： [statdist, statpar] = distrstat(@rand, {});

资料库

matlab随机过程.pdf

相关推荐

开发技术

热门标签

最新资料