中国科学技术大学：《实用统计软件》课程课件讲义（统计计算与软件）第四讲 随机数产生方法

Lecture4:随机数产生方法 张伟平 Sunday 11th October,2009

Lecture 4: ëÅÍ)ê{ ‹ï² Sunday 11th October, 2009

Contents 1 Methods for Generating Random Variables y 1.1 Generating Uniform(0,1)random number 44 1.2 Random Generators of Common Probability Distribution 3 l.2.1 The Inverse Transform Method...···.···· 6 l.2.2 The Acceptance-Rejection Method...·.·.··· 16 1.2.3 Transformation Methods................20 1.2.4 Sums and Mixtures..................·23 l.3 Multivariate Distribution···· 29 1.3.1 Multivariate Normal Distribution... 29 1.3.2 Mixtures of Multivariate Normals. 44 36 1.3.3 Vishart Distribution....·.。...。·· 38 l.3.4 Uniform Distribution on the d-Sphere...··.· 39 1.4 Stochastic Process.···.···.··.·····.·.···42 Previous Next First Last Back Forward 1

Contents 1 Methods for Generating Random Variables 1 1.1 Generating Uniform(0,1) random number . . . . . . . . . . 1 1.2 Random Generators of Common Probability Distribution in R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 The Inverse Transform Method . . . . . . . . . . . . 6 1.2.2 The Acceptance-Rejection Method . . . . . . . . . . 16 1.2.3 Transformation Methods . . . . . . . . . . . . . . . . 20 1.2.4 Sums and Mixtures . . . . . . . . . . . . . . . . . . . 23 1.3 Multivariate Distribution . . . . . . . . . . . . . . . . . . . 29 1.3.1 Multivariate Normal Distribution . . . . . . . . . . . 29 1.3.2 Mixtures of Multivariate Normals . . . . . . . . . . . 36 1.3.3 Wishart Distribution . . . . . . . . . . . . . . . . . . 38 1.3.4 Uniform Distribution on the d−Sphere . . . . . . . . 39 1.4 Stochastic Process . . . . . . . . . . . . . . . . . . . . . . . 42 Previous Next First Last Back Forward 1

Chapter 1 Methods for Generating Random Variables 统计计算中一个基础问题就是从特定的概率分布中产生随机变量（随机数）.在 最简单的情形，从一个有限总体中产生一个随机观测，就需要一种从离散均匀 总体中产生随机观测的方法.从而，一个合适的均匀（伪）随机数产生器是从根 本上所需要的，从其他概率分布中产生随机数都依赖于均匀随机数产生器， 1.1 Generating Uniform(0,1)random number 线性同余发生器： Xn+1 aXn +c(mod)m 这里： Previous Next First Last Back Forward 1

Chapter 1 Methods for Generating Random Variables ⁄OOé•òáƒ:ØK“¥lA½V«©Ÿ•)ëÅC˛(ëÅÍ). 3 Å{¸ú/, lòákÅoN•)òáëÅ*ˇ, “Iáò´ll—˛! oN•)ëÅ*ˇê{. l , òá‹·˛!(ñ)ëÅÍ)Ï¥lä ˛§Iá, lŸ¶V«©Ÿ•)ëÅÍ—ù6u˛!ëÅÍ)Ï. 1.1 Generating Uniform(0,1) random number Ç5”{u)Ï: Xn+1 = aXn + c (mod) m ˘p: Previous Next First Last Back Forward 1

(1)m:模=231-1(32位计算机可以直接表示的最大素数). (②)a:乘数，要小心选择 (3)c:增量（可以）=0. (4)Xo:初始值（种子） ·Xn∈{0,1，.，m-1},Un=Xn/m. ·使用奇数作为种子 ·代数和群理论有助于a的选择.。 ·希望生成器的周期很大 ·不要产生多与m/1000个数. 在R中，使用帮助主题来了解，Random.seed或者RNGkind关于缺省均 匀随机数产生器的详细信息.各种不同类型的随机数产生器及其性质可以参 考一些数值计算方面的资料 Previous Next First Last Back Forward 2

