33
龚光鲁，钱敏平著 应用随机过程教程及其在 算法与智能计算中的应用
清华大学出版社,2003
第２章 随机样本生成法
1 一维随机数
随机变量（或随机向量）的样本简称为 随机数,由于在统计中常用的是独立样本列,
所以我们不妨假设随机数之间都是独立的．生成随机数的方法，也称为随机数的取样法
(Sampling）,
1,1 均匀随机变量的计算机模拟
定义２,１ 在 [0,1]上均匀分布的随机变量的独立样本称为 均匀随机数 ( ]1,0[U 随机数 ),
在计算机上产生的称之为,伪随机数,的数列,是一种具有非常长周期的,且能通过数理统计中的独立性与均匀性假设检验的数列,实践证明伪随机数是均匀随机数的一种可行的近似,这种伪随机数虽然并不是独立同分布的 ]1,0[U 随机变量的样本,而是在 [0,1]中取值的周期数列,但是由于它可以像均匀随机数一样地通过数理统计中的独立性与均匀性假设检验,而且它的周期非常长,以至在计算机实际运算过程中不会出现重复,所以在实际计算中它能很好地替代均匀随机数,
最普遍用以产生伪随机数的方法是同余法,典型的例子如下,
）周期约为 1036036131 102(2,1,)2(mod5===?+ nnnn yxyyy ；
）周期约为 1242042171 10(2,1,)2(mod5?+?=== nnnn yxyyy ；
）周期约为 1444104412 1041(2,1,0),2(mod===+=?++ nnnnn yxyyyyy,
（关于伪随机数,可参见：现代数学手册，随机数学卷，第 10 篇，孙嘉阳,石坚，丛树铮，徐映波编 蒙特卡罗法．华中科技大学出版社,2000 年）,
1,2 分布函数 )( xF 的随机数
命题２,２ （反函数方法） 分布函数为 )( xF 的独立随机变量列的样本，称为 )( xF 随机
数,若 )( xF 严格单调递增,x 是 均匀 随机数,则 ( )x1?F 是 )( xF 随机数,其中 1?F 为 F 的反函数,
证明 )())(())(( 1 xFxFPxFP =≤=≤? xx,?
命题２,３ 设随机变量 x 只取有限个值，其分布为
n
n
pp
xx
L
L
1
1~x,把 [0,1]
分为 n 个不交 子 区间,使第 i 个区间 iJ 的长度为 ip,任取均匀随机数 U,则
34
∑
=
n
i
Ji UIx i
1
)(
就是一 个 x 随机数（它的意思是：如果 U 落入 iJ,就取 ix=x ）,
在统计再抽样中的应用
在样本组中再抽样，或者由样本作的参数估计代替分布中的未知参数 后,所 得 到的分布的随机取样，统一称为 Bootstrap 方法,具体地说 有如下两种方法
（１） 非参数 Bootstrap 方法．设自一个未知方差的分布取样 nXX,,1 L （不是计算机模 拟取样，而是人工取样）如果要作方差的区间估计，就需要 知道 方差 估计
2∧
s 的方差
)(
2∧
sVar,一般 )(
2∧
sVar 很不好求，需要对它用 再 抽样进行估计．为此 可将 样本分布
nn
XX n
11
1
L
L
作为离散随机变量的分布,独立地取样 N 次，每次独立地取样 m 个,设 从第 k 次的 m 个样本值得到方差的估计 )(
2
Nkk ≤∧s,将此 N 个的平均记为
2∧
s,最后 用
2
22
1
2
)(11 ∧∧
=
∧?
= ∑
∧
sss k
N
kN
估计 )(
2∧
sVar,
此法可以用于一般未知参数的方差估计,
（２）参数 Bootstrap 方法,设自一个带有未知参数 ),( 1 lJJ L 的分布 ),,( 1 lxp JJ L
的样本 nXX,,1 L （不是计算机模拟取样，而是人工取样）得到未知参数的估计 )( l∧∧ JJ,，１ L
后，对分布 ),lxp ∧∧ JJ,，（ １ L 用计算机模拟取样．独立地取样 N 次，每次独立地取样 m 个．其它与（１）相同,
注意,计算机模拟取样只能对已知的分布施行，对于含未知参数的分布，只能作普通的人工取样．以上的两种再抽样方法，补充了人工取样采样量的限制．因为计算机模拟取样 既快速又经济,
1,3 正态随机数
)1,0(N 随机数称为 标准正态随机数,生成标准正态随机数有一个比反函数的方法更为简单的实践方法,就是利用中心极限定理,设 121,,hh L 为均匀随机数（它们是独立的）,由中心极限定理，可以认为 )1,0(6121 N≈?++= hhx L,即用 6121?++= hhx L 近似地作为标准正态随机数,在实际计算中 ih )121( ≤≤ i 们还应该用伪随机数代替,
命题２,４ (生成标准正态随机数的 Box - Muller 方法 ) 取两个独立的均匀随机数
35
21,hh,令
)2cos(ln2 211 phhx?=,
)2sin(ln2 212 phhx?=,
则 21,xx 为相互独立的标准正态随机数,
证明 令 ]1,0[~,21 Uhh 且独立,则 21,1 hh? 也是独立的 ]1,0[U 随机变量,于是由命题
２,2,~ln2)1( 111 hhr?=?= F 分布函数 2
2
1)(
r
erF=,]20[~2 2 pphJ,U?=,且相互独立,而 JrxJrx sin,cos 21 ==,易见它们是独立的标准正态随机数,
命题２,５ （生成标准正态随机数的 Marsaglia 方法 ） 设 ),( YX 为单位圆上的均匀
随机数,则
+
+?=
Y
X
YX
YX
22
22 )ln(2
h
x
10
01,
0
0~ N,
（提示 将直角坐标 ),( YX 转换为极坐标 ),( JR ）,
一般正态随机数的生成 若 x 为标准正态随机数,则显见 msx + 为 ),( 2smN 随机数,
1,4 Poisson 随机数
下述结论给出了利用伪随机数生成 Poisson 随机数的方法。
命题２,６ 设 L,,21 UU 是相互独立的 [0,1]均匀随机数,若
∏ ∏
+
= =
λ? ≤<
1
1 1
n
k
n
k
kk UeU,则定义 nN =,那么 λPoissonN ~,
证明 令 1exp~ln kk UT?=?,在指数流与 Poisson 过程的关系 ( 参见第 3 章 ) 中取参数为 1,取时间 t 为 l 即得,
1,5 混合分布随机数
对 于权重为 npp,,1 L (和为 1 的 n 个正数 ) 的 混合分布随机数,我们有
命题２,７ 设 )(,10],1,0[~ 10 nipttttUU iiin ≤=?=<<=?L,))(( nixFi ≤ 为分布函数,那么
∑
=
n
i
tt
i
i
i UIp
tUF
ii
1
],(
11 ~)()(
1 分布函数为 ii
n
i
Fp∑
=1
的混合分布,
证明
36
=≤?
∑
=
))()((
],(
1
11
1
xUIptUFP
ii tt
n
i i
i
i ]),(,)()(( 1],(
11
1
1 iitt
i
i
i
n
i
ttUxUIptUFP
ii?
=
∈≤?
∑
)())((
11
11 xFpxFptUtP
n
i
ii
n
i
iiii ∑∑
==
=+≤<=,?
在实际计算中,应该用伪随机数来取代均匀随机数 U,如果取到的伪随机数落在 ],( 1 ii tt?
中,则取 )( 11
i
i
i p
tUF,这样得到的数就是 ∑
=
n
i
iiFp
1
混合分布随机数,
这个方法常用在排队论中,在那里的典型情形是 混合指数分布 (有的书上称为 超指数分布 ),即 xi iexF l= 1)(,此时简单地有 )1ln(1)(1 yyF
i
i?λ?=
,于是计算变得非常简单而有效 (当然也 可 利用命题 2.2 通过反函数来得到混合指数分布随机数,但是计算量会增加很多,因为这个反函数并不简单 ),
1,6 Von Neuman 取舍原则
Von Neuman 取舍原则,
假定我们要生成密度为 )( xp 的随机数,为此取一个参考分布密度 )(0 xp,使它满足,
(1) )(0 xp 随机数容易生成,例如 )(0 xp 为正态密度,均匀密度,指数密度,及它们的混合密度；
(2) )(0 xp 与 )( xp 的取值范围差不多,且存在 C,使 )()( 0 xpCxp?≤,
那么，我们有
命题 2,８ 设随机变量 h 具有密度 )(0 xp,而随机变量 ]1,0[~ UU 且与 h 独立,则
∫
∞?
=≥≤
x
dvvpUCppxP )())( )(|(
0 h
hh,
证明 对 h 的取值用推广了的全概率公式（ ∫ == dxxgxAPAP )()|()( h ），得到
右左 ==
≤
≤
=
≥
≥≤
=
∫
∫
∫
∫
∞
∞?
∞?
∞
∞?
∞?
dyypC
dyypC
dyypyCp ypUP
dyypyCp ypUP
UCppP
UCppxP
xx
)(1
)(1
)())()((
)())()((
))( )((
))()(,(
0
0
0
0
0
0
h
h
h
hh
,?
取舍原则 (Rejection Principle)的具体作法是,
(1) 独立地生成 n 个独立的 )(0 xp 随机数 nhh,,1 L 与 n 个与之独立的 ]1,0[U 随机数
nUU,,1 L,
37
(2) 对于 L,2,1=i,如果有 i
i
i U
Cp
p ≥
)(
)(
0 h
h,就保留
ih,否则就舍弃 ih 。
由命题 2.８,所有这样保留下来的那些 ih 们就成为一系列独立的 )( xp 随机数 (当然个数会比 n 小很多 ),这种取舍方法称为 Von Neuman 取舍原则,
取舍原则可以改良为,
如果 )()( xhxp?= g,只要存在 C,使 )()( 0 xpCxh?≤,那么我们可以在取舍原则中用 )(xh 代替 )( xp,得到 )( xp 随机数,具体为：独立地生成 n 个独立的 )(0 xp 随机数
nhh,,1 L 与 n 个与之独立的独立 ]1,0[U 随机数 nUU,,1 L,如果 i
i
i U
Cp
h ≥
)(
)(
0 h
h,则保留
ih,
否则舍弃 ih,那么所有保留下的是相互独立的 )( xp 随机数,
证明 只需注意到这时有 )()( 0 xpCxp≤ g,并且 =? )()(
0 xCp
xp
g )(
)(
0 xCp
xh 即可,
注１ 一般地 g 需通过 ∫= dxxh )(1g 计算,其中的积分不易计算,但是上面的事实说明不必计算 g,即可以忽视这个常数因子．这就使取舍原则变得非常好用．取舍原则的优点是简单易行,可以适用于非常复杂的分布,
注２ 如果 )( xp 只在有限区间 ],[ ba 上不等于零,而且有界,那么 )(0 xp 就可取均匀分布 ],[ baU ; 如果 )( xp 只在右半直线不等于零,那么指数分布就可 以 是 )(0 xp 的一个选择 ;
如果 )( xp 在实直线上不等于零，且分布密度的尾部不大,则正态分布就可 以 是 )(0 xp 的一个选择 ; 如果 )( xp 在实直线上不等于零，且分布密度的尾部较大（重尾分布）,则 t 分布就可 以 是 )(0 xp 的一个选择 ; 如果 )( xp 具有多个峰,则混合正态分布或混合指数分布就可以是 )(0 xp 的一个选择,可见适当精心地选取 )(0 xp 是使计算省时的关键,
注３ 原则上取舍原则也适用于离散分布和多维密度，但是在多维密度的情形,)(0 xp
的选取并不容易,
注４ 样本生成的一个重要应用,就 是对于函数的一些积分,用 Monte Carlo 方法 (随机模拟算法 )给出估计,直观地可以猜测,采用的随机数的密度 )(0 xp 的形状与被积函数越像,
则估计的方差会越小,即效果越好,这种取样法称为 重要度采样 (Impotance Samling)( 确切的定义与方差的最小性证明可参见第 8 章 ),
38
1,7 Gamma 随机数与 Beta 随机数的生成
在 1≥a 时,设 ][a 是正数 a 的整部,可以证明 ),( laΓ 分布的密度与 )],([ laΓ 分布的密度的比例是有界函数,由此可以设计一个由指数分布随机数,通过 Von Neuman 取舍 原则得到 ),( laΓ 分布的随机数的方法,
在 1,≥ba 时，设 21,hh 分别为独立的 )1,(),1,( ba ΓΓ 随机数,则可以证明
21
1
hh
hV
+=
为 ),( baB 随机数,
关于 10 <<a 时 ),( laG 分布的取样方法，可以参见文献 [LK],
2,多维随机数
2,1 连续型多维随机数
对于已知的分布密度,我们可以利用条件密度,把生成多维随机数归结为生成一系列一维随机数,设随机向量 ),,( 1 dXX L 的密度为 ),,,( 21 dxxxf L,那么，我们有表达式
),,|()|()(),,,( 1112121
1?
= ddXd xxxfxxfxfxxxf LLL,
其中
1X
f 为 1X 的边缘密度,),,|( 11?kk xxxf L 为在已知 1111,, == kk xXxX L 条件下,
kX 的条件密度,于是可以先取一个 1Xf 随机数 1x ; 然后，在 1x 固定的情形下,生成一个
)|( 1xf? 随机数 2x ; 再在 21,xx 固定的情形下,生成一个 ),|( 21 xxf? 随机数 3x ; … 最后,
在 11,,?dxx L 固定的情形下,生成一个 ),,|( 11 dxxf L 随机数 dx,这样得到的
),,( 1 dxx L 就是随机向量 ),,( 1 dXX L 的一个随机数,当然在生成各个条件密度的随机数时,仍然可以使用 Von Neuman 取舍原则,注意在用取舍原则于一维分布取样时，可以忽略一个因子,
然而，在 d 较大时,更为常用的是 Gibbs 取样法 (Gibbs Sampler),它是一种基本的动态
Monte Carlo 方法,即 Markov 链 Monte Carlo 方法,在本书第 8 章中我们将介绍这种方 法,
在很多实际情形中,多维密度常常并不知道,更无法知道各个条件密度,这时最自然而粗糙的想法是用条件频率来近似条件密度,
2,2 离散型多维随机数
只要用概率函数 ),,(),,( 111 ddd xXxXPxxp ===? LL 代替密度函数,2,1段中的办法就自然适用,
2,3 多维正态随机数
39
设 ),,( 1 dXX L 服从 ),( ΣmN,其密度函数为 )()(2
1
2
1
2
1
||)2(
1)( mSm
Sp
= xx
d
T
exp,
其中 ( ) djiij ≤=,sS 为对称正定矩阵,由线性代数我们知道,可以找到下三角形矩阵 C,
),0( jicij <?=,使 TCC=Σ,作独立 )1,0(N 随机数 dxx,,1 L,令 Td ),,( 1 xxx L?=,则
),(~ Σ+ mmx NC,可见 mx +C 为 ),( ΣmN 随机数,于是求多维正态随机数的快速生成问题,就归结为快速计算下三角形矩阵 C 的问题,
2,4 多维 Beta 随机数 (Dirichlet 随机数 ) 的生成
在 1,,1 ≥kaa L 时,设 khh,,1 L 分别为独立的?ΓΓ )1,(,),1,( 1 kaa L 随机数,则可以
证明
,,(
1
1 L
L khh
h
++ )1 k
k
hh
h
++L
就是一个 Dirichlet ),,( 1 kaa L 随机数,
3．附录 － 用 Matalab 生成随机数
3,1 Matlab 语言的简单提示
Matlab 提供了强大功能的软件包，其用法可以参阅 Matlab 6.1(或 Matlab 6.5) 的 Help．为了使用方便，我们在本段列出一些基本的记号与最常见的语句,
一,Matlab 中的矩阵,
1．矩阵 A的分量用方括号括起来．同行的分量用逗号或空格分开,不同行的分量用分号或回车 (Enter)分开,空的方括号 [ ]表示没有分量的矩阵，即空集,
2,A(3,6)表示矩阵 A 的 (3,6)分量．反之,给出了它的值也可以调用 A (在 A 的行列不足 43× 时将自动增加行列,即除了在该分量处赋要求的值以外，将其它分量增补为 0),
３．矩阵运算除＋,－外,还有,*（乘）,' (转置 ),^ (幂 ),
\ (左除,BA\ 即 BA 1? ),/ (右除,BA/ 即 1?BA ),
常数视其情形可以看成任意阶的矩阵，如 A+a (A 的所有分量都加常数 a）,
４．可以用冒号对矩阵进行剪裁,如 A[3,:]表示取 A 的第 3 行,A[:,3]表示取 A 的第 3
列,A[1:3,:] 表示取 A 的 1,2,3 行,A[1:2:3,:] 表示取 A 的 1,3 两行．又用 A[3,:]=[ ]表示删除 A 的第 3 行,
５．两个矩阵可以拼接成一个,
６．特殊的矩阵（在相应的语句后加 ),( nm 表示其阶）有,zeros(0 矩阵）,ones（分量全
1 的矩阵）,eyes（对角为 1，其它为 0 的矩阵）,rand ([0,1]均匀随机数的矩阵 )等,
二 Matlab 中的数组
１,与矩阵的差别是数组不加括号．用冒号表示等差数组，即 ndm,,表示公差为 d,
40
起始于 m,终止于 n 的等差数组,而将 nm,1,简记为 nm,
２,数组间在运算符号前加一个小点表示按分量运算,有,,
*（乘）,.^（幂）,./ (除（即右除 )）,.\ (反除（即左除），即除的倒数 ),
３,几个数组可以用方括号拼接成矩阵,
４,单引号括起来的字符串可以作为数组,
又如语句 ),,( nbalinspace 生成起始于 a,终止于 b 的 n 个等差数组,
三．逻辑关系
１．一个关系的逻辑值只取两个值,0（表示＂错误＂），或 1（表示＂正确＂ )．所以，
关系 x==0 的逻辑值，就是 x=0 的示性函数,
２．主要的逻辑关系与逻辑运算有,
＝＝（逻辑相等））,> =（不小于）,< =（不大于）,~ =（不等于）,＆（与），｜（或）,
～（非）,
四．语句，函数与操作命令
１．一般的语句的形式是
变量 =表达式
(如用分号结束，则表示只运算不显示．如用回车结束则显示它的计算结果．如只有表达式，
则自动生成变量 ans，并显示 ans =),
２．两个语句间用逗号或分号,
３．特殊变量有 pi (圆周率 p ),i 或 j (虚数单位）,Inf （无穷大）,
NaN （不定值,0／ 0）,eps (主机中的最小的浮点数,有时用它代替０可以防止溢出，即用表达式 (x==0)*eps 代替 0,它表示在 x=0 时用 eps 代替 x),
４．函数的变量置于圆括号之中．常见的函数除对数，指数，三角函数外，还有,
atan（反正切）,abs（绝对值）,sqrt（开方）,round(4 舍 5 入取整 ),floor(负无穷方向取整 ),ceil(正无穷方向取整 ),fix(0 方向取整 ),imag（虚部）,rats(最近有理数 )
length（数组的长度）,prod（乘）,sort(数组的分量按递增排序 ),trace(求迹 ),triu (矩阵取上三角 ),tril (矩阵取下三角 ),size（矩阵的大小）,norm（模）,inv（取逆）,eig（求特征值）,poly(求特征多项式 ),expm(求矩阵的指数函数 ),cond(条件的数目 ),lu(矩阵的 LU 分解 ),
gr(矩阵的正交分解 ),svd(矩阵的奇异值分解 ),
feval（赋值语句，如 feval('atan',x) 等价于 atan(x),
５．用于循环的常见语句有,
pause（暂停）,
all,any,
for,end
while,end,
if,elseif,else,end
switch,case,otherwise,end
break(跳出循环 ),
６．操作命令（很多与 Unix 命令相同）,
help,lookfor,format short（显示到小数点后４位）,whos（显示当前工作区的变量）,
disp(x)(显示 x 的内容 ),clear,save,load,diary filename（建立新文件）,diary off（停止记录）,
path,dir,type,delete,quit 等,
41
７．画图命令有,
plot(x,y)（一个图）,plot(x,y1,x,y2)（两个图）,或 fplot('图形名 ',[xmin xmax ymin ymax])
或在图上加,hold on,第２个图的语句,hold off
线图,-,,(点线 ),-,(虚点线 ),－－ (波折 )
点 图,．,+,*
颜色,y（黄）,r（红）,g（绿）,b（兰）,w（白）,k(黑 ),m(紫 ),c(青 cyan)
例如,plot(x,y1,'b:',x,y2,'c-')
控制坐标范围用,axis auto,axis off
在图上标字：用 title('字串 '),或用 gtext('字串 ')再用鼠标点
分区画多个图：用 subplot(m,n,k) (在 nm× 个区中的第 k 个区中画图）
曲面：［ X,Y］＝ meshgrid(x,y)；（把 (x,y)改写为数组［ X,Y］，以便用于表达式中）；
Z＝ (X,Y)的表达式 ; mesh(X,Y,Z)
空间曲线图例如,plot3(sin(t),cos(t),t)
等高线图,contour(X,Y,Z,m) (m 条等高线 )。
３．２ Matlab 生成随机数的语句
第一种方法是用 random 语句，其一般形式为
y = random('分布的英文名 ',A1,A2,A3,m,n),
表示生成 m 行 n 列的 nm× 个参数为 ),,( 321 AAA 的该分布的随机数。例如,
(1) R = random('Normal',0,２,2,4),生成期望为 0,标准差为 1 的 (2 行 4 列 )2× 4 个正态随机数
(2) R = random('Poisson',1:6,1,6),依次生成参数为 1 到 6 的 (1 行 6 列 )6 个 Poisson 随机数
第二种方法是针对特殊的分布的语句,
一,几何分布随机数 （下面的 P,m 都可以是矩阵）
R = geornd(P) （生 成参数为 P 的几何随机数）
R = geornd(P,m) （生成参数为 P 的 m×１ 个几何随机数）
R = geornd(P,m,n) （生成参数为 P 的 m 行 n 列的 nm× 个几何随机数）
例如
(1) R = geornd(1,/ 2.^(1:6)) (生成参数依次为 1/2,1/2^2,到 1/2^6 的 6 个几何随机数 )
(2) R = geornd(0.01,[1 5]) (生成参数为 0.01 的（１行５列） 5 个几何随机数 ),
二,Beta 分布随机数
R = betarnd(A,B) （生成参数为 A,B的 Beta 随机数）
R = betarnd(A,B,m) （生成 m×１ 个数为 A,B的 Beta 随机数）
R = betarnd(A,B,m,n) （生成 m 行 n 列的 nm× 个数为 A,B的 Beta 随机数）,
三．正态随机数
R = normrnd(MU,SIGMA) （生成均值为 MU，标准差为 SIGMA 的正态随机数）
R = normrnd(MU,SIGMA,m) （生成 m×１ 个正态随机数）
R = normrnd(MU,SIGMA,m,n) （生成 m 行 n 列的 nm× 个正态随机数）
例如
(1) R = normrnd(0,1,[1 5]) 生成 5 个正态 (0,1)随机数
(2) R = normrnd([1 2 3;4 5 6],0.1,2,3) 生成期望依次为 [1,2,3;4,5,6],方差为 0.1的 6 个正态随
42
机数,
四．二项随机数：类似地有
R = binornd(N,P) R = binornd(N,P,m) R = binornd(N,p,m,n)
例如
n = 10:10:60; r1 = binornd(n,1./n) 或 r2 = binornd(n,1./n,[1 6]) （都生成参数分别为
)601,60(,),101,10( L 的６个二项随机数,
五．自由度为 V 的 2c 随机数,
R = chi2rnd(V) R = chi2rnd(V,m) R = chi2rnd(V,m,n)
六．期望为 MU 的指数随机数（即
MU
Exp 1 随机数）,
R = exprnd(MU) R = exprnd(MU,m) R = exprnd(MU,m,n)
七．自由度为 V1,V2 的 F 分布随机数,
R = frnd(V1,V2) R = frnd(V1,V2,m) R = frnd(V1,V2,m,n)
八,),( lAΓ 随机数,
R = gamrnd（ A,lambda） R = gamrnd（ A,lambda,m) R = gamrnd（ A,lambda,m,n)
九．超几何分布随机数,
R = hygernd(N,K,M) R = hygernd(N,K,M,m) R = hygernd(N,K,M,m,n)
十．对数正态分布随机数
R = lognrnd(MU,SIGMA) R = lognrnd(MU,SIGMA,m) R = lognrnd(MU,SIGMA,m,n)
十一．负二项随机数,
R = nbinrnd(r,p) R = nbinrnd(r,p,m) R = nbinrnd(r,p,m,n)
十二,Poisson 随机数,
R = poissrnd(lambda) R = poissrnd(lambda,m) R = poissrnd(lambda,m,n)
例如，以下 3 种表达有相同的含义,lambda = 2; R = poissrnd(lambda,1,10)
（或 R = poissrnd(lambda,[1 10]) 或 R = poissrnd(lambda(ones(1,10)))
十三,Rayleigh 随机数,
R = raylrnd(B) R = raylrnd(B,m) R = raylrnd(B,m,n)
十四,V 个自由度的 t 分布的随机数,
R = trnd(V) R = trnd(V,m) R = trnd(V,m,n)
43
十五．离散的均匀随机数,
R = unidrnd(N) R = unidrnd(N,m) R = unidrnd(N,m,n)
十六,[A,B]上均匀随机数
R = unifrnd(A,B) R = unifrnd(A,B,m) R = unifrnd(A,B,m,n)
例如 unifrnd(0,1:6)与 unifrnd(0,1:6,[1 6])都依次生成 [0,1]到 [0,6]的６个均匀随机数．,
十七,Weibull 随机数
R = weibrnd(A,B) R = weibrnd(A,B,m) R = weibrnd(A,B,m,n)
习题 2
1,完成命题２．５的证明,
2,证明
2
1
0 )()( ll
ll
a
la
x
xxp
+=
是一个分布密度,设 U 是一个均匀随机数,取
lah
1
1
= U
U,证明它是一个 )(
0 xp 随机数,
3,设 21,hh 分别为独立的 )1,(),1,( ba ΓΓ 随机数,证明
21
1
hh
hV
+= 为 ),( baB 随机数,
4,移位指数分布的密度为 )()( ),[)( xIexp aax ∞?λ?λ= 。 如何得到它的样本？
５．若随机变量 )1,0(~ Nx,)(~ 2 nch 且与 x 独立,则
n
h
x
服从 )(nt 分布（其密度为
2
12
)1(
)2(
)2 1(
),(
+?
+
+
=
n
n
x
nn
n
xnt
Gp
G
）．如何构造 )(nt 随机数？
6,如何得到 ),(),,( mnBn lΓ 随 机数？
7,Weibull 分布的密度的另一个形式为 )()()( ),0[)(1 xIeaxamxf
m
a
x
m
∞
=
．如何得到其随机数？
8,证明若 x 为 lexp 随机数,则 rae
lx
h = 为 Pareto 随机数,
9,若 x 为 Weibull 分布 ),,( aW l 随机数，证明 mxlbh += )ln( a 为极值分布的随机数,
44
10,求密度为 2
12
2
1
)1(
)2(
)2 1( +
+
Γ
+Γ
k
ak
xa
kk
k
p
的 ),( akt 分布的方差,又问它的随机数应该如何得到?
(注意 )1,(kt 分布就是 )(kt 分布 ),
11,给出下列随机数的取样程序,二项分布 )31,10(B,参数为 )31,5( 的负二项分布,)5,7,10(),,( =MKN 的超几何分布,各生成 50 个 独立样本,并求其均值,方差,直方图,
12,C,G,Park,T,Park 和 D.W,Shin 在 1996 的一个构造相关系数为给定值 r 的二维二值随机向量
),( 21 xx 的样本的方法如下,令
j
PoissonX j l~ )3,2,1( =j,且相互独立,适当地选取
)3,2,1( =jjl 可以使 )2,1)(( 3)0( =+= iXXI iix 满足要求,为此只需证明
(1) )2,1(1 10~ =
ipp iiix,其中 )2,1(
)( 3 == +? iep i
i
ll,
(2) )1(),( 32121?= lxx eppCov,
最后,由 rpp,,21 解出 )3,2,1( =jjl,请对此设计得到 ),( 21 xx 随机数的一个程序,
13,设 ][a 是正数 a 的整部,证明 ),( laΓ 分布的密度与 )],([ laΓ 分布的密度的比例是有界函数,由此
设计一个由指数分布随机数,通过 Von Neuman 的 取舍 原则得到 ),( laΓ 分布的随机数的方法,
14,设 ][a 是正数 a 的整部,证明 ),( baB 分布的密度与 ])[],([ baB 分布的密度的 比例是有界函数,
由此设计一个由指数分布随机数,通过 Von Neuman 的 取舍 原则得到 ),( baB 分布的随机数的方法,
15,若 U 是［ 0,1 ］均匀随机数,则 [ p)-ln(1lnU ] +1 （整数部分）是参数为 p 的几何分布的随机数。
龚光鲁，钱敏平著 应用随机过程教程及其在 算法与智能计算中的应用
清华大学出版社,2003
第２章 随机样本生成法
1 一维随机数
随机变量（或随机向量）的样本简称为 随机数,由于在统计中常用的是独立样本列,
所以我们不妨假设随机数之间都是独立的．生成随机数的方法，也称为随机数的取样法
(Sampling）,
1,1 均匀随机变量的计算机模拟
定义２,１ 在 [0,1]上均匀分布的随机变量的独立样本称为 均匀随机数 ( ]1,0[U 随机数 ),
在计算机上产生的称之为,伪随机数,的数列,是一种具有非常长周期的,且能通过数理统计中的独立性与均匀性假设检验的数列,实践证明伪随机数是均匀随机数的一种可行的近似,这种伪随机数虽然并不是独立同分布的 ]1,0[U 随机变量的样本,而是在 [0,1]中取值的周期数列,但是由于它可以像均匀随机数一样地通过数理统计中的独立性与均匀性假设检验,而且它的周期非常长,以至在计算机实际运算过程中不会出现重复,所以在实际计算中它能很好地替代均匀随机数,
最普遍用以产生伪随机数的方法是同余法,典型的例子如下,
）周期约为 1036036131 102(2,1,)2(mod5===?+ nnnn yxyyy ；
）周期约为 1242042171 10(2,1,)2(mod5?+?=== nnnn yxyyy ；
）周期约为 1444104412 1041(2,1,0),2(mod===+=?++ nnnnn yxyyyyy,
（关于伪随机数,可参见：现代数学手册，随机数学卷，第 10 篇，孙嘉阳,石坚，丛树铮，徐映波编 蒙特卡罗法．华中科技大学出版社,2000 年）,
1,2 分布函数 )( xF 的随机数
命题２,２ （反函数方法） 分布函数为 )( xF 的独立随机变量列的样本，称为 )( xF 随机
数,若 )( xF 严格单调递增,x 是 均匀 随机数,则 ( )x1?F 是 )( xF 随机数,其中 1?F 为 F 的反函数,
证明 )())(())(( 1 xFxFPxFP =≤=≤? xx,?
命题２,３ 设随机变量 x 只取有限个值，其分布为
n
n
pp
xx
L
L
1
1~x,把 [0,1]
分为 n 个不交 子 区间,使第 i 个区间 iJ 的长度为 ip,任取均匀随机数 U,则
34
∑
=
n
i
Ji UIx i
1
)(
就是一 个 x 随机数（它的意思是：如果 U 落入 iJ,就取 ix=x ）,
在统计再抽样中的应用
在样本组中再抽样，或者由样本作的参数估计代替分布中的未知参数 后,所 得 到的分布的随机取样，统一称为 Bootstrap 方法,具体地说 有如下两种方法
（１） 非参数 Bootstrap 方法．设自一个未知方差的分布取样 nXX,,1 L （不是计算机模 拟取样，而是人工取样）如果要作方差的区间估计，就需要 知道 方差 估计
2∧
s 的方差
)(
2∧
sVar,一般 )(
2∧
sVar 很不好求，需要对它用 再 抽样进行估计．为此 可将 样本分布
nn
XX n
11
1
L
L
作为离散随机变量的分布,独立地取样 N 次，每次独立地取样 m 个,设 从第 k 次的 m 个样本值得到方差的估计 )(
2
Nkk ≤∧s,将此 N 个的平均记为
2∧
s,最后 用
2
22
1
2
)(11 ∧∧
=
∧?
= ∑
∧
sss k
N
kN
估计 )(
2∧
sVar,
此法可以用于一般未知参数的方差估计,
（２）参数 Bootstrap 方法,设自一个带有未知参数 ),( 1 lJJ L 的分布 ),,( 1 lxp JJ L
的样本 nXX,,1 L （不是计算机模拟取样，而是人工取样）得到未知参数的估计 )( l∧∧ JJ,，１ L
后，对分布 ),lxp ∧∧ JJ,，（ １ L 用计算机模拟取样．独立地取样 N 次，每次独立地取样 m 个．其它与（１）相同,
注意,计算机模拟取样只能对已知的分布施行，对于含未知参数的分布，只能作普通的人工取样．以上的两种再抽样方法，补充了人工取样采样量的限制．因为计算机模拟取样 既快速又经济,
1,3 正态随机数
)1,0(N 随机数称为 标准正态随机数,生成标准正态随机数有一个比反函数的方法更为简单的实践方法,就是利用中心极限定理,设 121,,hh L 为均匀随机数（它们是独立的）,由中心极限定理，可以认为 )1,0(6121 N≈?++= hhx L,即用 6121?++= hhx L 近似地作为标准正态随机数,在实际计算中 ih )121( ≤≤ i 们还应该用伪随机数代替,
命题２,４ (生成标准正态随机数的 Box - Muller 方法 ) 取两个独立的均匀随机数
35
21,hh,令
)2cos(ln2 211 phhx?=,
)2sin(ln2 212 phhx?=,
则 21,xx 为相互独立的标准正态随机数,
证明 令 ]1,0[~,21 Uhh 且独立,则 21,1 hh? 也是独立的 ]1,0[U 随机变量,于是由命题
２,2,~ln2)1( 111 hhr?=?= F 分布函数 2
2
1)(
r
erF=,]20[~2 2 pphJ,U?=,且相互独立,而 JrxJrx sin,cos 21 ==,易见它们是独立的标准正态随机数,
命题２,５ （生成标准正态随机数的 Marsaglia 方法 ） 设 ),( YX 为单位圆上的均匀
随机数,则
+
+?=
Y
X
YX
YX
22
22 )ln(2
h
x
10
01,
0
0~ N,
（提示 将直角坐标 ),( YX 转换为极坐标 ),( JR ）,
一般正态随机数的生成 若 x 为标准正态随机数,则显见 msx + 为 ),( 2smN 随机数,
1,4 Poisson 随机数
下述结论给出了利用伪随机数生成 Poisson 随机数的方法。
命题２,６ 设 L,,21 UU 是相互独立的 [0,1]均匀随机数,若
∏ ∏
+
= =
λ? ≤<
1
1 1
n
k
n
k
kk UeU,则定义 nN =,那么 λPoissonN ~,
证明 令 1exp~ln kk UT?=?,在指数流与 Poisson 过程的关系 ( 参见第 3 章 ) 中取参数为 1,取时间 t 为 l 即得,
1,5 混合分布随机数
对 于权重为 npp,,1 L (和为 1 的 n 个正数 ) 的 混合分布随机数,我们有
命题２,７ 设 )(,10],1,0[~ 10 nipttttUU iiin ≤=?=<<=?L,))(( nixFi ≤ 为分布函数,那么
∑
=
n
i
tt
i
i
i UIp
tUF
ii
1
],(
11 ~)()(
1 分布函数为 ii
n
i
Fp∑
=1
的混合分布,
证明
36
=≤?
∑
=
))()((
],(
1
11
1
xUIptUFP
ii tt
n
i i
i
i ]),(,)()(( 1],(
11
1
1 iitt
i
i
i
n
i
ttUxUIptUFP
ii?
=
∈≤?
∑
)())((
11
11 xFpxFptUtP
n
i
ii
n
i
iiii ∑∑
==
=+≤<=,?
在实际计算中,应该用伪随机数来取代均匀随机数 U,如果取到的伪随机数落在 ],( 1 ii tt?
中,则取 )( 11
i
i
i p
tUF,这样得到的数就是 ∑
=
n
i
iiFp
1
混合分布随机数,
这个方法常用在排队论中,在那里的典型情形是 混合指数分布 (有的书上称为 超指数分布 ),即 xi iexF l= 1)(,此时简单地有 )1ln(1)(1 yyF
i
i?λ?=
,于是计算变得非常简单而有效 (当然也 可 利用命题 2.2 通过反函数来得到混合指数分布随机数,但是计算量会增加很多,因为这个反函数并不简单 ),
1,6 Von Neuman 取舍原则
Von Neuman 取舍原则,
假定我们要生成密度为 )( xp 的随机数,为此取一个参考分布密度 )(0 xp,使它满足,
(1) )(0 xp 随机数容易生成,例如 )(0 xp 为正态密度,均匀密度,指数密度,及它们的混合密度；
(2) )(0 xp 与 )( xp 的取值范围差不多,且存在 C,使 )()( 0 xpCxp?≤,
那么，我们有
命题 2,８ 设随机变量 h 具有密度 )(0 xp,而随机变量 ]1,0[~ UU 且与 h 独立,则
∫
∞?
=≥≤
x
dvvpUCppxP )())( )(|(
0 h
hh,
证明 对 h 的取值用推广了的全概率公式（ ∫ == dxxgxAPAP )()|()( h ），得到
右左 ==
≤
≤
=
≥
≥≤
=
∫
∫
∫
∫
∞
∞?
∞?
∞
∞?
∞?
dyypC
dyypC
dyypyCp ypUP
dyypyCp ypUP
UCppP
UCppxP
xx
)(1
)(1
)())()((
)())()((
))( )((
))()(,(
0
0
0
0
0
0
h
h
h
hh
,?
取舍原则 (Rejection Principle)的具体作法是,
(1) 独立地生成 n 个独立的 )(0 xp 随机数 nhh,,1 L 与 n 个与之独立的 ]1,0[U 随机数
nUU,,1 L,
37
(2) 对于 L,2,1=i,如果有 i
i
i U
Cp
p ≥
)(
)(
0 h
h,就保留
ih,否则就舍弃 ih 。
由命题 2.８,所有这样保留下来的那些 ih 们就成为一系列独立的 )( xp 随机数 (当然个数会比 n 小很多 ),这种取舍方法称为 Von Neuman 取舍原则,
取舍原则可以改良为,
如果 )()( xhxp?= g,只要存在 C,使 )()( 0 xpCxh?≤,那么我们可以在取舍原则中用 )(xh 代替 )( xp,得到 )( xp 随机数,具体为：独立地生成 n 个独立的 )(0 xp 随机数
nhh,,1 L 与 n 个与之独立的独立 ]1,0[U 随机数 nUU,,1 L,如果 i
i
i U
Cp
h ≥
)(
)(
0 h
h,则保留
ih,
否则舍弃 ih,那么所有保留下的是相互独立的 )( xp 随机数,
证明 只需注意到这时有 )()( 0 xpCxp≤ g,并且 =? )()(
0 xCp
xp
g )(
)(
0 xCp
xh 即可,
注１ 一般地 g 需通过 ∫= dxxh )(1g 计算,其中的积分不易计算,但是上面的事实说明不必计算 g,即可以忽视这个常数因子．这就使取舍原则变得非常好用．取舍原则的优点是简单易行,可以适用于非常复杂的分布,
注２ 如果 )( xp 只在有限区间 ],[ ba 上不等于零,而且有界,那么 )(0 xp 就可取均匀分布 ],[ baU ; 如果 )( xp 只在右半直线不等于零,那么指数分布就可 以 是 )(0 xp 的一个选择 ;
如果 )( xp 在实直线上不等于零，且分布密度的尾部不大,则正态分布就可 以 是 )(0 xp 的一个选择 ; 如果 )( xp 在实直线上不等于零，且分布密度的尾部较大（重尾分布）,则 t 分布就可 以 是 )(0 xp 的一个选择 ; 如果 )( xp 具有多个峰,则混合正态分布或混合指数分布就可以是 )(0 xp 的一个选择,可见适当精心地选取 )(0 xp 是使计算省时的关键,
注３ 原则上取舍原则也适用于离散分布和多维密度，但是在多维密度的情形,)(0 xp
的选取并不容易,
注４ 样本生成的一个重要应用,就 是对于函数的一些积分,用 Monte Carlo 方法 (随机模拟算法 )给出估计,直观地可以猜测,采用的随机数的密度 )(0 xp 的形状与被积函数越像,
则估计的方差会越小,即效果越好,这种取样法称为 重要度采样 (Impotance Samling)( 确切的定义与方差的最小性证明可参见第 8 章 ),
38
1,7 Gamma 随机数与 Beta 随机数的生成
在 1≥a 时,设 ][a 是正数 a 的整部,可以证明 ),( laΓ 分布的密度与 )],([ laΓ 分布的密度的比例是有界函数,由此可以设计一个由指数分布随机数,通过 Von Neuman 取舍 原则得到 ),( laΓ 分布的随机数的方法,
在 1,≥ba 时，设 21,hh 分别为独立的 )1,(),1,( ba ΓΓ 随机数,则可以证明
21
1
hh
hV
+=
为 ),( baB 随机数,
关于 10 <<a 时 ),( laG 分布的取样方法，可以参见文献 [LK],
2,多维随机数
2,1 连续型多维随机数
对于已知的分布密度,我们可以利用条件密度,把生成多维随机数归结为生成一系列一维随机数,设随机向量 ),,( 1 dXX L 的密度为 ),,,( 21 dxxxf L,那么，我们有表达式
),,|()|()(),,,( 1112121
1?
= ddXd xxxfxxfxfxxxf LLL,
其中
1X
f 为 1X 的边缘密度,),,|( 11?kk xxxf L 为在已知 1111,, == kk xXxX L 条件下,
kX 的条件密度,于是可以先取一个 1Xf 随机数 1x ; 然后，在 1x 固定的情形下,生成一个
)|( 1xf? 随机数 2x ; 再在 21,xx 固定的情形下,生成一个 ),|( 21 xxf? 随机数 3x ; … 最后,
在 11,,?dxx L 固定的情形下,生成一个 ),,|( 11 dxxf L 随机数 dx,这样得到的
),,( 1 dxx L 就是随机向量 ),,( 1 dXX L 的一个随机数,当然在生成各个条件密度的随机数时,仍然可以使用 Von Neuman 取舍原则,注意在用取舍原则于一维分布取样时，可以忽略一个因子,
然而，在 d 较大时,更为常用的是 Gibbs 取样法 (Gibbs Sampler),它是一种基本的动态
Monte Carlo 方法,即 Markov 链 Monte Carlo 方法,在本书第 8 章中我们将介绍这种方 法,
在很多实际情形中,多维密度常常并不知道,更无法知道各个条件密度,这时最自然而粗糙的想法是用条件频率来近似条件密度,
2,2 离散型多维随机数
只要用概率函数 ),,(),,( 111 ddd xXxXPxxp ===? LL 代替密度函数,2,1段中的办法就自然适用,
2,3 多维正态随机数
39
设 ),,( 1 dXX L 服从 ),( ΣmN,其密度函数为 )()(2
1
2
1
2
1
||)2(
1)( mSm
Sp
= xx
d
T
exp,
其中 ( ) djiij ≤=,sS 为对称正定矩阵,由线性代数我们知道,可以找到下三角形矩阵 C,
),0( jicij <?=,使 TCC=Σ,作独立 )1,0(N 随机数 dxx,,1 L,令 Td ),,( 1 xxx L?=,则
),(~ Σ+ mmx NC,可见 mx +C 为 ),( ΣmN 随机数,于是求多维正态随机数的快速生成问题,就归结为快速计算下三角形矩阵 C 的问题,
2,4 多维 Beta 随机数 (Dirichlet 随机数 ) 的生成
在 1,,1 ≥kaa L 时,设 khh,,1 L 分别为独立的?ΓΓ )1,(,),1,( 1 kaa L 随机数,则可以
证明
,,(
1
1 L
L khh
h
++ )1 k
k
hh
h
++L
就是一个 Dirichlet ),,( 1 kaa L 随机数,
3．附录 － 用 Matalab 生成随机数
3,1 Matlab 语言的简单提示
Matlab 提供了强大功能的软件包，其用法可以参阅 Matlab 6.1(或 Matlab 6.5) 的 Help．为了使用方便，我们在本段列出一些基本的记号与最常见的语句,
一,Matlab 中的矩阵,
1．矩阵 A的分量用方括号括起来．同行的分量用逗号或空格分开,不同行的分量用分号或回车 (Enter)分开,空的方括号 [ ]表示没有分量的矩阵，即空集,
2,A(3,6)表示矩阵 A 的 (3,6)分量．反之,给出了它的值也可以调用 A (在 A 的行列不足 43× 时将自动增加行列,即除了在该分量处赋要求的值以外，将其它分量增补为 0),
３．矩阵运算除＋,－外,还有,*（乘）,' (转置 ),^ (幂 ),
\ (左除,BA\ 即 BA 1? ),/ (右除,BA/ 即 1?BA ),
常数视其情形可以看成任意阶的矩阵，如 A+a (A 的所有分量都加常数 a）,
４．可以用冒号对矩阵进行剪裁,如 A[3,:]表示取 A 的第 3 行,A[:,3]表示取 A 的第 3
列,A[1:3,:] 表示取 A 的 1,2,3 行,A[1:2:3,:] 表示取 A 的 1,3 两行．又用 A[3,:]=[ ]表示删除 A 的第 3 行,
５．两个矩阵可以拼接成一个,
６．特殊的矩阵（在相应的语句后加 ),( nm 表示其阶）有,zeros(0 矩阵）,ones（分量全
1 的矩阵）,eyes（对角为 1，其它为 0 的矩阵）,rand ([0,1]均匀随机数的矩阵 )等,
二 Matlab 中的数组
１,与矩阵的差别是数组不加括号．用冒号表示等差数组，即 ndm,,表示公差为 d,
40
起始于 m,终止于 n 的等差数组,而将 nm,1,简记为 nm,
２,数组间在运算符号前加一个小点表示按分量运算,有,,
*（乘）,.^（幂）,./ (除（即右除 )）,.\ (反除（即左除），即除的倒数 ),
３,几个数组可以用方括号拼接成矩阵,
４,单引号括起来的字符串可以作为数组,
又如语句 ),,( nbalinspace 生成起始于 a,终止于 b 的 n 个等差数组,
三．逻辑关系
１．一个关系的逻辑值只取两个值,0（表示＂错误＂），或 1（表示＂正确＂ )．所以，
关系 x==0 的逻辑值，就是 x=0 的示性函数,
２．主要的逻辑关系与逻辑运算有,
＝＝（逻辑相等））,> =（不小于）,< =（不大于）,~ =（不等于）,＆（与），｜（或）,
～（非）,
四．语句，函数与操作命令
１．一般的语句的形式是
变量 =表达式
(如用分号结束，则表示只运算不显示．如用回车结束则显示它的计算结果．如只有表达式，
则自动生成变量 ans，并显示 ans =),
２．两个语句间用逗号或分号,
３．特殊变量有 pi (圆周率 p ),i 或 j (虚数单位）,Inf （无穷大）,
NaN （不定值,0／ 0）,eps (主机中的最小的浮点数,有时用它代替０可以防止溢出，即用表达式 (x==0)*eps 代替 0,它表示在 x=0 时用 eps 代替 x),
４．函数的变量置于圆括号之中．常见的函数除对数，指数，三角函数外，还有,
atan（反正切）,abs（绝对值）,sqrt（开方）,round(4 舍 5 入取整 ),floor(负无穷方向取整 ),ceil(正无穷方向取整 ),fix(0 方向取整 ),imag（虚部）,rats(最近有理数 )
length（数组的长度）,prod（乘）,sort(数组的分量按递增排序 ),trace(求迹 ),triu (矩阵取上三角 ),tril (矩阵取下三角 ),size（矩阵的大小）,norm（模）,inv（取逆）,eig（求特征值）,poly(求特征多项式 ),expm(求矩阵的指数函数 ),cond(条件的数目 ),lu(矩阵的 LU 分解 ),
gr(矩阵的正交分解 ),svd(矩阵的奇异值分解 ),
feval（赋值语句，如 feval('atan',x) 等价于 atan(x),
５．用于循环的常见语句有,
pause（暂停）,
all,any,
for,end
while,end,
if,elseif,else,end
switch,case,otherwise,end
break(跳出循环 ),
６．操作命令（很多与 Unix 命令相同）,
help,lookfor,format short（显示到小数点后４位）,whos（显示当前工作区的变量）,
disp(x)(显示 x 的内容 ),clear,save,load,diary filename（建立新文件）,diary off（停止记录）,
path,dir,type,delete,quit 等,
41
７．画图命令有,
plot(x,y)（一个图）,plot(x,y1,x,y2)（两个图）,或 fplot('图形名 ',[xmin xmax ymin ymax])
或在图上加,hold on,第２个图的语句,hold off
线图,-,,(点线 ),-,(虚点线 ),－－ (波折 )
点 图,．,+,*
颜色,y（黄）,r（红）,g（绿）,b（兰）,w（白）,k(黑 ),m(紫 ),c(青 cyan)
例如,plot(x,y1,'b:',x,y2,'c-')
控制坐标范围用,axis auto,axis off
在图上标字：用 title('字串 '),或用 gtext('字串 ')再用鼠标点
分区画多个图：用 subplot(m,n,k) (在 nm× 个区中的第 k 个区中画图）
曲面：［ X,Y］＝ meshgrid(x,y)；（把 (x,y)改写为数组［ X,Y］，以便用于表达式中）；
Z＝ (X,Y)的表达式 ; mesh(X,Y,Z)
空间曲线图例如,plot3(sin(t),cos(t),t)
等高线图,contour(X,Y,Z,m) (m 条等高线 )。
３．２ Matlab 生成随机数的语句
第一种方法是用 random 语句，其一般形式为
y = random('分布的英文名 ',A1,A2,A3,m,n),
表示生成 m 行 n 列的 nm× 个参数为 ),,( 321 AAA 的该分布的随机数。例如,
(1) R = random('Normal',0,２,2,4),生成期望为 0,标准差为 1 的 (2 行 4 列 )2× 4 个正态随机数
(2) R = random('Poisson',1:6,1,6),依次生成参数为 1 到 6 的 (1 行 6 列 )6 个 Poisson 随机数
第二种方法是针对特殊的分布的语句,
一,几何分布随机数 （下面的 P,m 都可以是矩阵）
R = geornd(P) （生 成参数为 P 的几何随机数）
R = geornd(P,m) （生成参数为 P 的 m×１ 个几何随机数）
R = geornd(P,m,n) （生成参数为 P 的 m 行 n 列的 nm× 个几何随机数）
例如
(1) R = geornd(1,/ 2.^(1:6)) (生成参数依次为 1/2,1/2^2,到 1/2^6 的 6 个几何随机数 )
(2) R = geornd(0.01,[1 5]) (生成参数为 0.01 的（１行５列） 5 个几何随机数 ),
二,Beta 分布随机数
R = betarnd(A,B) （生成参数为 A,B的 Beta 随机数）
R = betarnd(A,B,m) （生成 m×１ 个数为 A,B的 Beta 随机数）
R = betarnd(A,B,m,n) （生成 m 行 n 列的 nm× 个数为 A,B的 Beta 随机数）,
三．正态随机数
R = normrnd(MU,SIGMA) （生成均值为 MU，标准差为 SIGMA 的正态随机数）
R = normrnd(MU,SIGMA,m) （生成 m×１ 个正态随机数）
R = normrnd(MU,SIGMA,m,n) （生成 m 行 n 列的 nm× 个正态随机数）
例如
(1) R = normrnd(0,1,[1 5]) 生成 5 个正态 (0,1)随机数
(2) R = normrnd([1 2 3;4 5 6],0.1,2,3) 生成期望依次为 [1,2,3;4,5,6],方差为 0.1的 6 个正态随
42
机数,
四．二项随机数：类似地有
R = binornd(N,P) R = binornd(N,P,m) R = binornd(N,p,m,n)
例如
n = 10:10:60; r1 = binornd(n,1./n) 或 r2 = binornd(n,1./n,[1 6]) （都生成参数分别为
)601,60(,),101,10( L 的６个二项随机数,
五．自由度为 V 的 2c 随机数,
R = chi2rnd(V) R = chi2rnd(V,m) R = chi2rnd(V,m,n)
六．期望为 MU 的指数随机数（即
MU
Exp 1 随机数）,
R = exprnd(MU) R = exprnd(MU,m) R = exprnd(MU,m,n)
七．自由度为 V1,V2 的 F 分布随机数,
R = frnd(V1,V2) R = frnd(V1,V2,m) R = frnd(V1,V2,m,n)
八,),( lAΓ 随机数,
R = gamrnd（ A,lambda） R = gamrnd（ A,lambda,m) R = gamrnd（ A,lambda,m,n)
九．超几何分布随机数,
R = hygernd(N,K,M) R = hygernd(N,K,M,m) R = hygernd(N,K,M,m,n)
十．对数正态分布随机数
R = lognrnd(MU,SIGMA) R = lognrnd(MU,SIGMA,m) R = lognrnd(MU,SIGMA,m,n)
十一．负二项随机数,
R = nbinrnd(r,p) R = nbinrnd(r,p,m) R = nbinrnd(r,p,m,n)
十二,Poisson 随机数,
R = poissrnd(lambda) R = poissrnd(lambda,m) R = poissrnd(lambda,m,n)
例如，以下 3 种表达有相同的含义,lambda = 2; R = poissrnd(lambda,1,10)
（或 R = poissrnd(lambda,[1 10]) 或 R = poissrnd(lambda(ones(1,10)))
十三,Rayleigh 随机数,
R = raylrnd(B) R = raylrnd(B,m) R = raylrnd(B,m,n)
十四,V 个自由度的 t 分布的随机数,
R = trnd(V) R = trnd(V,m) R = trnd(V,m,n)
43
十五．离散的均匀随机数,
R = unidrnd(N) R = unidrnd(N,m) R = unidrnd(N,m,n)
十六,[A,B]上均匀随机数
R = unifrnd(A,B) R = unifrnd(A,B,m) R = unifrnd(A,B,m,n)
例如 unifrnd(0,1:6)与 unifrnd(0,1:6,[1 6])都依次生成 [0,1]到 [0,6]的６个均匀随机数．,
十七,Weibull 随机数
R = weibrnd(A,B) R = weibrnd(A,B,m) R = weibrnd(A,B,m,n)
习题 2
1,完成命题２．５的证明,
2,证明
2
1
0 )()( ll
ll
a
la
x
xxp
+=
是一个分布密度,设 U 是一个均匀随机数,取
lah
1
1
= U
U,证明它是一个 )(
0 xp 随机数,
3,设 21,hh 分别为独立的 )1,(),1,( ba ΓΓ 随机数,证明
21
1
hh
hV
+= 为 ),( baB 随机数,
4,移位指数分布的密度为 )()( ),[)( xIexp aax ∞?λ?λ= 。 如何得到它的样本？
５．若随机变量 )1,0(~ Nx,)(~ 2 nch 且与 x 独立,则
n
h
x
服从 )(nt 分布（其密度为
2
12
)1(
)2(
)2 1(
),(
+?
+
+
=
n
n
x
nn
n
xnt
Gp
G
）．如何构造 )(nt 随机数？
6,如何得到 ),(),,( mnBn lΓ 随 机数？
7,Weibull 分布的密度的另一个形式为 )()()( ),0[)(1 xIeaxamxf
m
a
x
m
∞
=
．如何得到其随机数？
8,证明若 x 为 lexp 随机数,则 rae
lx
h = 为 Pareto 随机数,
9,若 x 为 Weibull 分布 ),,( aW l 随机数，证明 mxlbh += )ln( a 为极值分布的随机数,
44
10,求密度为 2
12
2
1
)1(
)2(
)2 1( +
+
Γ
+Γ
k
ak
xa
kk
k
p
的 ),( akt 分布的方差,又问它的随机数应该如何得到?
(注意 )1,(kt 分布就是 )(kt 分布 ),
11,给出下列随机数的取样程序,二项分布 )31,10(B,参数为 )31,5( 的负二项分布,)5,7,10(),,( =MKN 的超几何分布,各生成 50 个 独立样本,并求其均值,方差,直方图,
12,C,G,Park,T,Park 和 D.W,Shin 在 1996 的一个构造相关系数为给定值 r 的二维二值随机向量
),( 21 xx 的样本的方法如下,令
j
PoissonX j l~ )3,2,1( =j,且相互独立,适当地选取
)3,2,1( =jjl 可以使 )2,1)(( 3)0( =+= iXXI iix 满足要求,为此只需证明
(1) )2,1(1 10~ =
ipp iiix,其中 )2,1(
)( 3 == +? iep i
i
ll,
(2) )1(),( 32121?= lxx eppCov,
最后,由 rpp,,21 解出 )3,2,1( =jjl,请对此设计得到 ),( 21 xx 随机数的一个程序,
13,设 ][a 是正数 a 的整部,证明 ),( laΓ 分布的密度与 )],([ laΓ 分布的密度的比例是有界函数,由此
设计一个由指数分布随机数,通过 Von Neuman 的 取舍 原则得到 ),( laΓ 分布的随机数的方法,
14,设 ][a 是正数 a 的整部,证明 ),( baB 分布的密度与 ])[],([ baB 分布的密度的 比例是有界函数,
由此设计一个由指数分布随机数,通过 Von Neuman 的 取舍 原则得到 ),( baB 分布的随机数的方法,
15,若 U 是［ 0,1 ］均匀随机数,则 [ p)-ln(1lnU ] +1 （整数部分）是参数为 p 的几何分布的随机数。
