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: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
③%三维随机游走ranwalk3d 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
(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);
②二维Brown
npoints = 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;
③三维Brown
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;
④%泊松过程的模拟、检验及参数估计syms Un X S;
n=10;%生成n*n个随机数
r=1;%参数
temp=0;
tem=0;
Un=rand(n,1);%共产生n*n个随机数for i=1:1:n
X(i)=-log(Un(i))/r;
end
X=subs(X);
for i=1:1:n
for j=1:1:i
temp=temp+X(j);
end
S(i)=temp;
temp=0;
end
S=subs(S);
%检验泊松过程使用第四条for i=1:1:n
tem=tem+S(i);
end
sigmaN=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)