(1) m: = 231 − 1(32†OéÅå±ÜL´ÅåÉÍ). (2) a: ¶Í, á%¿J (3) c: O˛(å±)=0. (4) X0: –©ä(´f). • Xn ∈ {0, 1, . . . , m − 1}, Un = Xn/m. • ¶^¤Íäè´f. • ìÍ⁄+nÿkœua¿J. • F")§Ï±œÈå. • ÿá)ıÜm/1000áÍ. 3R•, ¶^êœÃK5 ).Random.seed½ˆRNGkind'u"é˛ !ëÅÍ )Ïç[&E. à´ÿ”a.ëÅÍ)Ï9Ÿ5üå±Î ò ÍäOéê°]. Previous Next First Last Back Forward 2

1.2 Random Generators of Common Probability Distribution in R 在R中均匀伪随机数的产生器是ruif.产生伪均匀随机数方式为 runif(n} #产生0到1上的长度为n的向量 runif(n,a,b】#产生a到b上的长度为n的向量 matrix(run1f(n*m),nrow=n,ncol=m)#产生0到1上的nxm的矩阵 Example 常用的一元概率分布的概率质量函数(pmf)或者概率密度函数(pdfD,累 积分布函数(cdf,分位数函数(quantile function)以及随机数产生器函数 在R中已经集成，比如对二项分布，可以参看帮助文档 dbinom(x,size,prob,log FALSE) pbinom(q,size,prob,lower.tail TRUE,log.p FALSE) Previous Next First Last Back Forward

1.2 Random Generators of Common Probability Distribution in R 3 R •˛!ñëÅÍ)Ï¥runif. )ñ˛!ëÅÍê™è ↑Example runif{n} #)01˛›ènï˛ runif{n,a,b} #)ab˛›ènï˛ matrix(runif(n*m),nrow=n,ncol=m) #)01˛nxm› ↓Example ~^òV«©ŸV«ü˛ºÍ(pmf)½ˆV«ó›ºÍ(pdf), \ »©ŸºÍ(cdf), ©†ÍºÍ(quantile function)±9ëÅÍ )ÏºÍ 3R•Æ²8§, 'XÈë©Ÿ, å±Îwꜩ ↑Code dbinom(x, size, prob, log = FALSE) pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE) Previous Next First Last Back Forward 3

qbinom(p,size,prob,lower.tail TRUE,log.p FALSE) rbinom(n,size,prob) Code R中常见的一元分布函数 分布 cdf 随机数产生器参数 beta pbeta rbeta shapel,shape2 二项分布 pbinom rbinom size,prob X2分布 pchisq rchisq df 指数分布 pexp rexp rate F分布 pf f dfl.df2 gamma pgamma rgamma shape,rate or scale 几何分布 pgeom rgeom prob 对数正态分布 plnorm rlnorm meanlog,sdlog 负二项分布 pnbinom rnbinom size,prob 正态分布 pnorm rnorm mean,sd Poisson分布 ppois rpois lambda t分布 pt rt df 均匀分布 punif runif min,max Previous Next First Last Back Forward 4

qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE) rbinom(n, size, prob) ↓Code R•~Ñò©ŸºÍ ©Ÿ cdf ëÅÍ)Ï ÎÍ beta pbeta rbeta shape1, shape2 ë©Ÿ pbinom rbinom size, prob χ 2©Ÿ pchisq rchisq df çÍ©Ÿ pexp rexp rate F©Ÿ pf rf df1,df2 gamma pgamma rgamma shape,rate or scale A¤©Ÿ pgeom rgeom prob ÈÍ©Ÿ plnorm rlnorm meanlog,sdlog Kë©Ÿ pnbinom rnbinom size, prob ©Ÿ pnorm rnorm mean,sd Poisson©Ÿ ppois rpois lambda t©Ÿ pt rt df ˛!©Ÿ punif runif min,max Previous Next First Last Back Forward 4

使用samplet函数从一个有限离散总体中抽样： sample可以实现从有限总体中有放回或者不放回两种抽样方式进行抽样. TExample #掷硬币 sample(0:1,size=10,replace=TRUE) #字母a-z的一个置换 sample(letters) #多项分布抽样 x<-sample(1:3,size=100,replace=TRUE,prob=c(0.2,0.3,0.5)) #table(x) #x #123 #123058 Example Previous Next First Last Back Forward 5

