班级:统计1702学号:17271119姓名:成长锦实用统计软件——R语言R语言随机数的产生与应用一、一般随机数的产生、分布列分布函数产生随机数二、圆域、三角形区域均匀分布随机数三、逆变换、合成法、剔除法一、一般随机数的产生、分布列分布函数产生随机数■R语言具有产生不同分布的随机数的函数,包括概率中常用的分布:正太、泊松、二项……,以下做一个简单归纳:分布随机数函数分布函数分位数函数概率函数正态分布rnorm(n,mean,sd)dnormpnormqnorm指数分布rexp(n,rate)dexppexpqexp均匀分布runif(n,min,max)dunifpunifqunifF分布rf(n,df1,df2)dfpfqft分布rt(n,df)dtptqt𝜒2分布rchisq(n,df)dchisqpchisqqchisqBeta分布rbeta(n,shape1,shape2)dbetapbetaqBetagamma分布rgamma(n,shape,scale=1)dgammapgammaqgamma二项分布rbinom(n,size,prob)dbinompbinomqbinom泊松分布rpois(n,lambda)dpoisppoisqpois几何分布rgeom(n,prob)dgeompgeomqgeom■R输出举例:1.均匀分布2.正太分别3.指数分布4.二项分布随机数产生并绘图。■set.seed()函数与随机数产生●set.seed()只对运行该命令后的第一次随机产生结果有效,设置相同的种子数产生的随机数一致,一些实际概率编程求解问题能使用到set.seed()函数,比如大数定律的模拟。从下图能很好理解该函数的意义。■sample()抽样函数与一般随机数产生●使用sample()抽样函数能够从某些数据集中等可能的产生随机数或指定概率或比例产生随机数数据集中等可能的产生有用的随机数据指定概率或比例产生随机数■举例:设随机变量X的密度函数为推导:产生的随机数用直方图和密度函数图展示■根据分布函数产生随机数,若分布函数严格单调增则有,通过分布函数反解,得到,即为该分布函数下的随机数。𝐹(𝑥)~𝑈(0,1)()xirxU𝑓(𝑥)=100𝑥2,𝑥1000 ,𝑥≤1002100100100()1~(0,1)1001xFxUxxXU(逆变换法)需要注意有:1.密度函数的分布函数严格单调增2.一般能很好反解出x关于U的表达式■根据分布列产生随机数,此时需要用到sample()函数𝑿 −𝟐 −𝟏 𝟎 𝟏 𝟐 𝟒 𝑷 𝟎.𝟐 𝟎.𝟏 𝟎.𝟑 𝟎.𝟏 𝟎.𝟐 𝟎.𝟏 ■举例:设随机变量X的分布列为产生的随机数局部展示和表的形式展示如下:二、圆域、三角形区域均匀分布随机数■例7(3)(P88)产生圆内的随机数二维随机变量(X,Y)服从圆域的均匀分布,试求p(X1/2|Y=1/2)221xy方法一:剔除法该圆域随机数如右上图所示:剔除法产生随机数比较简单,但其效率往往不高,对实际问题求解需要产生更多随机数才能得到更多有效随机数。效果如右下图所示本题的随机数10000个只有7788个是可用的随机数,效率为0.7788。服从单位圆的均匀分布随机数点可以通过如下产生:■例7(3)(P88)产生圆内的随机数二维随机变量(X,Y)服从圆域的均匀分布,试求p(X1/2|Y=1/2)221xy方法二:公式法该圆域随机数如右图所示:公式法产生随机数高效,要多少有多少,十分精确,对实际问题求解产生的随机数都是有效随机数。本题的随机数10000个全有效,效率为1。cos,sinxryr■产生圆内的随机数扩展221xy服从单位圆的均匀分布随机数点可以通过如下产生:该一般的单位圆的极坐标公式是不带根号的,但为什么圆类随机数的产生需要对r开根号呢?我们试着不开根号的极坐标公式产生随机数结果由上图所示:发现随机数离原点越近越密集。由于圆的圆周与距离r成正比,因此r的概率密度函数也应与r成比例:这里r=1,故有因此,应该遵循均匀分布的r的平方。故产生单位圆的随机数公式需要对r开平方cos,sinxryr()fxr2()Fxr2()Fxr222()()()~prRprRFRRURU221xy剔除法进一步改进,矩阵思想简化代码■产生圆内的随机数扩展产生的随机数如右图所示:用矩阵直接可以讲由均匀分布产生的x,y直接放入一个10000*2的矩阵,编写f函数作为判断准则,再用apply函数对矩阵直接操作即可。方法一:剔除法■例8(P93)产生三角形内的随机数随机变量(X,Y)的概率密度函数为1,||,01(,)0,yxxfxyotherwise产生的随机数如右图所示,剔除法剔除规则在这道题下为不满足|y|x的剔除,但效率不高,仅仅约0.5剔除法进一步改进:指定产生有效随机数个数■例8(P93)产生三角形内的随机数随机变量(X,Y)的概率密度函数为1,||,01(,)0,yxxfxyotherwise采用while语句可指定生成有效随机数的个数,但注意whlie的判断条件与剔除规则恰好相反(即不是有效随机数则重新产生随机数直至其有效),产生的随机数最终如右图所示。这种改进使得剔除法更加灵活,产生的最终随机数即为指定的n,该题中即为1000,但牺牲了算法效率设三角形顶点A,B,C,在该三角形区域的均匀分布随机数点,可以通过如下产生:■例8(P93)产生三角形内的随机数随机变量(X,Y)的概率密度函数为1,||,01(,)0,yxxfxyotherwise方法二:公式法1122112(1)[(1)](),~(0,1)PrArrBrrCrrU产生的随机数如右图所示:采用该方法效率很高,产生的随机数都是有效的,效率达到1■例8(P93)产生三角形内的随机数随机变量(X,Y)的概率密度函数为1,||,01(,)0,yxxfxyotherwise方法三:合成法合成法构思:•如果X的密度函数难于抽样,而X关于Y的条件密度函数以及Y的密度函数均易于抽样,则X的随机数可如下产生:•可以证明由此得到X的服从。本题中可以求得x的边际密度函数及边际分布函数如下,通过逆变换法得到X的随机数;并发现Y关于X的条件密度函数形如均匀分布的密度函数,故可以通过条件密度函数对应的分布函数产生Y的随机数。𝑓𝑋(𝑥)=2𝑥⇒𝐹(𝑥)=𝑥2𝑓𝑌|𝑋(𝑦|𝑥)=𝑓(𝑥,𝑦)𝑓𝑋(𝑥)=12𝑥■例7(P115)设(x,y)服从区域上的均匀分布,设区域B为产生该区域D均匀分布的随机数方法一:剔除法2,|01Dxyyx𝑥,𝑦|𝑦≥𝑥2从D区域的图像可以知道X,Y的范围在[-1,1]*[0,1]的矩形之间,故可以通过产生该区域随机数再通过判断准则剔除。随机数如下图所示:进一步还可以画出B∩D区域的随机数如右下图:■例7(P115)设(x,y)服从区域上的均匀分布,设区域B为产生该区域D均匀分布的随机数方法二:合成法2,|01Dxyyx𝑥,𝑦|𝑦≥𝑥2先求f(x,y)联合密度,在去看边际分布函数与条件密度函数,从一下推导后进行编程产生随机数发现x与U的解不好求,故用polyroot求解,因为是一元三次方程存在三个解,只能取其中在-1~1之间的一个解。随机数如右下图:232|23(,)3/4;()(,)(1)431()()~(0,1)441(|)|~(0,1)1XXXYXfxyfxfxyxFxfxxxUfyxYXUxx■例7(P115)设(x,y)服从区域上的均匀分布,设区域B为产生该区域D均匀分布的随机数2,|01Dxyyx𝑥,𝑦|𝑦≥𝑥2进一步解释为什么可以保证-1~1之间的解有且只有一个,因为F(x)=U的解的范围从右图中可以看出在[-1,1]中只可能存在一个,在[-1,1]外一定存在两个。故求解的方程解有这一取值:a-Re(polyroot(c(4*runif(1)-2,-3,0,1)))X[i]-a[a1&a-1]■可进一步改进本程序的代码如下:实际只需两行代码搞定随机数效果如右下图所示方法二:合成法