随机过程上机实验报告讲解.pdf

合集下载

2021年随机过程实验报告

2021年随机过程实验报告

过程试验汇报班级: 通信1004班姓名: 杨靖学号: U13098试验目:了解数产生, 而且利用数来模拟均匀分布、 正态分布、 指数分布、 泊松分布而且计算均值和自相关序列。

试验工具:C++编程模拟试验原理:数产生原理: 经过数学算法产生伪数来, 模拟数产生。

数序列含有循环周期性。

能够证实, 任何产生伪数算法总会进入循环, 这么为了确保随机数序列不产生反复数据, 就要求循环周期足够长。

均匀分布产生原理:利用线性同余法(1)设置y0, 即设置种子(2)yn=kyn-1(mod N), un=yn/N泊松分布产生原理: 从泊松分布分布律可知, 采取前述方法很不适用。

因为: 所以, 采取递推法组成泊松分布: (1)产生均匀分布数u; (2) (3)若u<F, 令X=i, 停止; (4) (5)转向(3)。

正态分布产生原理:标准正态变量分布函数 反函数不存在显式, 所以也不能用逆变法产生。

故采取以下方法:设Ui ~U(0, 1), i=1,2,…,n, 且相互独立, 由中心极限定理可知, 当n 较大时设Ui ~U(0, 1), i=1,2,…,n, 且相互独立, E(Ui)=1/2, D(Ui)=1/12, 当n 较大时有:取n=12, 近似有:也就是说, 只要产生12个伪数u1,u2,…u12, 将它们加起来, 再减去6, 就能近似得到标准正态变量样本值。

