Яࠔእ̮వ˺௧ʷవ3ឦᝓ֗ѬౢᄊК᫃ஔెὊ॰ऀຒᤉnjງКูѣὊඈ˔ᅼគགࡊ᧚̰ࠄᬅᄊऄၹವΓѣԧὊ̿᫈ᮥ˞ՔὊښᝍх᫈ᮥ˗ߦ˸ፒᝠவขnj3ឦᝓᄊ۳వΎၹ̿ԣᎄሮࢼǍవ˺Яࠔᄦ3ፇnjѦˁ͖ӑnjચನവલnjፒᝠѬౢnjϜೝᰎnjڀॆѬౢnjፒᝠፋڏ֗3ӊ҄ͻЯࠔǍవ˺ᄊࠀͯ௧˞ˊႍѬౢ̡րnjፃเኮေዝnjӞߦᄊߦၷଢΙவข֗ሮऀʽᄊԠᏦὊښиͻሮ˗ࡊ᧚ѻԝඋᣗေᄊߦԔေὊᤈನᑟܵࣟҰ឴Ꮸᣐʽߦ˸ǍళፃԻἻˀ४̿͊͵வरܭ҄ੋᜃవ˺˨ᦊѬੋЛᦊЯࠔǍྠిਫ਼దἻΦిॹቃǍڏ˺ښྠᎄᄬἷ$*1Ἰ3ѬౢவขˁವΓហᝍவӗӯὊసथࣱὊކՁᮻᎄᗃúӒ̛ႃߕࢺˊѣྠᇫὊ*4#/Ĉŀ3ĀĉŀவĀŁసĀłކĀĊŀሮऀឦᝓὋሮऀᝠὋஔెċŀ51˗ڎྠవڏ˺ᯞ$*1ߚኄՂ᠊͊ᎄᣤसథᖱӿNjNj҅ʼෲ࣊ԥӿ҅ᜉᝡదᬍНՃᜉNjNjᝡʼෲ࣊ԥӿ҅ᜉᝡదᬍНՃѣྠԧᛡႃߕࢺˊѣྠᇫӒ̛࣊๒ӝʺηኸNjᥫᎄनNjNjవfNjNjӿसߚӢߚྠNjNjࣲథኄྠӿNjNjࣲథኄӿ҅ࠀNjNj͉Ћјਫ਼᠔˼ႃߕࢺˊѣྠᇫڏ˺దᎥ૯᫈ᮥὊឰՔ᠔˼˺इូ૱Ǎᔪ˺इᎥὊឰˁవᇫԧᛡᦊᐏጇὊᐏጇԣᥫ᠔ႃភǍ᠏᧚આងឰԧᥫ͈ᒰ[MUT!QIFJDPNDOὊᄧྠΦిˡઑឰԧᥫ͈ᒰECRR!QIFJDPNDOǍҬབྷጳǍ第4章随机数与抽样模拟我们在做数据分析时,往往需要用到随机数与抽样模拟,这些在R语言环境下都能够轻松实现,因为R内部自带了大量的相关工具和方法。4.1 一元随机数的产生4.1.1 均匀分布随机数均匀分布随机数是最简单的随机数,也是其他分布随机数的基础。R语言生成均匀分布随机数的函数是runif(),其句法是:runif(n,min=0,max=1)n表示生成的随机数数量,min表示均匀分布的下限,max表示均匀分布的上限,若省略参数min、max,则默认生成[0,1]上的均匀分布随机数。例如:runif(3,1,3)#生成3个[1,3]的均匀分布的随机数[1]1.2041.3592.653runif(5)#默认生成5个[0,1]上的均匀分布的随机数[1]0.27840.77550.41070.83920.7455书籍2.indb542015/1/2814:32:42第4章 随机数与抽样模拟 55R提供了多种随机数生成器(randomnumbergenerators,RNG),默认采用Mersennetwister方法产生随机数,该方法是由MakotoMatsumoto和TakujiNishimura于1997年提出的一种随机数生成器,其循环周期是1993721−。R里面还提供了Wichmann-Hill、Marsaglia-Multicarry、Su-per-Duper、Knuth-TAOCP-2002、Knuth-TAOCP和L'Ecuyer-CMRG等几种随机数生成方法,用RNGkind()函数更改。例如要改为Wichmann-Hill方法:RNGkind(kind=Wich)runif()默认每次生成的随机数是不一样的,有时我们在做模拟时,为了比较不同的方法,就要求生成的随机数都一样,即重复生成同样的随机数,这时候可以使用set.seed()设定随机数种子,其参数为整数。set.seed(1)#种子取一样,生成的随机数相同runif(5)[1]0.26550870.37212390.57285340.90820780.2016819接下来,检验一下runif()生成的随机数的性质。通过直方图和散点图以及自相关系数图来检验独立同分布。从图4-1可以看出,基本上满足独立同分布性质。Nsim=10^3x=runif(Nsim)x1=x[-Nsim]#因为要求自相关系数,去掉最后一个数x2=x[-1]#去掉第一个数par(mfrow=c(1,3))hist(x,prob=T,col=gray(0.3),main=uniformon[0,1])#直方图curve(dunif(x,0,1),add=T,col=red)#添加均匀分布密度函数线plot(x1,x2,col=red)acf(x)#画自相关系数图图4-1 模拟均匀分布曲线图书籍2.indb552015/1/2814:32:4256 R数据分析——方法与案例详解4.1.2 正态分布随机数正态分布是古典统计学的核心,它涉及两个参数:位置参数均值μ和尺度参数标准差σ。正态分布的图形如倒立的钟型,对称分布。现实生活中,很多分布是服从正态分布的,例如,人类的智商IQ,大致服从均值为100,标准差为16的正态分布。所有的正态分布都可以通过标准化成均值为0、标准差为1的标准正态分布。假设1U和2U是[0,1]上均匀分布随机数,做如下变换。1122log()cos(2)XUU=−p,2122log()sin(2)XUU=−p这样得到的随机数独立同分布于N(0,1)N。正态分布随机数的生成函数是rnorm(),其句法是:rnorm(n,mean=0,sd=1),其中,n表示生成的随机数数量;mean是正态分布的均值,默认为0;sd是正态分布的标准差,默认时为1。正态分布随机数R提供了Kinderman-Ramage、Ahrens-Dieter、Box-Muller和逆变换法。默认的是逆变换法。例如:rnorm(5,10,5)#产生5个均值为10,标准差为5的正态分布随机数[1]3.17214.7057.1735.8428.879rnorm(5)#默认5个生成标准正态分布随机数[1]-0.582040.046060.96016-0.68698-0.35504下面随机产生100个正态分布随机数,作它们的概率直方图,然后再添加正态分布的密度函数线。程序如下:x=rnorm(100)hist(x,prob=T,main=normalmu=0,sigma=1)#作直方图curve(dnorm(x),add=T)#在直方图上添加标准正态分布密度函数线最终得到的模拟正态分布曲线图如图4-2所示。图4-2 模拟正态分布曲线图书籍2.indb562015/1/2814:32:42第4章 随机数与抽样模拟 574.1.3 指数分布随机数指数分布是统计分布中非常重要的一个分布,它有着诸多应用。例如,可以用来描述电子产品的寿命,一个灯泡的平均寿命是2500小时,我们可以认为这个灯泡的寿命服从均值为2500的指数分布。如果一个变量x服从指数分布,记为x~exp(λ),其中λ等于均值的倒数。R生成指数分布随机数的函数是rexp(),其句法是:rexp(n,lamda=1),n表示生成的随机数个数,lamda=1/mean。例如:x=rexp(100,1/10)#生成100个均值为10的指数分布随机数hist(x,prob=T,col=gray(0.9),main=均值为10的指数分布随机数curve(dexp(x,1/10),add=T)#添加指数分布密度函数线得到的模拟指数分布曲线图如图4-3所示。图4-3 模拟指数分布曲线图例4.1 rexp()和逆变换法比较。前面我们介绍了可以通过逆变换方法得到指数分布随机数,下面比较一下rexp()和逆变换法生成~(1)XExp的随机数,程序如下:Nsim=10^4U=runif(Nsim)X=-log(U)Y=rexp(Nsim)par(mfrow=c(1,2))hist(X,freq=F,main=ExpfromUniform)curve(dexp(x,1),add=T,col=red)hist(Y,freq=F,main=ExpfromR)curve(dexp(x,1),add=T,col=red)得到的指数分布随机数图如图4-4所示。书籍2.indb572015/1/2814:32:4258 R数据分析——方法与案例详解ExpfromUniformXDensity02468100.00.10.20.30.40.50.6ExpfromRYDensity02468120.00.10.20.30.40.50.6图4-4 rexp()生成随机数和逆变换法的比较从图4-4可以看出,这两者生成的随机数和指数分布~(1)XExp都很接近。4.1.4 离散分布随机数的生成设随机变量X的分布列{}iiPXxp==,1,2,i=。记(0)(0)0pPX=≤=,()1(),1,2,,iiijjpPXxpi==≤==∑()1(),1,2,,iiijjpPXxpi==≤==∑,设r是[0,1]区间上的均匀分布随机数。当且仅当(1)()iiprp−时,令iXx=,则(1)()()(1){}{}iiiiiiPprpPXxppp−−===−=,1,2,i=二项分布随机数二项分布是指n次独立重复贝努力试验(Bernoullitrials)成功次数的分布,每次贝努力试验的结果只有两个,成功和失败,记成功的概率为p。如果一个变量x服从二项分布,记为x~B(n,p),n表示试验次数,p表示成功概率。R生成二项分布随机数的函数是rbinom(),其句法是:rbinom(n,size,prob),n表示生成的随机数数量,size表示进行贝努力试验的次数,prob表示一次贝努力试验成功的概率。首先,生成二点分布(一次贝努力试验)的随机数。size=1;p=0.5rbinom(10,size,p)[1]000100110ExpfromUniformExpfromR书籍2.indb582015/1/2814:32:43第4章 随机数与抽样模拟 59接下来,生成服从B(10,0.5)的二项分布随机数。size=10;p=0.5rbinom(5,size,p)#生成5个服从B(10,0.5)的二项分布随机数[1]56333二项分布是离散分布,但随着试验次数n的增加,二项分布越来越接近于正态分布程序。下面将分别产生100个n为10、15、50,概率p为0.25的二项分布随机数:par(mfrow=c(1,3))p=0.25for(ninc(10,20,50)){x=rbinom(100,n,p)hist(x,prob=T,main=paste(n=,n))xvals=0:npoints(xvals,dbinom(xvals,n,p),type=h,lwd=3)}par(mfrow=c(1,1))结果如图4-5所示。从图中可以看出,随着实验次数n的增大,二项分布越来越接近于正态分布。图4-5 模拟二项分布曲线图4.1.5 常见分布函数表除了生成前面介绍的几种分布的随机数外,还可以生成poisson分布、t分布和F分布等多种分布的随机数,只要在相应的分布名前加“r”就可以,在此不一一赘述。现把常见分布归纳成表4-1,供读者参考。书籍2.indb592015/1/2814:32:4360 R数据分析——方法与案例详解表4-1 常见分布函数表分 布中文名称R中的表达参 数Beta贝塔分布beta(a,b)shape1、shape2Binomial二项分布binom(n,p)size、probCauchy柯西分布cauchy()location、scaleChi-square卡方分布chisq(df)dfExponential指数分布exp(lamda)rateFF分布f(df1,df2)df1、df2Gamma伽玛分布gamma()shape、rateGeometric几何分布geom()probHypergeometric超几何分布hyper()m、n、kLogistic逻辑