随机过程-Matlab

[e-mail]200831312008前言随机模拟若想检验数学模型是否反映客观现实，昀自然的方法是比较由模型计算的理论概率和由客观试验得到的经验频率。不幸的是，这两件事都往往是费时的、昂贵的、困难的，甚至是不可能的。此时，计算机模拟在这两方面都可以派上用场：提供理论概率的数值估计与接近现实试验的模拟。模拟的第一步自然是在计算机程序的算法中如何产生随机性。程序语言，甚至计算器，都提供了“随机”生成[0，1]区间内连续数的方法。这些数被称为伪随机数，因为每次运行程序常常生成相同的“随机数”。尽管如此，对于多数的具体问题这样的随机数已经够用。我们将假定计算机已经能够生成[0，1]上的均匀随机数。也假定这些数是独立同分布的，尽管它们常常是周期的、相关的、……。。。。。。。。。。。。本讲义的安排如下，第一章是Matlab简介，从实践动手角度了解并熟悉Matlab环境、命令，这将方便于Matlab的初学者。第二章是简单随机变量的模拟。只是给了常用的Matlab模拟语句，没有堆砌同一种变量的多种模拟方法。对于没有列举的随机变量的模拟，以及有特殊需求的读者应该由这些方法得到启发，或者参考更详细的其他文献资料。第三章是简单随机过程的模拟。主要是简单独立增量过程的模拟，多维的推广是直接的。第四章是Markov过程的模拟。包括服务系统，生灭过程、简单分支过程等。第五章包括这些模拟的应用，计算概率、估计积分、模拟现实、误差估计，以及减小方差技术。第六章给读者提供了一些经典的模拟问题，通过这些问题的时机模拟将会更加牢固地掌握实际模拟的步骤。平稳过程的模拟、以及利用平稳过程来预测的内容并没有包含在本讲义之内，但这丝毫不影响该内容的重要性，这也是将会增补进来的主要内容之一。希望读者碰到类似的问题时能够查阅相关资料解决之。各章的内容包括了模拟的基本思路和Matlab代码。源程序包展示了对各种随机过程与随机机制的有效模拟和可视化的基本技术。试图强调matlab自然处理矩阵和向量的方法。目标是为涉及应用随机模拟的读者在准备自己的程序代码时找到灵感和想法。建议读者在了解了模拟的基本方法之后就着手解决自己感兴趣的实际问题。对实际具体问题的解决有助于更深刻理解模拟的思想、也会在具体应用中拓展现有的模拟方法。第一章Matlab简介若你想在计算机上运行Matlab，点击：开始/程序/MATLAB6.5，这样将会出现有三个窗口的交互界面。如果你是初学者，可以先浏览一下Matlab的指导材料，点击：Help/MATLABHelp,打开窗口左边的“MATLAB”一节即可看到相关的内容。就自己而言，我学习Matlab更喜欢的方式是：输入并运行一些命令、观察出现的结果，然后查阅想了解的帮助文件。这也正是本节的方法。在“commandwindow”窗口中显示有提示符“”,在提示符后输入下面的命令，按回车键即可运行并显示相应的结果。当然，不要输入行号、也不必输入后面的注释。在这个部分讨论的Matlab文件有:rando.m，vrando.m，show.m。一、Matlab初步1：2*9Matlab当作计算器用。2：sin(1)Matlab仅显示四位小数，但保存的更多！3：formatlong显示更多位小数。4：sin(1)5：2^9996：formatshort7：x=sin(1)将计算结果存在变量中。8：x显示x的值。9：x=rand(10,1)x是包含有10个[0，1]上均匀分布随机数的集合，它是一个列向量或者是10×1的矩阵。10：x+5x的每个分量都加5。11：1000*xx的每个分量都乘以1000。12：x=rand(10,7)10×7的随机数矩阵。若想重复此命令或其他命令，按住向上的光标键直至看到想重复的命令。13：x=rand(1000,1)将1000换成更大的数试试。14：x=rand(1000,1);用分号“；”可以不显示结果。15：help显示标准的帮助列表。16：helpelmat显示关于初等矩阵的帮助，包括命令“rand”。17：helprand直接显示“rand”的帮助。18：x(1:20,1)取出x的第一列中的1-20个数。19：helppunctMatlab中关于标点符号的用法。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/MATLABHelp，打开左边的“MATLAB”节，选择“Functions–CategoricalList”即可。据我所知，这是寻求帮助的昀好方法。三、画出数据点33：plot(x(1:10),’*’);用“*”描出x的前10个点。注意两个单引号为英文的单引号，下同。34：plot(x-0.5);向下平移0.5，描出述据点，且将其连成线。35：holdon将下面的图形加到上面的图形中。36：plot(cumsum(x-0.5),’r’);将这个结果图画到上面的图形中。“’r’”表示用红色的线绘出，而缺省的颜色为蓝色。37：zoomon用鼠标点击可放大图形，双击回到原始的尺寸。38：clf清除当前的图形。39：z=randn(1000,1);生成1000个标准正态分布随机数。40：w=z+randn(1000,1);生成依赖z的随机数。41：plot(z,w,’*’);作出（z,w）的图形。42：axis([-33-44]);显示x在[-3，3]与y在[-4，4]范围的图形。.四、作图函数43：clf44：ezplot(’sin(x)’,[03*pi]);画出正弦函数的图形。45：holdon46：t=0:0.01:3*pi;定义一个时间点向量，间隔为0.01。47：tt为一行向量。48：t=t’现在t为一列向量。49：plot(t,sin(5*t),’r’);用红色画sin(t)关于t的函数。显然，函数ezplot不能设置图形的颜色。50：title(’sin(x)andsin(5x)’)给图形加上更恰当的标题。五、运行现有的Matlab程序51：上网下载或者拷贝一些编辑好的Matlab程序到自己的电脑中。52：如果在你电脑的某个文件夹中有现成的Matlab程序（*.m）,可以设置“CurrentDirectory”(CommandWindow窗口的上面)为该文件夹即可运行这些程序。53：如果在你电脑里的几个文件夹里都有Matlab程序，点击菜单中：File/SetPath/AddFolder,加入所有这些文件夹，昀后选择“Save”。当你在CommandWindow窗口键入一命令后，Matlab会在所有的这些文件夹中查找这个命令名。六、抛硬币54：35不等式满足结果为1。55：53结果是0。56：x=rand(20,1)前面已输入过类似的命令。输入“x=”，然后用向上的光标键往回翻看，找到后将1000改为20。57：x0.5对x的所有分量检查该不等式。58：z=1+(x0.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=[111111]/6掷色子的结果概率是一个行向量（或者1×6矩阵）。63：sum(d)确认它们的和为1！64：rando(d)用这些概率去模拟掷色子的每个结果。用向上的光标键重复这个命令几次。模拟掷色子的另一个简单的方法是放大均匀分布随机数后取整，floor(1+6*rand(1))。65：vrando([1111]/4,20)程序rando的向量版本。每个数是等概率出现的。66：show(vrando([1111]/4,100),’BGSU’)随机地生成字符B、G、S和U。出现BUGS之前，BGSU出现了吗？七、写一个Matlab程序你将创建一个新的Matlab程序，名字为mywalk.m，用它来模拟100步的随机游动。在“file”菜单下有一个空白的按钮，按下它即打开一个新的编辑窗口。在那个窗口里，分行输入下面的命令，然后保存该程序为mywalk.m。如果你保存在新的文件夹里，请确认这个文件夹是否已加入到Path中或者改变为CurrentDirectory。67：n=100;选取步数。68：x=rand(n,1);生成均匀分布随机数。69：y=2*(x0.5)-1;转换这些数到为-1和+1。70：z=cumsum(y);计算y的累积和。71：clf72：plot(z)画出z的第1,2,3,...等的值。在commandwindow窗口中输mywalk，运行（按回车）该程序，然后用光标键多次重复它。如果有错误提示，检查你的输入是否是正确的。73：运行几次后，你或许想一次就生成一个更长的字符串。到此目的的一个好的方法是将mywalk.m改为带参数的Matlabfunction，这样就可以调用它。74：在你的程序中，将行“n=100;”替换为function[z]=mywalk(n)这样，mywalk是一个带参数n的函数(生成序列的长度)，返回变量z。函数里面的变量（比如y）是内部变量，它的值不被带到函数外面。就像sin和rand一样，函数mywalk返回一个值（向量z）。回到commandwindow窗口输入：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.20.60.20]选择下一个字母。G出现后，又以另一不同的概率向量出现下一个字母，以此类推。为此，我们将创建名字为BGSU_markov.m的新程序。打开一个新的编辑窗口，输入下面的命令，然后再命令窗口输入BGSU_markov运行之。82：P=[[0.20.60.20];[00.20.60.2];[0.200.20.6];[0.60.200.2]];P是一个4×4矩阵。每一行表明将以多大的概率选择下一个字母。第一行即是数字1之后（对应字母B）的概率，第二行是数字2之后（G）的概率等等。83：x(1)=rando([1111]/4);随机地选择第一个状态。84：fori=1:399,85：x(i+1)=rando(P(x(i),:));