matlab讲义
随机过程实验讲义
刘继成
华中科技大学数学与统计学院
前言 (1)
第一章Matlab 简介 (2)
第二章简单分布的模拟 (6)
第三章基本随机过程 (9)
第四章Markov过程 (12)
第五章模拟的应用和例子 (16)
附录各章的原程序 (51)
参考文献 (75)
若想检验数学模型是否反映客观现实，最自然的方法是比较由模型计算的理论概率和由客观试验得到的经验频率。不幸的是，这两件事都往往是费时的、昂贵的、困难的，甚至是不可能的。此时，计算机模拟在这两方面都可以派上用场：提供理论概率的数值估计与接近现实试验的模拟。
模拟的第一步自然是在计算机程序的算法中如何产生随机性。程序语言，甚至计算器，都提供了“随机”生成[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程序
你将创建一个新的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的参数来生成下一个状态。
86：end
87：show(x,’BGSU’);
88：hist(x,4);
第二章简单分布的模拟
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是用来控制期望值的。在分布有重尾的情况下，若1
sample = ((1-rand(M, N)).^(-1/alpha)-1)./gamma;
二、一般离散分布 simdiscr.m（除了上面的外，还有Alias方法）
假设给出n个概率p=[p(1)... p(n)], 满足sum(p)=1且分量p(j)是非负的。为产生m个服从这个分布的随机数，可以想象将区间（0，1）以p(1)...,p(n)为长度间隔做一个划分.产生一个均匀随机数，如果该数落在第j个间隔中，赋予此离散分布值j，重复m次。
uni=rand(1,m);
cumprob=[0 cumsum(p)];
sample=zeros(1,m);
for j=1:n
ind=find((uni>cumprob(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文件有：ranwalk.m, brownian.m, poissonti.m, poissonjp.m, ranwalk2d.m, ranwalk3d.m, bm3plot.m, poisson2d.m, poisson3d.m
一、一维情形
1、随机游动
1). 简单随机游动 ranwalk.m
“从0开始，向前跳一步的概率为p，向后跳一步的概率为1－p”
p=0.5;
y=[0 cumsum(2.*(rand(1,n-1)<=p)-1)]; % n步。
plot([0:n-1],y); %画出折线图。
2).随机步长的随机游动
选取任一零均值的分布为步长, 比如，均匀分布。
x=rand(1,n)-1/2;
y=[0 cumsum(x)-1)];
plot([0:n],y);
2、布朗运动 brownian.m
这是连续情形的对称随机游动,每个增量W(s+t)-W(s)是高斯分布N(0, t)，不相交区间上的增量是独立的。典型的模拟它方法是用离散时间的随机游动来逼近。
n=1000;
dt=1;
y=[0 cumsum(dt^0.5.*randn(1,n))]; % 标准布朗运动。
plot(0:n,y);
3、泊松过程
产生随机事件，满足: (i) 事件彼此独立发生, (ii) 两次或更多事件不会同时发生, (iii) 事件以常数强度发生。[0，t]内事件发生的次数是期望值为lambda*t的泊松分布。计数过程N(t)是泊松过程。连续两次发生的时间间隔服从参数为lambda的指数分布。
1).固定步数 poissonjp.m
%模拟n个服从Exp(lambda)的间隔时间
interarr=[0 -log(rand(1, n))./lambda];
stairs(cumsum(interarr), 0:n); %stairs画出的是水平线条。
2).固定时间区间，一个过程
固定时间区间[ 0 tmax ] 。在该区间内事件发生的总数是期望值为lambda*tmax的泊松分布。在给定事件发生次数的条件下, 事件服从该区间上的均匀分布。
%总点数是服从泊松分布的。
npoints = poissrnd(lambda*tmax);
%在点数为N的条件下，点是均匀分布的。
if (npoints>0)
arrt = [0; sort(rand(npoints, 1)*tmax);
else
arrt = 0;
end
%画出计数过程
stairs(arrt, 0:npoints);
3).固定时间区间，N个过程 poissonti.m
它被复杂化为前面算法的向量形式。到达时间间隔为指数分布的更新过程也将使用相同的算法。
%将0赋给到达时间。
tarr = zeros(1, nproc);
%将指数分布的时间间隔求和作为矩阵的列。
i = 1;
while (min(tarr(i,:))<=tmax)
tarr = [tarr; tarr(i, :)-log(rand(1, nproc))/lambda];
i = i+1;
end
%画出计数过程
stairs(tarr, 0:size(tarr, 1)-1);
二、高维情形
1、二维随机游动 ranwalk2d.m
在(u, v)坐标平面上画出点(u(k), v(k)), k=1:n, 其中(u(k))和(v(k)) 是一维随机游动。例子程序是用四种不同颜色画了同一随机游动的四条轨道。
n=100000;
colorstr=['b' 'r' 'g' 'y'];
for k=1:4
z=2.*(rand(2,n)<0.5)-1;
x=[zeros(1,2); cumsum(z')];
col=colorstr(k);
plot(x(:,1),x(:,2),col);
hold on
end
grid
2、三维随机游动 ranwalk3d.m
三维空间和上面的一样。
p=0.5;
n=10000;
colorstr=['b' 'r' 'g' 'y'];
for k=1:4
z=2.*(rand(3,n)<=p)-1;
x=[zeros(1,3); cumsum(z')];
col=colorstr(k);
plot3(x(:,1),x(:,2),x(:,3),col);
hold on
end
grid
3、三维布朗运动 bm3plot.m
npoints = 5000;
dt = 1;
bm = cumsum([zeros(1, 3); dt^0.5*randn(npoints-1, 3)]);
figure(1);
plot3(bm(:, 1), bm(:, 2), bm(:, 3), 'k');
pcol = (bm-repmat(min(bm), npoints, 1))./ ...
repmat(max(bm)-min(bm), npoints, 1);
hold on;
scatter3(bm(:, 1), bm(:, 2), bm(:, 3), ...
10, pcol, 'filled');
grid on;
hold off;
4、二维和三维空间中的泊松点 poisson2d.m, poisson3d.m
这是在空间中随机、独立地放置点的通用模型。在任何给定的空间集合中，将放置强度与其容量成比例的泊松分布的点数。在任意两个不相交的集合中的点数是独立的。
%单位体积的泊松点数强度为lambda
lambda=100;
nmb=poissrnd(lambda)
x=rand(1,nmb);
y=rand(1,nmb);
z=rand(1,nmb);
grid
scatter3(x,y,z,5,5.*rand(1,nmb));
第四章Markov过程
Markov性是随机序列x(1), x(2)... 的一种特殊形式的依赖。如果我们知道过去
x(1)...,x(n-1) 和现在x(n)，这种信息可能会或者可能不会影响未来x(n+1), x(n+1)...。Markov过程(Markov链)没有从过去获得额外的信息。对于未来的一切都由我们现在的信息所决定。随机游动和泊松过程是特殊的简单Markov过程的例子。对于那些模拟，我们已经使用了他们结构中内蕴的特殊的独立性质。
在这个部分讨论的Matlab文件有：simgeod1.m, simmm1.m, simmd1.m, simmg1.m, simmginfty.m, simstmginfty.m, birthdeath.m, moran.m, galtonwatson.m
一、离散服务系统中的缓冲动力学simgeod1.m
假设时间是离散的，并且顾客按照一个独立序列a(1), a(2)...到达服务中心，其中a(k)是在第k期到达的数量。一名顾客被服务一期(单服务系统) 。其他的顾客在一个缓冲区域等候，直到可以被服务。因此，在k时刻系统中的顾客数量为n(k)=n(k-1)+a(k)-I{n(k-1)+a(k)>=1 }, k>=2。加上初始条件n(1)=0，递归定义一个Markov链 n(k)，k>=1. 试试参数p的不同取值。从长远看，会发生什么呢?
m=200;
p=0.2;
N=zeros(1,m); %初始化缓冲区
A=geornd(1-p,1,m); %生成到达序列模型, 比如，几何分布
for n=2:m
N(n)=N(n-1)+A(n)-(N(n-1)+A(n)>=1);
end
stairs((0:m-1),N);
二、M/M/1模型simmm1.m
这是一个连续时间的单服务缓冲模型。系统的到达由强度为lambda的泊松过程决定。服务员为每位顾客的服务时间服从指数分布，均值是1/mu。由此得到的系统规模N(t), t>=0，是一个连续时间的Markov 过程，其演变如下。从N(0)=n_0开始。等待强度为lambda+mu的指数分布时间(如果n_0=0，强度为lambda), 然后以可能性lambda/(lambda+mu)向前跳跃和以可能性
mu/(lambda+mu)向后跳跃。如此循环。
在N=0时动力学的改变是在开始用一个短循环来实现的：
if i==0
mutemp=0;
else
mutemp=mu;
end
主循环仅仅检查向前跳或者向后跳：
if rand<=lambda/(lambda+mutemp)
i=i+1; % &向前跳：一个客户到达
else
i=i-1; %向后跳：一个客户离开
end
x(k)=i; %在i时刻的系统大小
有一个避免所有循环的方法，见下面的M/G/1系统。
三、M/D/1系统simmd1.m
与M/M/1一样，这个系统的到达为泊松过程，但每次服务时间是固定的长度1(例如，在缓冲环节中，固定大小的数据包的传输时间) 。这不是Markov 过程。
可以证明, 顾客离开这个系统发生在u_k=k+max(t_1,t_2-
2+1..,t_k-k+1)时刻, 其中t_1, t_2... 是泊松过程的到达时刻。因此，系统规模过程N(t), t>=0,在t_k时向前跳跃，在u_k时向后跳跃。假设我们有长度为n到达时向量t_k。这里是得到离开时间的一种方法：arrsubtr=arrtime-(0:n-1)'; % t_k-(k-1)
arrmatrix=arrsubtr*ones(1,n);
deptime=(1:n)+max(triu(arrmatrix))
现在想画出 N(t):
B=[ones(n,1) arrtime ; -ones(n,1) deptime'];
Bsort=sortrows(B,2); %按次序将跳跃分类
jumps=Bsort(:,1);
jumptimes=[0;Bsort(:,2)];
X=[0;cumsum(jumps)];
四、M/G/1系统simmg1.m
这是将M/D/1推广到一般服务时间分布S，其均值为1/mu。一个相似的计算离开时间的技术在此情形也是可以的。特别地, 取Exp(mu)分布的服务时间就是M/M/1系统。因而，我们有了模拟M/M/1的另一种方法。如果lambda/mu<1，则系统处于稳定状态。
五、M/G/infinity系统simmginfty.m
这里每个顾客得到他自己的服务。没有排队。模拟比M/G/1简单。产生到达时间加上服务时间。然后，就像上面的M/D/1系统一样，不管系统规模的改变时间为何时，总标记+1或-1。在示例代码中，将演示如何得到Pareto分布的服务时间：
alpha = 1.5; % Pareto服务时间
servtimes = rand^(-1/(alpha-1))-1; % 平稳更新过程
servtimes = [servtimes; rand(npoints-1,1).^(-1/alpha)-1];
六、M/G/infinity系统, 平稳情形, 任意服务时间simstmginfty.m (依赖 distrmu.m, distrstat.m)
当我们从时间0观察一个平稳系统时，它在负时间一直是活跃到"永远"的。因此，在时刻0，服务时间服从平稳分布的系统中有参数为(lambda*mu)的泊松分布个顾客。这在simstmginfty.m
中得以实现，它是 simmginfty.m的一个修改版本。
simstmginfty.m 也允许用一个服从服务时间分布的随机数生成器作为输入参数。参数是带一个适当的外部函数和含有分布参数数组的MATLAB函数句柄, 参见作为外部参数的RNG。
例子．
产生平稳M/G/infinity队列中[0, 5)时间内的系统规模过程，到达强度为lambda=2，服务时间服从alpha=1.6, gamma=2的标准Pareto分布。
[jmptimes,syssize]=simstmginfty(5,2,@simparetonrm,{1.6, 2},1);
stairs(jmptimes,syssize);
加入simmginfty.m得到平稳版本的步骤：
1、产生在0时刻的"负" 到达和它们的平稳服务时间。用表查找函数distrstat.m来得到平稳分布随机数生成器的句柄。
2、%在时刻0，平稳服务时间的系统中有泊松顾客数。
3、nstart=poissrnd(lambda*servmu); % 泊松随机变量
4、if (nstart>0)
5、 [statdist,statpar]=distrstat(servdist,servpar); %平稳分布句柄
6、 rndpar1={nstart,1,statpar{:}}; %随机数生成器参数
7、 stattimes=feval(statdist, rndpar1{:}); %平稳服务时间
8、arrtimes=zeros(size(stattimes)); %在t＝0时刻前到达的顾客按到达时间为0来看待。
9、end
10、一旦创建了计数过程,就删除开始时"负"到达额外的0点。增加maxtime到跳跃，以得到正确的图形。
七、生灭过程
1、一般的生灭强度 birthdeath.m
作为例子，我们选择在水平i上出生强度为lambda_i=lambda/(1+i)，死亡强度为mu_i=mu*i 的模型,其中
lambda和mu为固定的常数。要求循环满足直到下次跳跃，跳跃强度和等待时间才被更新,即
lambda_i=lambda/(1+i);
if i==0
mu_i=0;
else
mu_i=mu*i;
end
time=-log(rand)./(lambda_i+mu_i);
2、Moran模型 moran.m
另一个起源于遗传学的生灭过程。有限状态空间, 吸收界限。
x=1:N+1;
lambda(x)=(x-1).*(1-(x-1)./N); % 出生率。
mu(x)=(x-1).*(1-(x-1)./N); % 死亡率。
q(x)=lambda(x)+mu(x); % 两次跳跃时间间隔的指数分布率。
八、分支过程
Galton-Watson过程 galtonwatson.m
离散时间，生命长度为1。死亡的每个个体产生随机的后代个数。函数offspring(k)给出从人口规模为k开始的祖先向量。
p=[1/2 0 1/2];
z=[cumsum(p)];
n=length(p); % 可能的子孙数量
offmu=dot(0:n-1,p); % 子孙的平均个数
u1=sort(rand(1,k));
for j=1:n
u(j)=length(find(u1 < z(j)));
end
u=diff([0 u]);
nu=u*(0:n-1)';
九、计数过程
计数过程N(t)记录的是实值随机点过程{T_k}在区间[ 0, t)内的点数。泊松过程及排队系统中遇到的系统规模过程都是计数过程。
十、更新过程
更新过程，是一列独立同分布的正随机变量的部分和序列。这个过程可以被想象为：当同种生物的一些个体生命期结束，同时他们也被新的生命所代替的时间点序列。更新计数过程记录的是在时间区间[0，t）内更新的次数。它是一个随机阶梯函数。
第五章模拟的应用和例子
大数定律表明：1、经验均值收敛到它的期望值；2、统计物理中，轨道平均与总体平均是渐进相同的；3、为随机模拟提供了理论基础，并建立了事件频率和概率的联系。
一、计算积分I=\int_a^b f(x)dx
1、I=(b-a)E(f(a+(b-a)U))，U\sim (0,1)。模拟X_j=E(f(a+(b-
a)U_j))，用平均
1/n\sum_{j=1}^n X_j逼近I/(b-a)。大数定律说明，可以用独立试验的频率来近似期望值。
2、选取一个包含函数f图形的矩形，比如[a,b][min{f},max{f}]。生成该矩形上n对均匀分布随机数(X,Y),记录事件“Y<="" p="">
3、I=\int_a^b f(x)dx=\int_a^b \frac{f(x)}{g(x)}g(x)dx=E(\frac{f(X)}{g(X)})，然后用1/N\sum_{j=1}^N f(X_j)/g(X_j)来近似I，其中X_j为独立的密度函数为g(x)的随机数，\int_a^b g(x)dx=1。通过选择合适的函数g可以减小方差,这被称为重要样本法。
这是Monte Carlo方法的基础，它是一类计算积分的概率方法。n次近似的方差的阶是
n^{-1/2}，这比光滑函数时的梯度法差些。但作为回报，该近似对维数、被积函数的光滑性不敏感。因此Monte Carlo方法应用于不正则区域上的多重积分非常有效。另外，Monte Carlo方法的可靠性、误差的上界依赖随机数生成器的质量。在需要大量随机化的问题中使用不知道的随机数生成器是很草率的。
第一章基础准备及入门 什么是MATLAB？ MATLAB是MathWorks公司于1984年推出的数学软件，是一种用于科学工程计算的高效率的高级语言。MATLAB最初作为矩阵实验室（Matrix Laboratory），主要向用户提供一套非常完善的矩阵运算命令。随着数值运算的演变，它逐渐发展成为各种系统仿真、数字信号处理、科学可是化的通用标准语言。 在科学研究和工程应用的过程中，往往需要大量的数学计算，传统的纸笔和计算机已经不能从根本上满足海量计算的要求，一些技术人员尝试使用Basic,Fortran,C\C++等语言编写程序来减轻工作量。但编程不仅需要掌握所用语言的语法，还需要对相关算法进行深入分析，这对大多数科学工作者而言有一定的难度。与这些语言相比， MATLAB的语法更简单，更贴近人的思维方式。用MATLAB编写程序，犹如在一张演算纸上排列公式和求解问题一样高效率，因此被称为“科学便笺式”的科学工程计算语言。 MATLAB由主包和功能各异的工具箱组成，其基本数据结构是矩阵。正如其名“矩阵实验室”，MATLAB起初主要是用来进行矩阵运算。经过MathWorks 公司的不断完善，时至今日，MATLAB已经发展成为适合多学科、多工作平台的功能强大的大型软件。 本章有两个目的：一是讲述MATLAB正常运行所必须具备的基础条件；二是简明系统地介绍高度集成的Desktop操作桌面的功能和使用方法。 本章的前两节分别讲述：MATLAB的正确安装方法和MATLAB 环境的启动。因为指令窗是MATLAB最重要的操作界面，所以本章用第 1.3、1.4 两节以最简单通俗的叙述、算例讲述指令窗的基本操作方法和规则。这部分内容几乎对MATLAB各种版本都适用。 MATLAB6.x 不同于其前版本的最突出之处是：向用户提供前所未有的、成系列的交互式工作界面。了解、熟悉和掌握这些交互界面的基本功能和操作方法，将使新老用户能事半功倍地利用MATLAB去完成各种学习和研究。为此，本章特设几节用于专门介绍最常用的交互界面：历史指令窗、当前目录浏览器、工作空间浏览器、内存数组编辑器、交互界面分类目录窗、M文件编辑/调试器、及帮助导航/浏览器。 本章是根据MATLAB6.5版编写的，但大部分内容也适用于其他6.x版。 1.1M ATLAB的安装和内容选择
matlab串口通信基础讲义 ①支持基于串行接口（RS-232、RS-422、RS-485）、GPIB总线（IEEE2488、HPIB标准）、VISA总线的通信； ②通信数据支持二进制和文本（ASCII）两种方式，文本方式支持SCPI（Standard Commands for Programmable Instruments）语言； ③支持异步通信和同步通信； ④支持基于事件驱动的通信。 从以上的Matlab设备控制工具箱的特点可以看到，Matlab完全可以满足我们实现串行通信的要求。 3.1 Matlab对串行口控制的基础知识 Matlab对串行口的编程控制主要分为四个步骤。 ①创建串口设备对象并设置其属性。 scom=serial('com1');%创建串口1的设备对象scom scom.Terminator='CR';%设置终止符为CR（回车符），缺省为LF（换行符） scom.InputBufferSize=1024;%输入缓冲区为256B，缺省值为512B scom.OutputBufferSize=1024;%输出缓冲区为256B，缺省值为512B scom.Timeout=0.5;%Y设置一次读或写操作的最大完成时间为0.5s,缺省值为10s s.ReadAsyncMode='continuous'(缺省方式);%在异步通信模式方式下,读取串口数据采用连续接收数据(continuous)的缺省方式,那么下位机返回的数据会自动地存入输入缓冲区中. 注意:在些属性只有在对象没有被打开时才能改变其值,如InputBufferSize、OutputBufferSize属性等。对于一个RS-232/RS-422/RS-485串口设备对象，其属性的缺省值为波特率9 600b/s，异步方式，通信数据格式为8位数据位，无奇偶校验位，1位停止位。如果要设置的串口设置对象的属性值与缺省值的属性值相同，用户可以不用另行设置。 另外，设置串口设置对象的属性也可以用一条指令完成，如：scom=serial('COM1','BaudRate',38400,'Parity','none','DataBits',8,'StopBits',1)。也可以用set命令，如set(scom,'BaudRate',19200,'Parity','even')。创建了对象后可以在Matlab命令窗口直接敲对象名并回车，看到其基本属性和当前状态。若需要知道其全部的属性，可以用get(scom)命令。
实验1 MATLAB环境及命令窗口的使用 目的和要求 （1）熟练掌握MATLAB的启动和退出。 （2）熟练MATLAB的命令窗口。A （3）熟练常用选单和工具栏。 （4）熟悉MATLAB桌面的其他窗口。 （5）使用“帮助”查找帮助信息。 内容和步骤 学习使用MATLAB必须先熟悉MATLAB的桌面环境。MATLAB的窗口包括命令窗口（Command Window）.历史命令窗口（Command History）.当前目录浏览窗口（Current Directory Browser）.工作空间窗口（Workspace Browser）和帮助导航/浏览窗口（Help Browser）.数组编辑器窗口（Array Editor）.交互界面分类目录窗口（Launch Pad）.M文件编辑/调试器窗口（Editor/Debugger）和程序性能剖析窗口(Profiler). 1.启动MATLAB 在安装完MATLAB后，就会在Windows的桌面上出现MATLAB的图标，双击该图标，启动MATLAB；或者通过Windows的“开始”按钮，在“程序”选单中选择“MATLAB”命令来启动。 2.使用命令窗口 在命令窗口中输入以下命令并查看运行结果； >> a=2.5 a= 2.5000 >>b=[1 2;3 4] b= 1 2 2 4 >>c=’a’ c= a >>d=sin(a*b*pi/180) d= 0.0436 0.0872 0.1305 0.1736 >> e=a+c e= 99.5000 (1)单独显示命令窗口。选择选单”desktop”-“desktop Layout”-“Command window only”
前言 MATLAB是电子信息工程、通信工程等专业的专业基础课程，是一个功能十分强大的数学应用软件，能够快速处理大量复杂的数学计算，如求矩阵的逆、矩阵的特征向量等等。学生熟练掌握MA TLAB，将能为后继课程的学习提供很好的计算工具和仿真平台。在经过全面的训练后，学生基本掌握MATLAB基本语法和基本函数的用法，利用MATLAB这门工具语言联系以前所学知识，突破数学计算方面的障碍，更好地理解本专业课程的基本概念、基本原理；本书让学生初步掌握MA TLAB的工具箱SIMULINK的使用，为后继课程提供方便；能根据需要选择参考书，查阅手册，通过独立思考，深入钻研有关问题，学会自己独立分析问题、解决问题，具有一定的创新能力。 本书根据理论课程的设置共编写了8个实验，内容丰富，涵盖了各个方面上机练习，同学们在上机时一定要大量的练习本书内容，才可以掌握相关知识。
实验一认识MATLAB6.0/7.0 一、实验目的 1.掌握MATLAB开发环境的启动和退出 2.了解MATLAB开发环境的基本组成 3.学会在MA TLAB开发环境中编辑简单程序 4. 了解MA TLAB6.0“在线帮助” 二、实验内容 （一）MATLAB开发环境 1.启动并熟悉MATLAB开发环境。 2.熟悉MA TLAB开发环境中操作桌面的组成:命令窗口、启动平台窗口、工作空间窗口、 命令历史窗口等。如图(1-1) 图(1-1) 3.练习MA TLAB开发环境中各窗口功能。（课本中涉及的内容逐一练习）。 例：①通过MATLAB参数设置使得编辑程序时关键词(Keyword)为红色、注释（Comment）为绿色等（随意编程者喜欢的颜色）如图（1-2）
chap6 习题 1. 请分别写出用for 和while 循环语句计算100000021000000 02.02.02.012 .0+++==∑=Λi i K 的 程序。此外，还请写出避免循环的数值、符号计算程序。 s1 = 1.2500 2. 编写一个函数M 文件，它的功能：没有输入量时，画出单位圆（见图p6-1）；输入量是 大于2的自然数N 时，绘制正N 边形，图名应反映显示多边形的真实边数（见图p6-2）；输入量是“非自然数”时，给出“出错提示”。此外，函数M 文件应有H1行、帮助说明和程序编写人姓名。 图 p6-1
图 p6-2 function prob_solve602(n) % prob_solve602(n) plot a circle or a polygon with n edges % prob_solve602 plot a circle % n 应为大于2的自然数 % By ZZY, 2006-2-15 if nargin==0 t=0:pi/100:2*pi; x=exp(i*t); str='Circle'; else if (nargin~=0)&(n<=2) error('输入量应是大于2的自然数') end ; if n-round(n)~=0 %检查非自然数 error('输入量应是大于2的自然数') end ; t=(0:n)/n*2*pi; x=exp(i*t); str=['Polygon with ', int2str(n),' edges']; % 合成字符串 end plot(real(x),imag(x),'r','LineWidth',4) title(str) axis square image off shg 3. 用泛函指令fminbnd 求|]sin[cos |)(x e x y x --=在x=0附近的极小值。fminbnd 的第 一个输入量要求使用匿名函数表达。 x = -0.8634 4. 在matlab 的 \toolbox\matlab\elmat\private 文件夹上有一个“烟圈矩阵”发生函数 smoke.m 。运行指令smoke(3,0,'double') ，将生成一个3阶伪特征根矩阵如下 A = -0.5000 + 0.8660i 1.0000 0 0 -0.5000 - 0.8660i 1.0000 1.0000 0 1.0000 现在的问题是：在MA TLAB 当前目录为\work 情况下，如何利用函数句柄调用smoke.m 函数，产生3阶伪特征根矩阵。请写出相应的程序或操作步骤。
MATLAB讲义 第一章MATLAB系统概述 1.1 MATLAB系统概述 MATLAB（MATrix LABoratory）矩阵实验室的缩写，全部用C语言编写。 特点： （1）以复数矩阵作为基本编程单元，矩阵运算如同其它高级语言中的语言变量操作一样方便，而且矩阵无需定义即可采用。 （2）语句书写简单。 （3）语句功能强大。 （4）有丰富的图形功能。如plot,plot3语句等。 （5）提供了许多面向应用问题求解的工具箱函数。目前，有20多个工具箱函数，如信号处理、图像处理、控制系统、系统识别、最优化、神经网络的模糊系统等。 （6）易扩充。 1.2 MATLAB系统组成 （1）MATLAB语言 MATLAB语言是高级的矩阵、矢量语言，具有控制流向语句、函数、数据结构、输入输出等功能。同时MATLAB又具有面向对象编程特色。MATLAB语言包括运算符和特殊字符、编程语言结构、字符串、文件输入/输出、时间和日期、数据类型和结构等部分。 （2）开发环境 MATLAB开发环境有一系列的工具和功能体，其中大部分具有图形用户界面，包括MATLAB桌面、命令窗口、命令历史窗口、帮助游览器、工作空间、文件和搜索路径等。 （3）图形处理 图形处理包括二维、三维数据可视化，图像处理、模拟、图形表示等图形命令。还包括低级的图形命令，供用户自由制作、控制图形特性之用。 （4）数学函数库 有求和、正弦、余弦等基本函数到矩阵求逆、求矩阵特征值和特征矢量等。 MATLAB数学函数库可分为基本矩阵和操作、基本数学函数、特殊化数学函数、线性矩阵函数、数学分析和付里叶变换、多项式和二重函数等。 （5）MATLAB应用程序接口（API） MATLAB程序可以和C/C++语言及FORTRAN程序结合起来，可将以前编写的C/C++、FORTRAN语言程序移植到MATLAB中。 1.3 MATLAB的应用范围 MATLAB的典型应用包括： ●数学计算 ●算法开发 ●建模、仿真和演算 ●数据分析和可视化 ●科学与工程绘图 ●应用开发（包括建立图形用户界面） 以矩阵为基本对象
2.3终值及其应用 2.3.1终值的概念 终值是与现值相对的概念，是指当前的一项现金流在未来某个时刻的价值。在求终值问题时应该考虑单利和复利的问题，一般如果没有特别的说明则都是按照复利（离散复利）进行计算。 在复利计息的情况下，当前的现金流PV在利率为r时到第t期期末的终值为： t FV) = 1(+ r PV 2.3.2终值的计算 在Matlab中，用来计算现金流的终值的函数有fvfix和fvvar两个。同样，-fix函数用来计算规则现金流的终值；而-var函数则用来计算不规则现金流的终值。 【例2.9】一投资者的储蓄账户初始余额为$1500，在随后的10年中，每月末都会收到$200并存入该账户，银行的年利率为9%。试计算其到期时的价值。 通过执行fvfix函数命令： FutureVal = fvfix(Rate, NumPeriods, Payment, PresentVal, Due) 即可计算出该固定收入现金流的的终值。 变量解释： Rate：周期性收支的利息率，以小数的形式输入； NumPeriods：周期性收支的次数； Payment：每期收支的现金流数额； PresentVal：初始余额 Due：收支被预定或确定的时间：0表示在期末收支（默认值），1表示在期初收支（任选）。 输入命令： >>FutureVal = fvfix(0.09/12, 12*10, 200, 1500, 0) 输出结果： FutureVal = 42379.89 即该现金流到期时的价值为42379.89$。 【例2.10】设某投资者期初投资为$10,000，在随后的5年投资期中每年产生的收入流依次为$2000、$1500、$3000、$3800、$5000，年利率为8%。试计算该现金流到期时的价值。 通过执行fvvar函数命令： FutureVal = fvvar(CashFlow, Rate, IrrCFDates) 即可求出这个规则（周期性的）现金流的终值。 输入命令： >>FutureVal = fvvar([-10000 2000 1500 3000 3800 5000], 0.08) 输出结果： FutureVal = 2520.47 即该现金流到期时的价值为2520.47$。 如果期初投资的$10,000产生的是一个不规则的现金流（如下所示），则计算时要将期初的投资和各个现金流发生的日期也考虑进去。利率为9%。
第二章 序列的付氏变换（DTFT ）和Z 变换 一、序列的DTFT ：X （e jw ） 序列x(n)的付氏变换如下： ∑+∞ ∞--=n j j e n x e X ωω )()( ∞≤≤∞-ω ωπ ωπ π ωd e e X n x n j j ?- = )(21)( 方法一:定义法 此法只为让大家熟悉DTFT 形成的过程 因为MATLAB 无法计算连续变量w,只能在πωπ<≤-（或考虑对称性，在0<≤-ωπ）范围内，把w 赋值为很密的、长度很长的向量来近似连续变量。通常最简单的就是赋以K 个等间隔的值： K k d k π ωω2=?= 1:0-=K k 它表示把基本数字频率的范围π2均分成K 份后，每一份的大小。 则频率向量表示为w=[w1,w2,…,wk]=k ?dw, 序列的位置向量为n=[n1,n2,…,nN] 则DTFT 的计算式可用一个向量与矩阵相乘的运算来实现： [X(w 1),X(w 2),…,X(w k )] =[x(n 1),x(n 2),…,x(n N )]?????? ?????????N k N N k k n j n j n j n j n j n j n j n j n j e e e e e e e e e ωωωωωωωωω 212222111211 =[x(n 1),x(n 2),…,x(n N )]? w *jn T e 又因为 (jn T *w)=jn T *kdw=j*dw*n ’*k 所以
Xw=xn*exp(j*dw*n ’*k) 例1：求有限长序列xn=[1 3 5 3 1]，n=-1:3的DTFT,画出它在w=-8~8rad/s 范围内的频率特性，讨论其对称性。再把xn 左右移动，讨论时移对DTFT 的影响。 解:根据DTFT 定义得 ωωωωωω 32353)()(j j f f n n j j e e e e e n x e X ---∞ -∞ =-++++== ∑ 88≤≤-ω 将w 在-8~8rad/s 之间分为1000份。 >> xn=[1 3 5 3 1];nx=-1:3; >> w=linspace(-8,8,1000);%设定频率向量? >> Xw=xn*exp(-j*nx'*w);%用定义计算DTFT >> subplot(3,5,1);stem(nx,xn);axis([-2,6,0,6]);title('?-ê?DòáD');ylabel('xn');%画序列图 >> subplot(3,5,2);plot(w,abs(Xw));title('·ù?è');% 画幅频曲线 >> subplot(3,5,3);plot(w,angle(Xw));title('?à??')% 画相频曲线 >> subplot(3,5,4);plot(w,real(Xw));title('êμ2?')% 画实频曲线 >> subplot(3,5,5);plot(w,imag(Xw));title('Dé2?')% 画虚频曲线 >> nx_2=nx+2;%使xn 右移2位 >> Xw_2=xn*exp(-j*nx_2'*w); >> subplot(3,5,6);stem(nx_2,xn);axis([-2,6,0,6]);ylabel('xn_2');title('óòò?2??') >> subplot(3,5,7);plot(w,abs(Xw_2)); >> subplot(3,5,8);plot(w,angle(Xw_2)); >> subplot(3,5,9);plot(w,real(Xw_2)); >> subplot(3,5,10);plot(w,imag(Xw_2)); >> nx_1=nx-1;%使xn 左移1位 >> Xw_1=xn*exp(-j*nx_1'*w); >> subplot(3,5,11);stem(nx_1,xn);axis([-2,6,0,6]);ylabel('xn_1');title('×óò?1??') >> subplot(3,5,12);plot(w,abs(Xw_1)); >> subplot(3,5,13);plot(w,angle(Xw_1)); >> subplot(3,5,14);plot(w,real(Xw_1)); >> subplot(3,5,15);plot(w,imag(Xw_1));
第一章 Matlab 基本介绍 一、数学建模常用软件简介 数值计算 Matlab 符号计算 Maple , Mathematica 统计软件 SPSS, SAS 优化软件 LINGO OFFICE 软件 Word , Excel 二、matlab 界面介绍 1、command window （命令窗口） 2、wordspace （工作空间） 3、command history （历史命令窗口） 4、菜单 (1)File->import data (2)View->desklayout->default (3)Help 三、一些常用命令 1、clc (清空命令窗口) 2、clear （清空工作空间变量） 3、save （保存工作空间中变量到指定文件） 4、load （导入文件中数据）（注：双击数据文件也可） 5、help （帮助） 6、doc （查询帮助文档） 第二章 数值计算 一、数据类型 1、主要四大类数据类型：数值型，字符串，符号型（代数式），逻辑型 1???????? 、浮点数值型长整2、整型短整 字符串 name=’lisan ’ a=’x ’ 符号型 用syms ， sym 定义 逻辑型 取值只能为0或1，即真或假 2 常用运算符 数值运算：+ , - ,* , / , \, ^, .*, ./, .^ 关系运算 （运算结果为逻辑型，即0或1） >,<,>=,<=,==,~=
逻辑运算（运算结果为逻辑型，即0或1） 与或非，&，|，~ any, all 基本数学函数 三角sin ,asin,cos,acos,tan,atan,cot,acot 指数exp,log,log2,log10,sqrt 其他abs,real,imag,sign,mod,floor,ceil 特殊字符, ; . : [] () @ % 2、变量命名规则 （1）以字母开头，可包含字母、数字、下划线，不超过31位字符。 （2）区分大小写。 3、常量 i, j 虚数单位 pi , 圆周率 eps, inf NaN 4、数字的输入输出格式 format 格式参数 short long rat 5、字符串 （1）字符串生成 name=’lisan’ a=char(‘l’,’i’) size(name) 查看字符串长度 length(name) 查看字符串长度 （2）字符与数组之间的转换 double 字符转换为ASC码 num2str 数字转化为字符 str2num 字符转位数字 a=’2’ b=a*2 b=double(a)*2 b=str2num(a)*2 （3）字符串操作相关函数。 5、结构体 定义方式（1）利用struct函数 （2）直接定义 https://www.doczj.com/doc/2419357717.html,=’marry’ a.length=170
通信原理实验的MATLAB仿真讲义 MATLAB原意为“矩阵实验室—MA-TrixLABoratory”，它是目前操纵界国际上最流行的软件，它除了传统的交互式编程之外，还提供了丰富可靠的矩阵运算、图形绘制、数据与图象处理、Windows编程等便利工具。MATLAB还配备了大量工具箱，特别是还提供了仿真工具软件SIMULINK。SIMULINK这一名字比较直观地说明了此软件的两个显著的功能：SIMU(仿真)与LINK(连接)，亦即能够利用鼠标在模型窗口上“画”出所需的系统模型，然后利用SIMULINK提供的功能对操纵系统进行仿真与线性化分析。 MATLAB在80年代一出现，首先是在操纵界得到研究人员的瞩目。随着MA T-LAB软件的不断完善,特别是仿真工具SIMULINK的出现，使MA TLAB的应用范围越来越广。 MATLAB的仿真环境(simulink)提供的系统模型库包含下列几个子模型库：Sources(输入源)、Sinks(输出源)、Discrete(离散时间系统)、Linear(线性环节)、Non-linear(非线性环节)、Connections(连接及接口)、Extras(其它环节)。打开子模型库,你会发现每个模型库都包含许多个子模块，比如Sources模型库里含有阶跃函数、正弦函数、白噪声函数、MATLAB空间变量、信号发生器等子模块。另外在Extras子模型库下还有一个BlockLibrary，集中了子模型库中最常用及其它常用的子模块,使用起来特别方便。通信系统通常都能够建立数学模型,在数学模型中,要紧包含乘法器、加法器、信号发生器、滤波器等，而这些在上述的simulink 系统模型库中通常都可找到，关于没有的模块(如伪随机信号发生器),可自己根据掌握的技术生成所需的子模块，随时调用。这样就可根据数学模型，建立通信系统的仿真模型。 应用MA TLAB下的SIMULINK仿真工具能够很方便地进行通信系统仿真，利用SIMULINK仿真工具下的现有子模块进行仿真。 实验一.利用SIMULINK仿真常规调幅AM 滤波调制与卷积定理 从信号与线性系统分析观点看，滤波如图 1 是系统的冲激响应 h(t)对输入信号 x(t)的卷积作用，即 y(t) = x(t) * h(t) (1) 对应的频域分析是 Y( ω) = X(ω )H(ω ) (2) 即时域卷积处理对应于频域内相乘，（1）与（2）式是时域卷积定理 图 1 滤波图 2 调制 再看调制(包含解调、混频等)如图2, 是两信号相乘即 y(t) = x(t)c(t) (3) 对应的频域分析是 Y(ω)=1/2π[X(ω)*C(ω)](4)
天津大学matlab讲义 第2章 MATLAB程序设计 MATLAB语言为解释型程序设计语言。在程序中可以出现顺序、选择、循环三种基本控制 结构，也可以出现对M-文件的调用(相当于对外部过程的调用)。 由于 MATLAB开始是用FORTRAN语言编写、后来用 C语言重写的，故其既有FORTRAN的 特征，又在许多语言规则方面与C语言相同。 2.1 顺序结构语句 在顺序结构语句中，包括表达式语句、赋值语句、输入输出语句、空语句等。 2.1.1 表达式语句格式： 表达式，％显示表达式值表达式；％不显示表达式值表达式％显示表达式值如： x + y, sin(x); –5 最后的表达式值暂保存在变量ans中。 2.1.2 赋值语句格式： v = 表达式，％结果送v并显示v v = 表达式；％结果送v不显示v v = 表达式％结果送v并显示v 2.1.3 空语句格式： ，； 2.1.4 输入语句 1、input语句（实际上是函数）格式1： input（提示字符串）功能： 显示提示字符串，可输入数字、字符串（两端用单引号括起）、或表达式格式2： input（提示字符串，'s'）功能： 显示提示字符串，并把输入视为字符串 2、yesinput语句格式： yesinput(提示字符串，缺省值，值范围) 功能：
显示提示字符串和缺省值，若只打入回车则以缺省值作为输入值，若输入的值不在指定 范围内则认为输入无效，B并等待用户重新输入。 如： t=yesinput('指定线的颜色',… 'red','red|blue|green') 运行结果如下：指定线的颜色(red):yellow %不在值内指定线的颜色(red):blue %重输 t = blue x=yesinput('输入元素个数',10,[1,20]) 运行结果如下：输入元素个数(10): x = 10 3、Keyboard语句格式： Keyboard 功能： 暂停M-文件的执行，并等待用户从键盘输入命令以查看或改变变量的值，直 到输入 return命令而返回相应的M-文件继续执行。本语句用于调试M-文件。 4、pause语句格式1： pause 功能：暂停，敲下任一键继续格式2： pause(n) 功能：暂停n秒格式3： puase on 功能：本命令后的pause语句有效格式4： pause off 功能：本命令后的pause语句无效 5、menu语句格式： menu('菜单名',S1,S2,…,Sk) 功能：生成一个按钮式菜单系统其中: 字符串S1,S2,…,Sk为菜单项(K≤32)。函数返回值为用户选中的菜单项号。如M-文件，menu_d.m如下： %选择一种颜色 while 1 k=menu('选择一种颜色','红色',… '黄色','兰色','绿色','白色','关闭'); if k = = 1 color = 'Red' elseif k = = 2 color = 'Yellow' elseif k = = 3 color = 'Blue' elseif k = = 4 color = 'Green' elseif k = = 5 color = 'White' elseif k = = 6 break end end
matlab讲义 随机过程实验讲义 刘继成 华中科技大学数学与统计学院 前言 (1) 第一章Matlab 简介 (2) 第二章简单分布的模拟 (6) 第三章基本随机过程 (9) 第四章Markov过程 (12) 第五章模拟的应用和例子 (16) 附录各章的原程序 (51) 参考文献 (75) 若想检验数学模型是否反映客观现实，最自然的方法是比较由模型计算的理论概率和由客观试验得到的经验频率。不幸的是，这两件事都往往是费时的、昂贵的、困难的，甚至是不可能的。此时，计算机模拟在这两方面都可以派上用场：提供理论概率的数值估计与接近现实试验的模拟。 模拟的第一步自然是在计算机程序的算法中如何产生随机性。程序语言，甚至计算器，都提供了“随机”生成[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仅显示四位小
Matlab程序设计及应用实验讲义 自编 电子科学与工程系
实验一 MATLAB 环境与命令窗口 一、实验目的 1）熟悉MATLAB 的操作环境及大体操作方式； 2）把握MATLAB 的搜索途径及其设置方式； 3）熟悉MATLAB 帮忙信息的查阅方式； 二、实验要紧仪器设备和材料 运算机PC 一台 2020a 软件 3、实验内容和原理 一、先成立自己的工作目录，再将自己的工作目录设置到MATLAB 的搜索途径下，再实验用help 命令可否查询到自己的工作目录。 二、在MATLAB 环境下验证下面几个例子，并总结MATLAB 的要紧优势。 1）绘制正弦曲线和余弦曲线 2）求方程432379230x x x ++-=的全数根 3）求积分()1 0ln 1x x dx +⎰ 4）求解线性方程组234832245917x y z x y z x y z -+=⎧⎪++=⎨⎪+-=⎩ 3、利用MATLAB 的帮忙功能别离查询inv 、plot 、 max 、round 等函数的功能与用法 4、完成以下操作： 1）在MATLAB 命令窗口输入以下命令： x=0:pi/10:2*pi;
y=sin(x); 2）在工作空间窗口选择变量y，再在工作空间窗口选择画图菜单命令或在工具栏中单击画图命令按钮，绘制变量的图形，并分析图形的含义。 五、访问MathsWorks公司的主页，查询有关MATLAB的产品信息。 试探与练习 一、如何启动和退出MATLAB的集成环境？ 二、简述MATLAB的要紧功能。 3、若是一个MATLAB命令包括的字符很多，需要分成多行输入，该如何处置？ 4、help命令和look for命令有何区别？ 五、在MATLAB环境下，成立一个变量fac,同时又在当前目录下成立了一个M 文件，若是需要运行文件，该如何处置？
第1章方程求解 本章学习目的： 1．复习求解方程的基本原理和方法，掌握解方程的迭代算法； 2．能利用MATLAB软件编写迭代算法程序，了解迭代过程的图形表示； 3．熟练掌握用MATLAB软件的函数来求解方程和方程组； 4．通过范例展现求解实际问题的初步建模过程和MATLAB程序设计。 §1.1 引言 “方程是很多工程和科学工作的发动机”。研究大型的土建结构、机械结构、输电网络、管道网络，研究经济规划、人口增长、种群繁殖等问题时，简单的分析可以直接归结为线性或非线性方程组，复杂一些要用到（偏）微分方程，求数值解时将转化为n非常大的方程组。若干世纪以来，工程师和科学家花了大量的时间用于求解方程（组），数学家研究各种各样的方程求解方法。 本章我们就是要学习求解线性方程组、非线性方程（组）的方法，以及利用数学软件利用计算机对方程和方程组进行求解。 §1.2 方程的求解方法 考虑求方程f(x)=0的解，我们通常采用这样的几种方法：因式分解法、图形放大法、数值迭代逼近法。
1．因式分解法： 这是我们最熟悉、常用的一种方法，这个方法的关键在分解因式，包括对多项式函数、三角函数和指数函数等的分解。但对于无法进行分解的函数则无能为力。 2．图形放大法： 由于计算机的广泛应用，可以非常方便地作出函数f(x)的图形（曲线），找出曲线与x轴的交点的横坐标值，就可求出f(x)＝0的近似根。这些值尽管不精确，但是直观，方程有多少个根、在什么范围，一目了然。并且可以借助于计算机使用图形局部放大功能，将根定位得更加准确。 3．数值迭代逼近法： 利用图形的方法或连续函数的零点存在性定理，可以推知f(x)在某一区间内有根，我们就可以用数值方法来求方程的根，这就是迭代逼近法。 迭代逼近法分为区间的迭代和点的迭代。 区间迭代又分为对分法和黄金分割法；点的迭代又分为简单迭代法、单点割线法、两点割线法、牛顿法等。迭代失败后又可以采用加速迭代收敛方法。各种迭代方法原理都很简单（数值分析课有详细介绍），请同学们自学。 1．2．1图形放大法 用图形放大法求解方程f(x)＝0的步骤： （1）建立坐标系，作曲线f(x)； （2）观察f(x)与x轴的交点； （3）将其中的一个交点进行局部放大； （4）该交点的横坐标值就是方程的一个根； （5）对所有的交点进行相同的处理，就得到方程的所有解。
《Matlab仿真应用详解》 一、基本概念 １.1、什么是计算机仿真 1、仿真定义 基本思想：仿真的基本思想是利用物理的或数学的模型来类比模仿现实过程，以寻求过程和规律。它的基础是相似现象，相似性一般表现为两类：几何相似性和数学相似性。当两个系统的数学方程相似，只是符号变换或物理含义不同时，这两个系统被称为“数学同构”。 仿真的方法可以分为三类： （1）实物仿真。它是对实际行为和过程进行仿真，早期的仿真大多属于这一类。物理仿真的优点是直观、形象,至今在航天、建筑、船舶和汽车等许多工业系统的实验研究中心仍然可以见到.比如：用沙盘仿真作战,利用风洞对导弹或飞机的模型进行空气动力学实验、用图纸和模型模拟建筑群等都是物理仿真.但是要为系统构造一套物理模型，不是一件简单的事，尤其是十分复杂的系统,将耗费很大的投资,周期也很长。此外，在物理模型上做实验,很难改变系统参数,改变系统结构也比较困难。至于复杂的社会、经济系统和生态系统就更无法用实物来做实验了. (2）数学仿真。就是用数学的语言、方法去近似地刻画实际问题，这种刻画的数学表述就是一个数学模型。从某种意义上,欧几里德几何、牛顿运动定律和微积分都是对客观世界的数学仿真.数学仿真把研究对象(系统）的主要特征或输入、输出关系抽象成一种数学表达式来进行研究。数学模型可分为: ●解析模型(用公式、方程反映系统过程); ●统计模型（蒙特卡罗方法）；（一种基于随机数的计算方法) ●表上作业演练模型。（用列表的方法求解线性规划问题中运输模型的计算方法。是指线性规划一种求解方法。当某些线性规划问题采用图上作业法难以进行直观求解时,就可以将各元素列成相关表，作为初始方案，然后采用检验数来验证这个方案，否则就要采用闭回路法、位势法或矩形法等方法进行调整，直至得到满意的结果。这种列表求解方法就是表上作业法。) 然而数学仿真也面临一些问题,主要表现在以下几个方面： ●现实问题可能无法用数学模型来表达，即刻画实际问题的表达式不存在或找不到； ●找到的数学模型由于太复杂而无法求解； ●求出的解不正确，可能是由模型的不正确或过多的简化近似导致的. （3)混合仿真.又称为数学—物理仿真，或半实物仿真，就是把物理模型和数学模型以及实物联合在一起进行实验的方法，这样往往可以获得较好的效果。 2、计算机仿真 计算机仿真也称为计算机模拟，就是利用计算机对所研究系统的结构、功能和行为以及参与系统控制的主动者——人的思维过程和行为，进行动态性的比较和模仿，利用建立的仿真模型对系统进行研究和分析,并可将系统过程演示出来. 1．2计算机仿真模型与方法 １、系统 系统是指相互联系又相互作用的元素之间的有机组合.这里所指的系统是广义的,它包含所有的工程系统和非工程系统.电气、机械和通信系统都是工程系统，而经济、交通、管理和生物系统等都是非工程系统。 任何系统都存在三方面需要研究的内容: 实体：组成系统的具体对象. 属性：实体的特性（状态和参数)。即实体、属性和活动。 由于组成系统的实体之间相互作用而引起实体属性的变化,通常用“状态"的概念来描述。研究系统就是研究系统状态的改变，即系统的转变。 系统具有四个特性: （1）目的性（2)集合性（３)相关性(４）环境适应性
MATLAB软件基础 §1MATLAB 概述 MATLAB 是MATrix LABoratory〔“矩阵实验室〞〕的缩写，是由美国MathWorks 公司开发的集数值计算、符号计算和图形可视化三大根本功能于一体的，功能强大、操作简单的语言。是国际公认的优秀数学应用软件之一。 20世纪80年代初期，Cleve Moler与John Little等利用C语言开发了新一代的MATLAB语言，此时的MATLAB语言已同时具备了数值计算功能和简单的图形处理功能。1984年，Cleve Moler与John Little等正式成立了Mathworks公司，把MA TLAB 语言推向市场，并开始了对MATLAB工具箱等的开发设计。1993年，Mathworks公司推出了基于个人计算机的MATLAB 4.0版本，到了1997年又推出了MATLAB 5.X版本〔Release 11〕，并在2000年又推出了最新的MATLAB 6版本〔Release 12〕。 现在，MATLAB已经开展成为适合多学科的大型软件，在世界各高校，MATLAB已经成为线性代数、数值分析、数理统计、优化方法、自动控制、数字信号处理、动态系统仿真等高级课程的根本教学工具。特别是最近几年，MATLAB在我国大学生数学
建模竞赛中的应用，为参赛者在有限的时间内准确、有效的解决问题提供了有力的保证。 概括地讲，整个MATLAB系统由两局部组成，即MATLAB 内核及辅助工具箱，两者的调用构成了MATLAB的强大功能。MATLAB语言以数组为根本数据单位，包括控制流语句、函数、数据结构、输入输出及面向对象等特点的高级语言，它具有以下主要特点： 1〕运算符和库函数极其丰富，语言简洁，编程效率高，MATLAB 除了提供和C语言一样的运算符号外，还提供广泛的矩阵和向量运算符。利用其运算符号和库函数可使其程序相当简短，两三行语句就可实现几十行甚至几百行C或FORTRAN的程序功能。 2〕既具有结构化的控制语句〔如for循环、while循环、break 语句、if语句和switch语句〕，又有面向对象的编程特性。 3〕图形功能强大。它既包括对二维和三维数据可视化、图像处理、动画制作等高层次的绘图命令，也包括可以修改图形及编制完整图形界面的、低层次的绘图命令。 4〕功能强大的工具箱。工具箱可分为两类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能、图
连续时间信号的分析 一、实验目的 1.学习使用MATLAB 产生基本的连续信号、绘制信号波形。 2.实现信号的基本运算，为信号分析和系统设计奠定基础。 二、实验原理 1、基本信号的产生 时间间隔代替连续信号。 连续指数信号的产生 连续矩形脉冲信号(门信号)的产生。 连续周期矩形波信号的产生。 2、信号的基本运算 相加、相减、相乘、平移、反折、尺度变换。 三、实验内容 1. 用MATLAB 编程产生正弦信号()sin(2),2,5Hz,3 f t K ft K f π πθθ=+===，并画图。 代码如下： clc clear f0=5; w0=2*pi*f0; t=0:0.001:1; x=2*sin(w0*t+pi/3); plot(t,x) title('正弦信号')
正弦信号 2. 用MATLAB 编程产生信号122 ()0t f t -<<⎧=⎨⎩其它，画出波形。 代码如下： clc clear f0=2; t=0:0.0001:2.5; y=square(w0*t,50); plot(t,y); axis([0 2.5 -1.5 1.5]) title('周期方波'); 图形如下： 单位阶跃信号 3. 分别画出2中()f t 移位3个单位的信号(3)f t -、反折后的信号()f t -、尺度变换后的信号(3)f t 。 代码如下： clc
clear t=-10:0.001:10; subplot(3,1,1) plot(t,f(t-3)) axis([-7 7 -2 2]) xlabel('t') ylabel('f(t-3)') title('移位') grid on subplot(3,1,2) plot(t,f(-t)) axis([-7 7 -2 2]) xlabel('t') ylabel('f(-t)') title('反折') grid on subplot(3,1,3) plot(t,f(3*t)) axis([-7 7 -2 2]) xlabel('t') ylabel('f(3t)') title('尺度变换') grid on 图形如下： x f (t ) x f (t -3) x f (-t ) x f (3*t ) 4. 用MATLAB 编程画出下图描述的函数。