¶^sampleºÍlòákÅl—oN•ƒ: sample å±¢ylkÅoN•kò£½ˆÿò£¸´ƒê™?1ƒ. ↑Example #ïM1 sample(0:1,size=10,replace=TRUE) #i1a-zòáòÜ sample(letters) #ıë©Ÿƒ x table(x) #x # 1 2 3 #12 30 58 ↓Example Previous Next First Last Back Forward 5

1.2.1 The Inverse Transform Method 口连续型场合 生成随机数的逆变换方法是基于以下熟知的定理 Theorem1(Probability Integral Transformation).若X为连续型随机变 量，其cd为Fx,则U=Fx'(X)~U(0,1). Proof.定义 Fx(u)=inf{z:Fx(x)=u},0<u<1. 若随机变量U~U(0,1),则对所有x∈R,有 P(Fx'(U)≤x)=P(inf{t:Fx(t)=U}≤x) P(U<Fx(x))=FU(Fx(x))=Fx(x). 因此F(U)和随机变量X同分布. ▣ Previous Next First Last Back Forward 6

1.2.1 The Inverse Transform Method ❏ ÎY.|‹ )§ëÅÍ_CÜê{¥ƒu±eŸ½n Theorem 1 (Probability Integral Transformation). eXèÎY.ëÅC ˛, Ÿcdf èFX, K U = F −1 X (X) ∼ U(0, 1). Proof. ½¬ F −1 X (u) = inf{x : FX(x) = u}, 0 < u < 1. eëÅC˛U ∼ U(0, 1), Kȧkx ∈ R, k P(F −1 X (U) ≤ x) = P(inf{t : FX(t) = U} ≤ x) = P(U ≤ FX(x)) = FU (FX(x)) = FX(x). œdF −1 X (U)⁄ëÅC˛X”©Ÿ. Previous Next First Last Back Forward 6

从而，若要产生X的一个随机观测x,可以 1.导出逆变换函数F灭'(u). 2.从均匀分布U(0.1)中产生一个随机数u,令x=Fx(u). 这种方法要求F的逆函数要容易求出。 例：使用逆变换方法产生连续型密度f(x)=3x2I(0<x<1)的随机观测 值。 此处Fx(x)=x3,(0<x<1),因此F(u)=2/3,因此 FCode n<-1000 u<-runif(n) x<-u^(1/3) hist(x,prob=TRUE,main=expression(f(x)==3*x*2)) y<-seq(0,1,0.01) Previous Next First Last Back Forward 7

l , eá)XòáëÅ*ˇx, å± 1. —_CÜºÍ F −1 X (u). 2. l˛!©ŸU(0, 1)•)òáëÅÍu, -x = F −1 X (u). ˘´ê{á¶F_ºÍáN¥¶—. ~: ¶^_CÜê{)ÎY.ó›f(x) = 3x 2 I(0 < x < 1)ëÅ*ˇ ä. d?FX(x) = x 3 , (0 < x < 1), œd F −1 X (u) = u 1/3 , œd ↑Code n<-1000 u<-runif(n) x<-u^(1/3) hist(x,prob=TRUE,main=expression(f(x)==3*x^2)) y<-seq(0,1,0.01) Previous Next First Last Back Forward 7

11nes(y,3*y^2) Code 绘图中数学符号的表示更多内容参看帮助文档？plotmath. 例：使用逆变换方法产生指数分布的随机数. X~Exp(A),因此对x>0,Fx(x)=1-e-A江，从而Fx(u)= 一log(1-u).注意到U和1一U同分布，因此产生参数是入的指数分布的长度 为n的随机数命令为 Code -log(runif(n))/1mbda ⊥Code 在R中也可以使用rexp来产生. Previous Next First Last Back Forward 8

lines(y,3*y^2) ↓Code ±„•ÍÆŒ“L´çıSNÎwꜩ ?plotmath. ~: ¶^_CÜê{)çÍ©ŸëÅÍ. X ∼ Exp(λ), œdÈx > 0, FX(x) = 1 − e−λx, l F −1 X (u) = − 1 λ log(1 − u). 5øU⁄1 − U”©Ÿ, œd)ÎÍ¥λçÍ©Ÿ› ènëÅÍ·-è ↑Code -log(runif(n))/lmbda ↓Code 3R•è屶^rexp5). Previous Next First Last Back Forward 8