{}!i i e p P X i i λλ-===11(1)!1i i i e p p i i λλλ+-+==++0,,;i p e F p λ-===/(1),,1;p p i F F p i i λ=+=+=+()~(0,1)n i i U nE U Z N -=∑~(0,1)ni n U Z N -=∑1216~(0,1)i i Z U N ==-∑指数分布产生原理:(1)产生均匀分布数{ui};(2)计算指数分布数: xi=-ln ui /λ试验代码:(1)数产生/*函数功效, 采取线性同余法, 依据输入种子数产生一个伪数, 假如种子不变,则将能够反复调用产生一个伪序列利用CMyRand类中定义全局变量: S, K, N, Y。

《随机过程》课程研究与实践

《随机过程》课程研究与实践

《随机过程》课程研究与实践本文将结合课堂教学和实践经验,以及相关理论知识,对《随机过程》课程进行研究和实践总结。

一、课程研究1.离散时间马尔可夫链:离散时间马尔可夫链是一种序列随机变量,其特点是过去的状态对下一个状态的概率分布没有影响,只与当前状态有关。

在研究中,我们要掌握马尔可夫链的转移矩阵、迭代计算等基本概念和方法。

2.连续时间马尔可夫链:与离散时间马尔可夫链相对应,连续时间马尔可夫链是一种连续随机变量,其特点是状态变量是连续的。

在研究中,我们要了解连续时间马尔可夫链的转移密度函数、平稳分布等基本概念和性质。

3.泊松过程:泊松过程是一种连续时间的独立增量随机过程,其特点是相邻两个事件的时间间隔服从指数分布。

在研究中,我们要掌握泊松过程的均值和方差计算、重要性质等基本理论和应用方法。

4.连续时间马尔可夫过程:连续时间马尔可夫过程是一种连续时间的马尔可夫链,其特点是状态变量是连续随机变量。

在研究中,我们要了解连续时间马尔可夫过程的转移密度函数、平稳分布等基本概念和性质。

二、课程实践1.实例分析:通过实践案例分析,将理论知识与实际问题相结合,加深学生对随机过程的理解和应用能力。

例如,通过分析金融领域中的股票价格变化,引导学生应用随机过程来建模并预测未来的价格变化。

2. 编程模拟:通过编程模拟的方式,让学生实际操作随机过程的理论知识,加深对概念和方法的理解。

例如,使用Python编写程序模拟随机过程并进行统计分析,进一步验证随机过程的性质和特点。

3.实验设计:设计实验,引导学生通过实验方法验证随机过程的基本特性。

例如,设计一组抛硬币实验,统计硬币正面朝上的次数,并利用极大似然估计法估计硬币正面朝上的概率。

4.论文研究:鼓励学生进行论文研究,深入探究随机过程在实际问题中的应用。

例如,学生可以选择在金融、通信或生物医学领域进行相关论文研究,分析和总结已有研究成果,并提出自己的见解。

三、总结与展望通过对《随机过程》课程的研究和实践,我们可以发现,这门课程的理论体系完备、应用广泛。

随机过程实验报告

随机过程实验报告

一、实验目的1. 理解随机过程的基本概念和性质。

2. 掌握随机过程的基本运算和性质。

3. 通过实验验证随机过程的性质和规律。

二、实验原理随机过程是指一系列随机变量按照一定规则排列而成的序列。

在现实生活中,随机过程广泛存在于自然界和人类社会,如股票价格、气象变化、生物进化等。

随机过程的研究有助于我们更好地理解和预测这些现象。

随机过程可以分为两类:离散随机过程和连续随机过程。

本实验主要研究离散随机过程。

三、实验设备与材料1. 计算机2. 随机过程模拟软件(如Matlab)3. 纸笔四、实验内容1. 随机过程的基本概念(1)随机变量的概念随机变量是指具有不确定性的变量,它可以取多个值。

在随机过程中,随机变量是基本的研究对象。

(2)随机过程的概念随机过程是由一系列随机变量按照一定规则排列而成的序列。

2. 随机过程的基本性质(1)无后效性无后效性是指随机过程的前后状态相互独立。

(2)无记忆性无记忆性是指随机过程的状态只与当前时刻有关,与过去时刻无关。

(3)马尔可夫性马尔可夫性是指随机过程的状态只与当前时刻有关,与过去时刻无关。

3. 随机过程的运算(1)随机过程的和设{Xn}和{Yn}是两个随机过程,则它们的和{Zn}定义为Zn = Xn + Yn。

(2)随机过程的差设{Xn}和{Yn}是两个随机过程,则它们的差{Zn}定义为Zn = Xn - Yn。

(3)随机过程的乘积设{Xn}和{Yn}是两个随机过程,则它们的乘积{Zn}定义为Zn = Xn Yn。

4. 随机过程的模拟利用随机过程模拟软件(如Matlab)模拟随机过程,观察其性质和规律。

五、实验步骤1. 初始化随机数生成器2. 定义随机过程(1)根据随机过程的基本性质,定义随机过程{Xn}。

(2)根据随机过程的运算,定义随机过程{Yn}。

3. 模拟随机过程(1)使用随机过程模拟软件(如Matlab)模拟随机过程{Xn}和{Yn}。

(2)观察模拟结果,分析随机过程的性质和规律。

随机过程上机实验报告讲解

随机过程上机实验报告讲解

2015-2016第一学期随机过程第二次上机实验报告实验目的:通过随机过程上机实验,熟悉Monte Carlo计算机随机模拟方法,熟悉Matlab的运行环境,了解随机模拟的原理,熟悉随机过程的编码规律即各种随机过程的实现方法,加深对随机过程的理解。

上机内容:(1 )模拟随机游走。

(2)模拟Brown运动的样本轨道。

(3)模拟Markov过程。

实验步骤:(1)给出随机游走的样本轨道模拟结果,并附带模拟程序。

①一维情形%—维简单随机游走% “从0开始,向前跳一步的概率为p,向后跳一步的概率为1-p”n=50;p=0.5;y=[0 cumsum(2.*(rand(1,n-1)v=p)-1)]; % n 步。

plot([0:n-1],y); %画出折线图如下。

w%一维随机步长的随机游动%选取任一零均值的分布为步长,比如,均匀分布。

n=50;x=rand(1,n)-1/2;y=[0 (cumsum(x)-l)];plot([0:n],y);②二维情形%在(u, v)坐标平面上画出点(u(k), v(k)), k=1:n,其中(u(k)) 和(v(k))是一维随机游动。

例%子程序是用四种不同颜色画了同一随机游动的四条轨道。

n=100000;colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(2,n)<0.5)-1;x=[zeros(1,2); cumsum(z')];col=colorstr(k);plot(x(:,1),x(:,2),col);③%三维随机游走 ranwalk3dp=0.5;n=10000; colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(3,n)v=p)-1; x=[zeros(1,3); cumsum(z')];col=colorstr(k);plot3(x(:,1),x(:,2),x(:,3),col);hold on end gridhold onendgrid4:04003?0-200-300-400-2OD20050、-100-200 -20D⑵给出一维,二维Brown运动和Poisson过程的模拟结果,并附带模拟程序,没有结果的也要把程序记录下来。

随机过程实验报告全

随机过程实验报告全

随机过程实验报告学院:专业:学号:姓名:一、实验目的通过随机过程的模拟实验,熟悉随机过程编码规律以及各种随机过程的实现方法,通过理论与实际相结合的方式,加深对随机过程的理解。

二、实验内容(1)熟悉Matlab工作环境,会计算Markov链的n步转移概率矩阵和Markov链的平稳分布。

(2)用Matlab产生服从各种常用分布的随机数,会调用matlab自带的一些常用分布的分布律或概率密度。

(3)模拟随机游走。

(4)模拟Brown运动的样本轨道的模拟。

(5)Markov过程的模拟。

三、实验原理及实验程序n步转移概率矩阵根据Matlab的矩阵运算原理编程,Pn = P ^n。

已知随机游动的转移概率矩阵为:P =0.5000 0.5000 00 0.5000 0.50000.5000 0 0.5000求三步转移概率矩阵p3及当初始分布为P{x0 = 1} = p{x0 = 2} = 0, P{x0 = 3} = 1 时经三步转移后处于状态3的概率。

代码及结果如下:P = [0.5 0.5 0; 0 0.5 0.5; 0.5 0 0.5] %一步转移概率矩阵P3 = P ^3 %三步转移概率矩阵P3_3 = P3(3,3) %三步转移后处于状态的概率1、两点分布x=0:1;y=binopdf(x,1,0.55);plot(x,y,'r*');title('两点分布');2、二项分布N=1000;p=0.3;k=0:N;pdf=binopdf(k,N,p);plot(k,pdf,'b*');title('二项分布');xlabel('k');ylabel('pdf');gridon;boxon3、泊松分布x=0:100;y=poisspdf(x,50);plot(x,y,'g.');title('泊松分布')4、几何分布x=0:100;y=geopdf(x,0.2);plot(x,y,'r*');title('几何分布');xlabel('x');ylabel('y');5、泊松过程仿真5.1 % simulate 10 timesclear;m=10; lamda=1; x=[];for i=1:ms=exprnd(lamda,'seed',1);x=[x,exprnd(lamda)];t1=cumsum(x);end[x',t1']5.2%输入:N=[];for t=0:0.1:(t1(m)+1)if t<t1(1)N=[N,0];elseif t<t1(2)N=[N,1];elseif t<t1(3)N=[N,2];elseif t<t1(4)N=[N,3];elseif t<t1(5)N=[N,4];elseif t<t1(6)N=[N,5];elseif t<t1(7)N=[N,6];elseif t<t1(8)N=[N,7];elseif t<t1(9)N=[N,8];elseif t<t1(10)N=[N,9];elseN=[N,10];endendplot(0:0.1:(t1(m)+1),N,'r-') 5.3% simulate 100 timesclear;m=100; lamda=1; x=[];for i=1:ms= rand('seed');x=[x,exprnd(lamda)];t1=cumsum(x);end[x',t1']N=[];for t=0:0.1:(t1(m)+1)if t<t1(1)N=[N,0];endfor i=1:(m-1)if t>=t1(i) & t<t1(i+1)N=[N,i];endendif t>t1(m)N=[N,m];endendplot(0:0.1:(t1(m)+1),N,'r-')6、泊松过程function I=possion(lambda,m,n)for j=1:mX=poissrnd(lambda,[1,n]); %参数为lambda的possion 过程N(1)=0;for i=2:nN(i)=N(i-1)+X(i-1);endt=1:n;plot(t,N)grid onhold onend7、布朗运动7.1一维布朗运动程序:function [t,w]=br1(t0,tf,h)t=t0:h:tf;t=t';x=randn(size(t));w(1)=0;for k=1:length(t)-1w(k+1)=w(k)+x(k);endw=sqrt(h)*w;w=w(:);end调用t0=1;tf=10;h=0.01;[t,w]=br1(t0,tf,h);figure;plot(t,w,'*');xlabel('t');ylabel('w');title('一维Brown运动模拟图'); 7.2二维布朗运动:function [x,y,m,n]=br2(x0,xf,y0,yf,h)x=x0:h:xf;y=y0:h:yf;a=randn(size(x));b=randn(size(y));m(1)=0;n(1)=0;for k=1:length(x)-1m(k+1)=m(k)+a(k);n(k+1)=n(k)+b(k);endm=sqrt(h)*m;n=sqrt(h)*n;end调用x0=0;xf=10;h=0.01;y0=0;yf=10;[x,y,m,n]=br2(x0,xf,y0,yf,h);figure;plot(m,n);xlabel('m');ylabel('n');title('二维Brown运动模拟图');7.3三维布朗运动:npoints =1000;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;8、马尔科夫链离散服务系统中的缓冲动力学m=200;p=0.2;N=zeros(1,m); %初始化缓冲区A=geornd(1-p,1,m); %生成到达序列模型, for n=2:mN(n)=N(n-1)+A(n)-(N(n-1)+A(n)>=1);endstairs((0:m-1),N);9、随机数游走9.1 100步随机游走n = 100; %选取步数。

随机过程实验

随机过程实验

齐次泊松过程的matlab数据分析一、参数设定:二、数据分析(一)通过分析统计Possion_data.txt求期望、方差。

1、数据导入:设置参数:lamda=5,仿真时间=10,样本函数数目=200;生成Possion_data.txt 文件,在matlab中使用“Import Data”功能,将txt文件所有行以“Matrix”格式导入Workplace空间,生成Possiondata变量。

如下图:图1. Possion_data原始数据图2. Possion_data数据导入图3. 生成Possiondata2、数据提取编程将Possiondata数组中的第3,6,9…300行提取出来形成一个新的数值poiss。

图4. 提取有效数据3、数据判别(1)、按照试验指导大纲,将时间间隔设成0.1,将每一条样本函数按照0.1的时间间隔进行统计,将在同一个0.1间隔内的数据归为一类。

得到“t1”图5. 数据判别归类(2)、采用“length”函数将上表中的数据进行计数得到以下参数:图6. 统计类中数量(3)、将“t2”进行累加,并同样的方法计算所有样本函数得到“s”:图7. 得到样本计数样本图8. 第一条样本计数过程图9. SJGC生成的第一条样本函数4、计算期望、方差使用“mean”函数计算期望值,使用“var”函数计算方差得到下图:图10. 均值_方差图从图中可以看出,泊松过程的均值与方差具有一致性。

图11. SJGC生成的均值函数图(二)求泊松过程的速率方法1根据所得到的的均值函数,使用“polyfit”函数采用一次函数模拟得到斜率4.9138,即为泊松速率。

方法2考虑到泊松事件的时间间隔是指数分布,且均值为泊松过程速率的倒数。

对样本函数进行处理,将两次到达时间相减得到每相邻两次事件发生的时间间隔,使用“expfit”函数得到估计的均值,对其求倒得到泊松速率。

使用循环语句得到每一条样本函数的速率,最后求平均得到要求的泊松速率。

随机过程-实验报告

随机过程-实验报告
1 / 2 P 1/3 1/ 3 1/3 1/3 1/ 2 1/ 6 1/3 1/ 6
(1) 计算 2 步转移概率;(2) 已知初始分布为 P 2 / 5, 2 / 5,1 / 5 ,求 X 2 的分布律 (3) 求平稳分布,要求给出程序与结果。 程序:
程序: p=[0.2 0.8 0;0.8 0 0.2; 0.1 0.3 0.6]; P2=p^2 a=[p'-eye(3);ones(1,3)];b=[0 0 0 1]';T=a\b 结果:
0.1389
0.0611
解:由题意可知,该问题的转移概率矩阵 P 为:
8
0 .2 P 0 .8 0 .1
9
实验三
实验题目 实验目的 实验地点及时间 模拟 Possion 流 用 Matlab 语言产生随机数,了解 Possion 流 信息楼 127 机房 2012 年 6 月 4 日
4
实验内容 用 Matlab 语言产生随机数,并编程实现 possion 流的模拟 程序: U=rand(1,20); a=2; X=-a^(-1)*log(U); S=zeros(1,22); d=zeros(1,22); S(1)=0;S(2)=X(1); for n=3:21 S(n)=S(n-1)+X(n-1); end for i=0:21 %--if 0<=i<S(2) d(i+1)=0; else for j=2:21 if (S(j)<=i)&(S(j+1)<i) d(i+1)=j; end end end end plot(d)
实验内容 判定一个 Markov 链是否是遍历的,若是遍历的,求其极限分布。并能从实际问 题中抽象出 Markov 链,并求出其极限分布,并理解其实际意义。 实验习题 1、已知齐次马氏链 X n , n 0,1, 2, 的状态空间 E 1, 2, 3 ,状态转移矩阵为

随机过程上机实验报告-华中科技大学--HUST

随机过程上机实验报告-华中科技大学--HUST
流程图:
6、/*函数功能,计算任意分布的随机过程的均值
*/
double CMyRand::Ex(void)
{
double Ex = 0;
//添加均值计算代码
int i;
double sum=0;
for(i=0;i<500;i++)
sum+=AverageRandom(0,2);
Ex=sum/i;
return Ex;
Miu为均值,sigma为标准差
*/
double CMyRand::NormalRandom(double miu, double sigma, double min, double max)
{
double dResult;
dResult = 0;
int i,n;
double sum=0.0;
n=200;
我们在示波器界面上点击一个按钮它就会执行这个按钮所对应功能比如点击正态分布它就会调用crandomdlg中的对应凼数在调用cmyrand中的产生正态分布的凼数再将结果送到cscope类中迚行显示最后我们可以在示波器上看到图形

班级:通信1301班
姓名:郭世康
学号:U201313639
指导教师:卢正新
dResult=dResult*(max-min)+min;//将0~1之间的均匀分布通过乘以倍数放大到到min~max
return dResult;
}
输入参数为min,max,即均匀分布的范围。输出参数为dResult,即为随机序列。
流程图:
3、
/*函数功能,根据大数定律,在min到max范围内产生正态分布的随机数
流程图:

哈尔滨工程大学2017年随机过程上机作业

哈尔滨工程大学2017年随机过程上机作业

物业化管理调研报告8篇(经典版)编制人:__________________审核人:__________________审批人:__________________编制单位:__________________编制时间:____年____月____日序言下载提示:该文档是本店铺精心编制而成的,希望大家下载后,能够帮助大家解决实际问题。

文档下载后可定制修改,请根据实际需要进行调整和使用,谢谢!并且,本店铺为大家提供各种类型的经典范文,如调研报告、总结报告、述职报告、心得体会、自我鉴定、条据文书、合同协议、教学资料、作文大全、其他范文等等,想了解不同范文格式和写法,敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you!Moreover, our store provides various types of classic sample essays for everyone, such as research reports, summary reports, job reports, reflections, self-evaluation, normative documents, contract agreements, teaching materials, essay summaries, and other sample essays. If you want to learn about different sample essay formats and writing methods, please stay tuned!物业化管理调研报告8篇认真撰写调研报告是对调研任务的最终总结,也是对相关领域的专业呈现,调研报告是一种详尽记录研究、分析和发现的文档,旨在提供全面的信息,以下是本店铺精心为您推荐的物业化管理调研报告8篇,供大家参考。

随机过程实验报告

随机过程实验报告

随机过程实验报告一、实验问题两赌徒模型对于上述模型现在假定赌徒甲的对手赌徒乙有N-i的初始财富,N为两个赌徒的总财富。

则赌徒甲破产的概率有多大?模拟之。

二、问题分析该问题实质上为带有两个吸壁的随机游动,我们可以仍可把它看作数学中的一个一维随机游动问题。

其马尔可夫链状态空间为{0,1,2,…,N},N为赌徒甲、乙的总财富。

类似于赌徒与游戏机模型,我们也可以把财富抽象地看成是一个质点。

可知求赌徒甲破产的概率转化为现在的问题就是求质点从i点出发到达0状态先于到达N状态的概率。

这里较赌徒与游戏机模型中多出一个条件,即:赌徒甲先于赌徒乙到达0状态。

我们不难得到这一模型的解:三、问题解决1、先讨论p=q的随机游动情况对于简单的随机游动,如果从0开始,向前跳一步的概率为p,向后跳一步的概率为1-p,则由计算机可以模拟此情形。

这只是许多模拟结果中的一种。

现在我们假设,有A、B两个赌徒,他们共同用于赌博的财富M=100(元),A、B输赢的概率(即赌博的技巧相同)时,他们破产的概率。

假设,共同的财富中A、B分别投入的资金如下表:运算结果如下:由上图可知,当赌徒甲、乙输赢的概率相等时,其中一人破产的概率与对方所拥有的财富成正比关系。

这样我们可以得出结论:在两人的赌博游戏中,如果赌徒甲、乙的赌博技术差不多即输赢概率相当的话,那么谁要想最终获胜的最好方法就是多带赌本。

2、下面讨论p!=q时随机游动情况我们不妨将之具体为p=0.4,q=0.6。

用计算机模拟上述数据。

可得图如下:由上图可知,在每次输赢都为1元时,就算甲90元、乙10元,甲也几乎不可能赢。

如果我们把每次下的赌注加大到5元,修改程序三,模拟之,又可得图如下:由上图我们可以更清晰地看出:在两人的赌博游戏中,如果赌徒甲的赌博技术比乙的赌博技术差的话,那么甲要想最终获胜就要带比乙多很多的赌本。

四、结果拓展现实中的赌博还可能有三人、四人甚至更多的人一起进行。

下面我们简单地讨论当赌徒输赢概率相等时的二维随机游动。

随机过程实验报告

随机过程实验报告

师范学院
实验报告
2)计算x的方差(var())
实 验 内 容 及 步 骤 (续)
3) 利用rnorm 生成100个随机数赋值给变量y ,计算x 、y 的协方差(cov()),相关系数(cor())。

4) 将变量x 、y 放到数据框A 中(用data.frame()函数),计算A 的协方差矩阵(cov())、相关系数矩阵(cor()).
2.利用rnorm 函数生成随机向量x (长度100),将x 转换为一个10阶矩阵A 。

rnorm
1) 计算A 的特征值和特征向量
2)计算A的逆
3)随机生成长度为10的向量b,求解线性方程组Ay b
3.随机生成随机向量x(长度10)。

将x的数据复制到文本文件data.txt中,并保存。

利用read.table()函数,将data.txt中的数据重新读入到R工作空间中,并命名为y。

实验心得
这次试验是我第一次接触R语言,刚开始遇到了很多困难,对于R语言一窍不通,后来经过老师的悉心指导,以及自己积极的去查找资料,对R 语言有了进一步的了解。

这次试验通过随机生成数字,并且对对其进行运、组合、求解,定义函数,以及生成表。

在此过程中,锻炼了自己的学习能力、动手操作能力,也让我对R语言产生了兴趣,发现R语言与matlab 相似之处,以及方便之处。

希望以后有机会可以更加系统的掌握、了解R语言,并达到熟练的应用。

学习并提升了R语言的使用能力,更好的实现课本与实际操作的结合,让我更好的学习并理解了应用随机过程这门课。

评语。

随机过程作业报告

随机过程作业报告

随机过程作业报告姓名:学号:班级:一、问题定义(1)问题描述:信号源发信号,只能发1和10,发送1的概率是0.7,发送10的概率是0.3。

信道上加一个高斯噪声,满足均值为0,方差为5的正态分布。

利用计算机仿真,得到正确检验概率。

(2)实验要求:·利用C++编程实现对信源发送信号并配合噪声干扰、门限计算和仿真检测的流程编译。

·根据MAP规划,做10000次仿真检测。

(3)输入/输出:无输入,输出成功检验概率。

(4)程序功能:对信号源随机发送的信号进行检验,并得到正确检验概率。

二、概要设计(1)解题思路:通过使用TR1扩展在c++中用随机数生成设施生成随机数(具体使用定义在头文件random中的随机数库中的随机数引擎类和随机数分布类来生成范围在(0,1)服从均匀分布的随机数和均值为0,方差为5的服从高斯分布的随机数)。

根据生成的均匀分布随机数来决定输入信号(范围在(0,0.7)为1,范围在(0.7,1)为10)。

计算门限值,将信号加上噪声后的值与门限相比,判断接收到的信号与发送信号是否相同,循环10000次,累积发送信号与接受信号相同的次数,得到检验结果。

(2)头文件:<time.h> /用于设置种子时,使用time函数<random> /用于使用random头文件中的随机引擎类和随机分布类<iostream> /输入输出操作是通过该头文件中的输入输出流cin和cout来实现的。

(3)命名空间:std::tr1(4)随机引擎类:Mersenne -Twister generator(5)随机分布类:uniform_real类生成服从均匀分布的浮点数normal_distribution类生成服从高斯分布的浮点数(6)主程序流程:定义命名空间中的引擎类和分布类、设置种子—>使用for循环,循环次数10000次—>在循环中调用定义过的随机数生成函数,生成随机数—>根据生成随机数范围确定发送信号—>将加上噪声后的信号与计算得到的门限值比较,判断检验是否正确—>输出成功检验概率三、详细设计(1)源程序:#include<time.h>#include<iostream>#include<random>using namespace std::tr1;int main(){std::tr1::mt19937 e; // 随机引擎类:Mersenne Twister generator std::tr1::uniform_real<double>unif(0, 1); //定义均匀分布类,范围为0-1 std::tr1::normal_distribution<double>normal(0, sqrt(5.0));//定义正态分布类,均值为0,方差为5//e.seed((unsigned int)time(NULL));//设置种子,调用time函数float counter = 0.0;for (int i = 0; i<10000; i++){double s = unif(e);double x, y;if (s>0.7) x = 10;else x = 1;y = x + normal(e); //信号和噪声叠加std::cout << y << std::endl;if (y > 5.9707) y = 10;else y = 1;if (y == x) counter++;}counter = counter / 10000;std::cout <<"正确检测概率="<< counter << std::endl;system("pause");return 0;}(2)结果:。

随机过程实验讲义

随机过程实验讲义

随机过程实验讲义刘继成华中科技大学数学与统计学院2011-2012年上半年为华中科技大学数学系本科生讲授随机过程课程参考资料前言 (1)第一章Matlab 简介 (2)第二章简单分布的模拟 (6)第三章基本随机过程 (9)第四章Markov过程 (12)第五章模拟的应用和例子 (16)附录各章的原程序 (51)参考文献 (75)若想检验数学模型是否反映客观现实,最自然的方法是比较由模型计算的理论概率和由客观试验得到的经验频率。

不幸的是,这两件事都往往是费时的、昂贵的、困难的,甚至是不可能的。

此时,计算机模拟在这两方面都可以派上用场:提供理论概率的数值估计与接近现实试验的模拟。

模拟的第一步自然是在计算机程序的算法中如何产生随机性。

程序语言,甚至计算器,都提供了“随机”生成[0,1]区间内连续数的方法。

因为每次运行程序常常生成相同的“随机数”,因此这些数被称为伪随机数。

尽管如此,对于多数的具体问题这样的随机数已经够用。

我们将假定计算机已经能够生成[0,1]上的均匀随机数。

也假定这些数是独立同分布的,尽管它们常常是周期的、相关的、……。

……本讲义的安排如下,第一章是Matlab简介,从实践动手角度了解并熟悉Matlab环境、命令、帮助等,这将方便于Matlab的初学者。

第二章是简单随机变量的模拟,只给出了常用的Matlab 模拟语句,没有堆砌同一种变量的多种模拟方法。

对于没有列举的随机变量的模拟,以及有特殊需求的读者应该由这些方法得到启发,或者参考更详细的其他文献资料。

第三章是基本随机过程的模拟。

主要是简单独立增量过程的模拟,多维的推广是直接的。

第四章是Markov过程的模拟。

包括服务系统,生灭过程、简单分支过程等。

第五章是这些模拟的应用。

例如,计算概率、估计积分、模拟现实、误差估计,以及减小方差技术,特别给读者提供了一些经典问题的模拟,通过这些问题的模拟将会更加牢固地掌握实际模拟的步骤。

平稳过程的模拟、以及利用平稳过程来预测的内容并没有包含在本讲义之内,但这丝毫不影响该内容的重要性,这也是将会增补进来的主要内容之一。

随机过程实验报告

随机过程实验报告

随机过程实验报告随机过程实验报告一、引言随机过程是概率论和数理统计中的一个重要分支,它研究的是随机事件随时间的演化规律。

在现实生活中,我们经常会遇到各种各样的随机过程,比如天气变化、股票价格波动、人口增长等等。

本次实验旨在通过实际观测和数据分析,探究随机过程的特性和规律。

二、实验目的本次实验的主要目的是研究和分析一个具体的随机过程,以加深对随机过程理论的理解。

通过实际观测和数据分析,我们将探究该随机过程的概率分布、平均值、方差等统计特性,并尝试利用数学模型对其进行建模和预测。

三、实验方法我们选择了一个经典的随机过程作为研究对象:骰子的投掷。

我们将进行多次骰子投掷实验,并记录每次投掷的结果。

通过统计分析这些结果,我们可以得到骰子的概率分布、平均值和方差等重要参数。

四、实验过程我们使用了一颗标准的六面骰子进行了100次投掷实验。

每次投掷后,我们记录了骰子的点数,并将这些数据整理成了一个数据集。

五、实验结果通过对实验数据的统计分析,我们得到了以下结果:1. 概率分布我们统计了每个点数出现的次数,并计算了它们的频率。

结果显示,每个点数的频率接近于1/6,符合骰子的均匀分布特性。

2. 平均值我们计算了所有投掷结果的平均值,发现它接近于3.5。

这是因为骰子的点数从1到6,平均为(1+2+3+4+5+6)/6=3.5。

3. 方差我们计算了所有投掷结果的方差,发现它接近于2.92。

方差是衡量随机变量离其均值的分散程度的指标,它的大小反映了骰子点数的变化范围。

六、讨论与分析通过对实验结果的分析,我们可以得出以下结论:1. 骰子的点数具有均匀分布的特性,每个点数出现的概率接近于1/6。

2. 骰子的平均值为3.5,这是由于骰子的点数从1到6,平均为(1+2+3+4+5+6)/6=3.5。

3. 骰子的方差为2.92,这意味着骰子的点数变化范围较大。

通过以上结果,我们可以看出骰子的投掷过程是一个典型的随机过程。

它符合随机过程的基本特性,即随机性和不可预测性。

概率论与随机过程上机实验报告

概率论与随机过程上机实验报告

概率论与随机过程上机实验报告题目一题目对二项分布事件的概率的精确计算与用泊松分布和中心极限定理的近似计算进行对比。

P变化n固定,进行比较n固定,p变化进行比较。

源代码运行结果黑星代表二项分布,蓝色是泊松分布绿线是中心极限定理小结n变化从50开始到150,中心极限定理的计算方法更加接近二项分布的精确计算,泊松分布于精确计算差距稍微增大但保持原有的变化趋势。

p改变时,p=0.5时取最大值,仍然是中心极限定理比泊松分布更加接近二项分布精确计算。

第二题题目对正态总体参数的区间估计,进行验证及区间长度的变化情况(注:对一个参数,验证一种情形即可)。

(a)样本容量固定,置信度变化;(b)置信度固定,样本容量变化。

源程序运行结果小结可以看出来,当样本容量不断增加时,区间估计的精度越来越高;同时,当置信度不断提高时,区间估计的精度也越来越高。

第三题题目自己选一个总体,验证样本k阶矩的观察值随样本容量的增大与总体k阶矩接近程度(对k=1,2进行验证)源代码运行结果小结使用自由度为10的卡方分布作为研究总体,取样本容量大小从1到10000。

图像表明,,随着样本容量的增加,样本观测值的一阶原点矩和二阶原点矩都越来越接近于总体的一二阶原点矩,即10和120。

第五题题目自己设计一种情形,当样本至少为多少时,产品的合格率才能符合给定的合格率源程序运行结果小结观察可知,卡方分布产生的500个随机数的统计直方图的形状与真实卡方分布曲线形状基本拟合。

个人感想之前大一在进行数学建模的时候通常要用到数理统计的相关知识,但由于没有系统的学习过,始终是一知半解。

经过一学期对概率论与随机过程的学习,掌握了很多统计学上的观点以及方法,这对之后的工作或是科研都有着很大的作用。

经过这次的上机实验,也能让我们从编程的角度更深入的理解一些方法在实践中的用法,受益匪浅。

最后,感谢老师一学期的辛勤教学,也希望老师之后身体健康工作顺利。

随机过程实验报告

随机过程实验报告

随机过程试验报告班级:信息与计算科学2010级1班姓名:李翠珍学号:20104609实验实验总结:本次试验熟练的掌握了三维图像的matlab 编程语句,最重要的是学习了 rnd 使x 为泊松随机数。

4实验二2o-1-228均值函数已知u=0,令自变量x的取值范围[-1,1]x=-1:0.01:1u=0;plot(x,u,'-+');方差函数var(x(t)) =2;在matlab 中用v 代替方差,同样令x=-1:0.01:1v=0;plot(x,v,'-+');-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 10.80.60.40.2-0.2-0.4-0.6-0.8-1 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8自相关函数令t 1 , t 2的范围为[0,2二].根据已知条件编写下程序: t1=0:0.01:2;t2=0:0.01:2; t=t1-t2;w=3;r=2*cos(w*t); plot3(r,t1,t2); axis square; grid on;实验总结:本次试验,主要是借鉴了课本 2.2的课上例题,利用了随机过程中的 中的相关公式求解。

3实验三实验总结:本试验主要锻炼了我们从大量信息内摘取有用信息的能力学习及理解 运用新知识。

我在这方面比较欠缺,以后一定要多加练习。

3030实验成绩 评阅时间 评阅教师 2520151055 10 15 20 25实验内容判定一个Markov 链是否是遍历的,若是遍历的,求其极限分布。

并能从实际问 题中抽象出Markov 链,并求出其极限分布,并理解其实际意义。

实验习题课本p125 5.7 将两个红球,四个白球分别放入甲乙两个盒子中。

每次从两个 盒子中各取一球交换,以X n 记第n 次交换后甲盒中红球数。

(1) 说明{ X n , n=0,1,…}是一 Markov 链并求转移矩阵P ;(2) 试证{ X n , n=0 , 1,…}是遍历的;(3) 求它的极限分布;(1) 设X n 记第n 次交换后甲盒中红球数,则易见{ X n ,n=0,1,…}是状态空间S 二{0,1,2}的Markov 链,一步转移概率矩阵为:123 8 0(2) 由于状态空间S 有限,且状态互通,故{ X n , n=0,1,-}不可约,从而 正常返,又状态1为非周期的,故{ X n , n=0 , 1,…}还是遍历链。

哈尔滨工程大学 研究生 随机过程课程实验报告~.

哈尔滨工程大学 研究生 随机过程课程实验报告~.

第一题:用PC机产生[0,1]均匀分布的白色序列{}kkX(=,3,2,1),2000(1) 打印出前50个数{}iiX=),,(3,2,150(2) 分布检验(3) 均值检验(4) 方差检验(5) 计算出相关函数{}±±=iB(±i),10,2,,0,1x源程序:clear;clc;x=rand(1,2000);fprintf('1.输出前50个数:');for i=1:5j=1:10;X(i,j)=x((i-1)*10+j);endX % 打印出前50个数y1=x(find(x>=0&x<0.1));t(1)=length(y1);y2=x(find(x>=0.1&x<0.2));t(2)=length(y2);y3=x(find(x>=0.2&x<0.3));t(3)=length(y3);y4=x(find(x>=0.3&x<0.4));t(4)=length(y4);y5=x(find(x>=0.4&x<0.5));t(5)=length(y5);y6=x(find(x>=0.5&x<0.6));t(6)=length(y6);y7=x(find(x>=0.6&x<0.7));t(7)=length(y7);y8=x(find(x>=0.7&x<0.8));t(8)=length(y8);y9=x(find(x>=0.8&x<0.9));t(9)=length(y9);y10=x(find(x>=0.9&x<1));t(10)=length(y10) ;fprintf('2.分布检验:');tsubplot(2,1,1);hist(x,10); % 分布检验fprintf('3.均值检验:');EX=mean(x) % 均值检验fprintf('4.方差检验:');DX=var(x) % 方差检验fprintf('5.计算相关函数:');for m=-10:1:10j=2000-abs(m);for i=1:jC(i)=(x(abs(m)+i)-EX).*(x(i)-EX);endB(m+11)=sum(C)/j;endfor i=1:3j=1:7;Bx(i,j)=B((i-1)*7+j);endBx % 计算相关函数subplot(2,1,2)m=-10:10;plot(m,B)1.输出前50个数:X =Columns 1 through 80.1315 0.6175 0.4759 0.0236 0.8753 0.0960 0.5479 0.07460.8483 0.4888 0.4260 0.5609 0.6730 0.1103 0.7614 0.49120.5077 0.5892 0.0702 0.0386 0.4879 0.3002 0.0358 0.79340.4440 0.4423 0.5000 0.0325 0.0196 0.2932 0.0558 0.72080.8507 0.1279 0.4534 0.6225 0.4175 0.6702 0.0820 0.8725 Columns 9 through 100.9542 0.25160.5314 0.59830.5083 0.01650.8429 0.54420.4153 0.55662.分布检验:t =210 192 197 202 197 214 198 191 188 211图(1)分布检验3.均值检验:理论值:EX =0.5实际值:EX =0.49914.方差检验:理论值:DX =1/12实际值:DX =0.0839均值和方差表:5.计算相关函数:Bx =0.0022 0.0011 -0.0010 -0.0014 -0.0013 0.0034 -0.0051-0.0026 0.0018 -0.0019 0.0838 -0.0018 0.0019 -0.0025 -0.0051 0.0033 -0.0014 -0.0015 -0.0013 0.0009 0.0020图(2)相关函数第二题:用PC机产生()1,0kkX(=),N分布的正态序列{},20003,2,1(1)打印出前50个数{}ii=X3,2,1,(50),(2)分布检验(3)均值检验(4)方差检验(5)计算出相关函数{}±±=iB(±i),10,2,,0,1x源程序:clear;clc;x=randn(1,2000);fprintf('1.输出前50个数:');for i=1:5j=1:10;X(i,j)=x((i-1)*10+j);endX % 打印出前50个数y1=x(find(x>=0&x<0.1));t(1)=length(y1);y2=x(find(x>=0.1&x<0.2));t(2)=length(y2);y3=x(find(x>=0.2&x<0.3));t(3)=length(y3);y4=x(find(x>=0.3&x<0.4));t(4)=length(y4);y5=x(find(x>=0.4&x<0.5));t(5)=length(y5);y6=x(find(x>=0.5&x<0.6));t(6)=length(y6);y7=x(find(x>=0.6&x<0.7));t(7)=length(y7);y8=x(find(x>=0.7&x<0.8));t(8)=length(y8);y9=x(find(x>=0.8&x<0.9));t(9)=length(y9);y10=x(find(x>=0.9&x<1));t(10)=length(y10) ;fprintf('2.分布检验:');tsubplot(2,1,1);hist(x,10); % 分布检验fprintf('3.均值检验:');EX=mean(x) % 均值检验fprintf('4.方差检验:');DX=var(x) % 方差检验fprintf('5.计算相关函数:');for m=-10:1:10j=2000-abs(m);for i=1:jC(i)=(x(abs(m)+i)-EX).*(x(i)-EX);endB(m+11)=sum(C)/j;endfor i=1:3j=1:7;Bx(i,j)=B((i-1)*7+j);endBx % 计算相关函数subplot(2,1,2)m=-10:10;plot(m,B)1.输出前50个数:X =Columns 1 through 8-1.0457 -1.0045 -0.7384 -0.9445 -0.1354 -0.4226 1.5979 -0.38110.3409 0.5486 -1.0160 -1.6335 -1.8104 -0.0349 0.6758 -0.8909-0.9381 -1.5436 0.1596 -0.3688 -1.0122 0.1134 0.8850 -0.5823 -0.3197 1.6065 1.0613 0.3005 0.3511 0.9522 -0.6329 -0.8587 -0.0243 0.9170 -0.5015 -0.2513 1.6728 -1.3644 -0.3351 1.2946 Columns 9 through 100.2348 -0.3093-1.8913 2.2175-0.7176 -0.67331.7461 -0.55610.4811 -0.25202.分布检验:t =71 81 70 78 66 62 71 71 58 49图(3)分布检验3.均值检验: 理论值:EX =0实际值:EX = -0.0054 4.方差检验: 理论值:DX =1实际值:DX = 0.9916 均值和方差表:5.计算相关函数: Bx =-0.0097 -0.0258 -0.0077 0.0131 0.0244 -0.0224 0.0590 0.0228 0.0272 0.0208 0.9911 0.0212 0.0271 0.0224 0.0588 -0.0223 0.0238 0.0125 -0.0083 -0.0255 -0.0082图(4)相关函数第三题:设{}1000,3,2,1),( =k k ε为正态白色序列,服从()1,0N 分布,()()()14-+=k k k X εε,1000,3,2,1 =k 求(1) ()()∑==1000110001k k X k EX (2) ()()∑==100012210001k k X k EX(3) ()()()[]22k EX k EX k DX -=(4) ()()[]()[]{}∑-=--+=mn xxx m n X m m n X m B 1000110001,10,,2,1,0±±±= m ,并画出x B (m)图源程序:clfclearp=randn(1,1001);k=2:1001;x=p(k)+4.*p(k-1);m=mean(x)m1=mean(x.^2)s=m1-m.^2for i=-10:10l=0;p=1000-abs(i);for k=1:pl=l+[x(k+abs(i))-m]*[x(k)-m];endb(i+11)=l/p;endi=-10:10;plot(i,b)1. 均值EX:理论值:EX =0实际值:EX =-0.02092.均方值:EX^2:理论值:EX^2= 17实际值:EX^2= 16.39993. 方差DX:理论值:DX = 17实际值:DX = 16.3995均值和方差表:4. 相关函数:B(m) =x0.1462 0.6689 -0.0319 0.1473 0.1085 -0.3709 0.2299 0.6721 0.1214 4.4418 16.8904 4.4418 0.1214 0.6721 0.2299 -0.3709 0.1085 0.1473 -0.0319 0.6689 0.1462图(5) 相关函数第四题:设{()k ξ,k=0,1,2,…}为N (0,1)正态白序列,()k ξ~N(0,1) 令()()()0.7071X k X k k ξ+-=,k =1,2,…,1000; ()10X -=。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

2015-2016第一学期随机过程第二次上机实验报告实验目的:通过随机过程上机实验,熟悉Monte Carlo计算机随机模拟方法,熟悉Matlab的运行环境,了解随机模拟的原理,熟悉随机过程的编码规律即各种随机过程的实现方法,加深对随机过程的理解。

上机内容:(1)模拟随机游走。

(2)模拟Brown运动的样本轨道。

(3)模拟Markov过程。

实验步骤:(1)给出随机游走的样本轨道模拟结果,并附带模拟程序。

①一维情形%一维简单随机游走%“从0开始,向前跳一步的概率为p,向后跳一步的概率为1-p”n=50;p=0.5;y=[0 cumsum(2.*(rand(1,n-1)<=p)-1)]; % n步。

plot([0:n-1],y); %画出折线图如下。

%一维随机步长的随机游动%选取任一零均值的分布为步长, 比如,均匀分布。

n=50;x=rand(1,n)-1/2;y=[0 (cumsum(x)-1)];plot([0:n],y);②二维情形%在(u, v)坐标平面上画出点(u(k), v(k)), k=1:n, 其中(u(k))和(v(k)) 是一维随机游动。

例%子程序是用四种不同颜色画了同一随机游动的四条轨道。

n=100000;colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(2,n)<0.5)-1;x=[zeros(1,2); cumsum(z')];col=colorstr(k);plot(x(:,1),x(:,2),col);hold onendgrid③%三维随机游走ranwalk3d p=0.5;n=10000;colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(3,n)<=p)-1;x=[zeros(1,3); cumsum(z')]; col=colorstr(k);plot3(x(:,1),x(:,2),x(:,3),col);hold onendgrid(2) 给出一维,二维Brown运动和Poisson过程的模拟结果,并附带模拟程序,没有结果的也要把程序记录下来。

①一维Brown% 这是连续情形的对称随机游动,每个增量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);②二维Brownnpoints = 5000;dt = 1;bm = cumsum([zeros(1, 3); dt^0.5*randn(npoints-1, 3)]); figure(1);plot(bm(:, 1), bm(:, 2), 'k');pcol = (bm-repmat(min(bm), npoints, 1))./ ...repmat(max(bm)-min(bm), npoints, 1);hold on;scatter(bm(:, 1), bm(:, 2), ...10, pcol, 'filled');grid on;hold off;③三维Brownnpoints = 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;④%泊松过程的模拟、检验及参数估计syms Un X S;n=10;%生成n*n个随机数r=1;%参数temp=0;tem=0;Un=rand(n,1);%共产生n*n个随机数for i=1:1:nX(i)=-log(Un(i))/r;endX=subs(X);for i=1:1:nfor j=1:1:itemp=temp+X(j);endS(i)=temp;temp=0;endS=subs(S);%检验泊松过程使用第四条for i=1:1:ntem=tem+S(i);endsigmaN=tem;T=S(n);alpha=0.05;%置信水平p=sigmaN/T;p1=(1/2)*(n-1.96*(n/3)^(1/2)); p2=(1/2)*(n+1.96*(n/3)^(1/2)); c1=subs(p-p1)c2=subs(p-p2)if (c1<=0&c2>=0)|(c1>=0&c2<=0)disp('这是一个泊松过程!')%参数估计使用极大似然估计r_=subs(n/T);if abs(subs(r_-r))<0.1disp('参数估计正确!')disp('参数估计值为:')r_end%绘制轨迹y=0;x=0:0.001:subs(S(1));plot(x,y)for k=1:1:ny=k;x=subs(S(k)):0.001:subs(S(k+1));hold onplot(x,y)endelsedisp('这不是一个泊松过程!') end⑤%二维poisson2d lambda=100;nmb=poissrnd(lambda)x=rand(1,nmb);y=rand(1,nmb);gridscatter(x,y,5,5.*rand(1,nmb));⑥%三维poisson3d%单位体积的泊松点数强度为lambda lambda=100;nmb=poissrnd(lambda)x=rand(1,nmb);y=rand(1,nmb);z=rand(1,nmb);gridscatter3(x,y,z,5,5.*rand(1,nmb));(3)Markov过程的模拟结果。

①离散服务系统中的缓冲动力学%离散服务系统中的缓冲动力学m=200;p=0.2;N=zeros(1,m); %初始化缓冲区A=geornd(1-p,1,m); %生成到达序列模型, 比如,几何分布for n=2:mN(n)=N(n-1)+A(n)-(N(n-1)+A(n)>=1);endstairs((0:m-1),N);②M/M/1模型% [tjump, systsize] = simmm1(n, lambda, mu) % Inputs: n - number of jumps% lambda - arrival intensity% mu - intensity of the service times % Outputs: tjump - cumulative jump times% systsize - system size% set default parameter values if ommited Ⅰ:nargin=0nargin=0;if (nargin==0)n=500;lambda=0.8;mu=1;endi=0; %initial value, start on level itjump(1)=0; %start at time 0systsize(1)=i; %at time 0: level ifor k=2:nif i==0mutemp=0;elsemutemp=mu;endtime=-log(rand)/(lambda+mutemp); % Inter-step times:%Exp(lambda+mu)-distributedif rand<lambda/(lambda+mutemp)i=i+1; %jump up: a customer arriveselsei=i-1; %jump down: a customer is departing end %ifsystsize(k)=i; %system size at time itjump(k)=time;end %for itjump=cumsum(tjump); %cumulative jump times stairs(tjump,systsize);Ⅱ:nargin不为0时nargin=2;if (nargin==0)n=500;lambda=0.8;mu=1;endi=0; %initial value, start on level itjump(1)=0; %start at time 0systsize(1)=i; %at time 0: level ifor k=2:nif i==0mutemp=0;elsemutemp=mu;endtime=-log(rand)/(lambda+mutemp); % Inter-step times:%Exp(lambda+mu)-distributedif rand<lambda/(lambda+mutemp)i=i+1; %jump up: a customer arriveselsei=i-1; %jump down: a customer is departingend %ifsystsize(k)=i; %system size at time itjump(k)=time;end %for itjump=cumsum(tjump); %cumulative jump timesstairs(tjump,systsize);③M/D/1系统% function [jumptimes, systsize] = simmd1(tmax, lambda)% SIMMD1 simulate a M/D/1 queueing system. Poisson arrivals% of intensity lambda, deterministic service times S=1.% [jumptimes, systsize] = simmd1(tmax, lambda)% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - time points of arrivals or departures % systsize - system size in M/D/1 queue% systtime - system times% Authors:% v1.2 07-Oct-02% set default parameter values if ommitedⅠ:nargin=0nargin=0;if (nargin==0)tmax=1500;lambda=0.95;endarrtime=-log(rand)/lambda; % Poisson arrivalsi=1;while (min(arrtime(i,:))<=tmax)arrtime = [arrtime; arrtime(i, :)-log(rand)/lambda];i=i+1;endn=length(arrtime); % arrival times t_1,...t_narrsubtr=arrtime-(0:n-1)'; % t_k-(k-1)arrmatrix=arrsubtr*ones(1,n);deptime=(1:n)+max(triu(arrmatrix)); % departure times%u_k=k+max(t_1,..,t_k-k+1)B=[ones(n,1) arrtime ; -ones(n,1) deptime'];Bsort=sortrows(B,2); % sort jumps in order jumps=Bsort(:,1);jumptimes=[0;Bsort(:,2)];systsize=[0;cumsum(jumps)]; % M/D/1 processsysttime=deptime-arrtime'; % system times% plot a histogram of system timeshist(systtime,30);Ⅱ:nargin不为0% function [jumptimes, systsize] = simmd1(tmax, lambda)% SIMMD1 simulate a M/D/1 queueing system. Poisson arrivals% of intensity lambda, deterministic service times S=1.%% [jumptimes, systsize] = simmd1(tmax, lambda)%% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - time points of arrivals or departures% systsize - system size in M/D/1 queue% systtime - system times% Authors:% v1.2 07-Oct-02% set default parameter values if ommitednargin=2;if (nargin==0)tmax=1500;lambda=0.95;endarrtime=-log(rand)/lambda; % Poisson arrivalsi=1;while (min(arrtime(i,:))<=tmax)arrtime = [arrtime; arrtime(i, :)-log(rand)/lambda];i=i+1;endn=length(arrtime); % arrival times t_1,...t_n arrsubtr=arrtime-(0:n-1)'; % t_k-(k-1)arrmatrix=arrsubtr*ones(1,n);deptime=(1:n)+max(triu(arrmatrix)); % departure times%u_k=k+max(t_1,..,t_k-k +1)B=[ones(n,1) arrtime ; -ones(n,1) deptime'];Bsort=sortrows(B,2); % sort jumps in order jumps=Bsort(:,1);jumptimes=[0;Bsort(:,2)];systsize=[0;cumsum(jumps)]; % M/D/1 processsysttime=deptime-arrtime'; % system times% plot a histogram of system timeshist(systtime,30);④M/G/infinity系统%function [jumptimes, systsize] = simmginfty(tmax, lambda)% SIMMGINFTY simulate a M/G/infinity queueing system. Arrivals are% a homogeneous Poisson process of intensity lambda. Service times% Pareto distributed (can be modified).% [jumptimes, systsize] = simmginfty(tmax, lambda) %% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - times of state changes in the system% systsize - number of customers in system% See SIMSTMGINFTY, SIMGEOD1, SIMMM1, SIMMD1, SIMMG1.% set default parameter values if ommitedⅠ: nargin=0nargin=0;if (nargin==0)tmax=1500;lambda=1;end% generate Poisson arrivals% the number of points is Poisson-distributednpoints = poissrnd(lambda*tmax);% conditioned that number of points is N,% the points are uniformly distributedif (npoints>0)arrt = sort(rand(npoints, 1)*tmax);elsearrt = [];end% uncomment if not available POISSONRND% generate Poisson arrivals% arrt=-log(rand)/lambda;% i=1;% while (min(arrt(i,:))<=tmax)% arrt = [arrt; arrt(i, :)-log(rand)/lambda];% i=i+1;% end% npoints=length(arrt); % arrival times t_1,...,t_n% servt=50.*rand(n,1); % uniform service times s_1,...,s_kalpha = 1.5; % Pareto service times servt = rand^(-1/(alpha-1))-1; % stationary renewal processservt = [servt; rand(npoints-1,1).^(-1/alpha)-1]; servt = 10.*servt; % arbitrary choice of meandept = arrt+servt; % departure times% Output is system size process N.B = [ones(npoints, 1) arrt; -ones(npoints, 1) dept];Bsort = sortrows(B, 2); % sort jumps in orderjumps = Bsort(:, 1);jumptimes = [0; Bsort(:, 2)];systsize = [0; cumsum(jumps)]; % M/G/infinity system size% process stairs(jumptimes, systsize);xmax = max(systsize)+5;axis([0 tmax 0 xmax]);gridⅡ: nargin不为0时%function [jumptimes, systsize] = simmginfty(tmax, lambda)% SIMMGINFTY simulate a M/G/infinity queueing system. Arrivals are% a homogeneous Poisson process of intensity lambda. Service times% Pareto distributed (can be modified).% [jumptimes, systsize] = simmginfty(tmax, lambda) %% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - times of state changes in the system% systsize - number of customers in system% See SIMSTMGINFTY, SIMGEOD1, SIMMM1, SIMMD1, SIMMG1.% set default parameter values if ommitednargin=2;if (nargin==0)tmax=1500;lambda=1;end% generate Poisson arrivals% the number of points is Poisson-distributednpoints = poissrnd(lambda*tmax);% conditioned that number of points is N,% the points are uniformly distributedif (npoints>0)arrt = sort(rand(npoints, 1)*tmax);elsearrt = [];end% uncomment if not available POISSONRND% generate Poisson arrivals% arrt=-log(rand)/lambda;% i=1;% while (min(arrt(i,:))<=tmax)% arrt = [arrt; arrt(i, :)-log(rand)/lambda];% i=i+1;% end% npoints=length(arrt); % arrival times t_1,...,t_n% servt=50.*rand(n,1); % uniform service times s_1,...,s_kalpha = 1.5; % Pareto service times servt = rand^(-1/(alpha-1))-1; % stationary renewal process servt = [servt; rand(npoints-1,1).^(-1/alpha)-1]; servt = 10.*servt; % arbitrary choice of meandept = arrt+servt; % departure times% Output is system size process N.B = [ones(npoints, 1) arrt; -ones(npoints, 1) dept];Bsort = sortrows(B, 2); % sort jumps in orderjumps = Bsort(:, 1);jumptimes = [0; Bsort(:, 2)];systsize = [0; cumsum(jumps)]; % M/G/infinity system size% process stairs(jumptimes, systsize);xmax = max(systsize)+5;axis([0 tmax 0 xmax]);grid实验总结:通过这次实验,加深了我对随机过程这门课程的理解与认识,对各种随机过程的模拟有了更深刻的了解,同时也认识到模拟编程的重要性。

相关文档
最新文档