统计 Monte Carlo,Jackknife,and Bootstrap 蒙特卡罗法等.pdf

Monte CarlOjJackknife,and Bootstrap第9组程艳超，董笠，李月飞，赵小丹蒙特卡罗法Monte Carlo董笠蒙特卡洛法基本原理 实际应用优点和缺点基本原理蒙特卡洛方法又称统计模拟法、随机抽样技术，是一种随机模拟方法，以概率和统 计理论方法为基础的一种计算方法，是使用 随机数来解决很多计算问题的方法。将所求 解的问题同一定的概率模型相联系，用电子 计算机实现统计模拟或抽样，以获得问题的 近似解。为象征性地表明这一方法的概率统 计特征，故借用赌城蒙特卡罗命名。基本思想当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时,通过某种“实验”的方法，以这种事件出 现的频率估计这一随机事件的概率，或者 得到这个随机变量的某些数字特征，并将 其作为问题的解。蒙特卡罗法可以用来求解两类问题。第一类问题称之为概率问题，用直接 模拟某种过程的方法解决。第二类问题，是所谓定数问题。在解 定数问题时，必须把问题化为相适应 的能作模拟的概率问题实际应用例一随机化检验通过观察结果多次随机化并重新计算给 定统计量，从而产生所有可能结果的分布，然后利用假设检验从计算出的分布的概率直接估计。有两组数据，现求出了其相关系数X I Xl X2 X3 X4 X5 X6 X7 X8 X9 X10Y Y1 丫2 Y3 丫4 丫5 Y6 丫7 Y8 Y9 I Y10R=Rreal由于有抽样误差，需进行相关系数 的假设检验,这里我们尝试用蒙特卡 罗法。随机排列X和Y的顺序，模拟9999次。每个组合都计算R值 R值出现次数符合正态分布,R0检验Rrea I是否显著。如 Rreal P+1.96 ORreal V 口 T.96。而说明显著实际应用uffon的投针实验目的：利用蒙特卡罗法计算”的值 实验方法：画一组水平线，取长度小于水平线间距离的针，做投针实验。记录针与水平线相交的 次数。重复投针实验N次。结果：能算出兀的近似值、/设 D：两水平线间的距离L：针长(LD)a：针与相交水平线右侧的夹角X：针中点到最近水平线的垂直距离0XD/2oan针与平行线相交：ox(lsina)/2D绘出如下坐标轴(Isina)/2以角度a为横轴.生相交时每个a对应的x瓯大为P为针与水平线相交的概率Sa为阴影部分的面积Sa为矩形的面积n为投针与水平线相交的次数N为投针的总次数则有P=Pa/PA二(fJ1-sin adn)/(-tt)=21/d n Jo 2 2TT=2l/dPP=n/NTT=2IN/dn关于Buffon投针实验，历史上的相关记录是:实验者年份投针次数相交次数近似的冗沃尔夫1850500025313.1608史密斯1855320412193.1554德摩根18606003833.137福克斯188410304893.1595拉泽里尼1901342818083.1415929赖纳192525308593.1975在不用计算机的情况下，精度比较低优点：借用计算机，特别是近年来出现的高速电子计算 机，使得用数学方法在计算机上大量、高效地模 拟试验成为可能；程序结构简单，易于实现；受几何条件限制小缺点：不借用计算机，需要大量的重复实验，工作量极 大和繁琐，且不够精准；即使使用计算机，仍然比较耗时Jackknife赵小丹引言假设在一个给定地区，估计一个物种的平均繁殖 输出。统计学中的总体是关于该地区每棵植物的 数值（如每株植物种子数）。平均繁殖输出是一个 参数，表示关于感兴趣总体的一些特征。如每棵植物都被调查，该参数可被精确计算，但在现实中，完全调查一个总体是不现实的。A统计学上，对于本例，人们只调查该物种的一个 随机样本，用样本繁殖输出的平均值来估计总体均值。A该估计的准确度如何?Sample 1Estimate 1Population DistributionSample 2 Sample 3Estimate 2 Estimate 3Distribution of estimates from repeated sampling植物体都不同的繁殖输出，不 同的样本会提供不同的平均繁殖输 出估计。有的估计低于总体的均值,另外一些会高于总体的均值。从所 有可能的植物样本得到的繁殖输出 均值构成样本均值的取样分布;这 个分布的特征可以用于表述一个统 计量的准确度。统计量的准确度准确度(accuracy)r偏差(bias)：总体均值与样本均值之 差。度量一个统计量是否恒偏高或恒 偏低。L精度(precision)：取决于样本分布的变异性，常用方差或标准误来度量。如果样本分布的中心与总体均值一致，这个统计 量就是无偏的；如果样本分布围绕的中心是一个与总体均值不同 的数值，一个精确的统计量可能有偏差。如果现在我们只有一个样本数据和一个估计值，那 又该怎么估计偏差或精度呢?重抽样(Resampling)What is a resampling?传统的参数检验是比较已知统计量与其理论 分布参数之间的关系。而重抽样是在已知抽样统计量（一个样本）的基础上，对抽样样本进行这样或那样的再抽样，从而获得更多的数据（样本）以 进行检验的方法。Why use resampling?抽样成本很高，很多情况下，只能获得一个样本的数据，不可能对一个总体抽取多个样本以取得样本估计值的分布。许多生态指数比样本均值复杂得多，样本分布难于用 数学方法计算。选择具有生态意义的系数比选择统计 上可解的指数更为重要。重抽样：A.刀切法B.自助法What is a jackknife?刀切法是回答我的估计有多精确这个问题的一种通用方法。它可以用于估计统计量的偏差或标准 误，但一般不能用于对置信区间的估计。它的基 本思路是通过重复计算这些数据的多个子集的统 计量对偏差和标准误进行估计。该法每次从原样本中剔除一个样品，得样本含量为 的新样本,称为Jackknife样本洪有n个,由每个样本计算估计值,称为Jackknife估计。A设x1,Kn是抽自某一总体X的随机样本，又设。是X的分布函数F（x；e）的一个参数，它的估计量为A根据0,我们构造一个jackknife样本估计量:1）去掉原样本中的第i个数据，用剩下的大小为n1的jackknife样本来求伊，并用记为此时的估计量 供有n个jackknife估计量）。2）每个扰动后的数据与原统计量结合起来，对n个数 据点的每一个计算一个伪值（pseudo value）:Pi=9*+（一1）（*-9*_,.）3）求伪值的平均值耳，得到0的jackknife估计量:1 n n n万二卒“*一丁卒*-，因此，由刀切法估计的偏差biaq为:bias.-0-pJ刀切法对于标准误SEj的估计为:SEj=E(a-p)2n-1Pi：伪值0*：原样本的参数估计量N：样本数目生态学实例：用基尼系数度量臭椿幼苗个体大小的等级臭椿协us a/ss/ma)是一个具有入侵性的引 入物种。Evans(1983)对该树种幼苗在竞争环境中 是否比在单独生长情况下更高的大小不等性进行了 研究。他做了两个实验：1)从大量种子中随机选 择6粒种在独立花盆里。2)选100粒种在一个种植 盘里。5个月后，测量每个成活植物的叶节(leaf nodes)数仁小叶数)。基尼系数(Gini coefficient)基尼系数是20世纪初意大利学者基尼（CorradoGini）根据劳伦茨曲线所定义的判断居民收入分配公平程度的指标。A基尼系数的取值范围在0T之间。基尼系数最大为“1”，表示居民之间的收入分配绝对不平均，即100%的收入被一个单位的人全部占有了;A最小等于“0”,表示居民之间 的收入分配绝对平均，即人 与人之间收入完全平等，没 有任何差异。生态指数一基尼系数G基尼系数G是关于植物大小不均衡性的量度。它的数值从0（即所有植物大小一样）到理论上的上限1（即一个植物非常巨大而所有其他植物非常小）。计算公式：n 一 1）为G=-n（几-DE%i=l其中n是植物的个数，为是第i个植物的大小,按从小到大的顺序排列，x1 x2 xnTable 14.1 Number of leaf nodes for 5-month-old Ailanthus altissima grown under two conditions:in individual pots and in a common flat6 plants grown individually:18 18 20 23 25 28Gl=0.11275 surviving plants grown together in a common flat:2627500 12 3 4 52 2 2 2 2 24 5 6 7 8 91112131415131415131415131315271718192021222324252718192022232518202223202323仅有一个值,无法进行精 度估计G2=0.115A对于种在独立花盆的6株臭椿，其原数据为:6 plants grown individually:181820232528用jackknife处理基尼系数:Table 142 Jackknife samples with Gini coefficients for individually grown AilanthusPseudovalueJackknife sample第一个数据被除去coefficient0.11281820232528|（s）z0.12418202325280.1100.12418182325280.1200.07018182025280.1240.05318182023280.1170.08918182023250.0910.216mean p:其中伪值(psuedo-value)：Pi=G+(n-1)(G-G.j)A伪值均值:P=Pi n _A由刀切法估计的偏差为：bias=G-pA刀切法对于标准误的估计：(-0.0008)sEj=J Z（，-万Y n（in 1）n:样本数 G:样本基尼系数p：刀切法伪值的均值Jackknife的优缺点优点：运用Jackknife法对参数进行精度估计,不要求样本服从特定分布，通用性好，结果可 靠，并能对样本数据尤其是单次测量样本数据 的准确性进行检验。缺点：一般不能被用于估计置信区间或假设检 验。自助法(Bootstrap)李月飞自助法简介Bootstrap是一个很通用的工具，用来估计偏差、标准误和置信区间。它是由Bradley Efron于1979年 提出。“Bootstrap”来自短语“to pull oneself up by one,s bootstraps（源自西方神话故事“TheAdventures of Baron Munchausen吹牛大王历险 记”，吹牛大王说他掉到了湖底，没办法上来，他被救后却说是自己拎着鞋带将自己提起来的）。意义:不靠外界力量,而靠提升自己的性能,达到某种目标,翻译为自助/自举。自助法基本思想Bootstr叩方法是一类非参数Monte Car Io方法，通过再抽样对总体分布进行估 计。再抽样：将观测到的样本视为一个有限 总体，从中进行随机再抽样来估计总体特 征以及对抽样总体作出统计推断。当目标总体分布不能确定，样本是唯一已有的信息时，Bootstrap是经常使用的方法。自助法步骤1）采用再抽样方法，从观测样本中抽取容量相同的样本（自助样本），此过程允许重复抽样。2）根据抽出的自助样本计算给定的统计量（如:G）3）重复上述步骤N次得到N个自助样本和N个统计量（如：Gi,i=1 j2v.,N）o4）通过N个样本统计量的分布计算统计量的偏差、标 准误和置信区间。自助法估计偏差和标准误偏差：自助法分布的平均值与原样本平均值之差。将N个自助样本的均数进行平均便得自助样本均数。令Gi为第i个自助样本均数,i=1,2,N,则自助样本均数,NG*=l/NfGii=l故:bias=G G*标准误:自助法估计标准误SE=Z(GlG*)21N 1在实际应用中，需要50100个自助样本来估计 标准误。自助法的三个问题1.如何趋近总体2.如何抽取自助法样本3.如何计算置信区间如何趋近总体参数自助法和非参数自助法参数自助法：通过一个特定分布去趋近总体值（如：对数正态分布或泊松分布），并通过样本去估计分布的参数。（对于非常小的样本此法更好）非参数自助法：通过非连续分布的观测值来趋近，即假定观测样本为总体，对其进行重复抽样。（生态学中大多应用此方法,它要求的关于总体的假定条件比较少）参数bootstrap vs.非参数的bootstrap取决于分布函数是否已知参数自助法利用分布函数，需要估计函数模 型中的参数。当函数模型正确时，参数自助法 能得到较精确的结果；非参数自助法不利用函 数就能得到较精确的结果。使用参数自助法与使用普通参数方法其效果 相当，并且操作简单。但实际应用中非参化自 助法更具代表性。如何抽取自助法样本一般的自助法:可替代地从，总体即已知样本代替）数值中随 机抽取样本。平衡的自助法:强迫，总体内每个数值具有相同被抽取的频率。移动区组的自助法:产生的自助法样本能在一定程度上保持其原始数据的相关性，是一种非参的方法。如何用自助法估计置信区间1.百分位数自助法2.偏差校正自助法3.基本自助法4.t检验化的自助法百分位数自助法本方法中，自助法分布的2.5和97.5的百分位数被用来估计95%置信区间的误差限。方法：把N个自助法样本中计算的估计值从小到大排序，如果取0.05的置信水平，那么从抽样分布的两尾处各数0Q25N个观察值，便可得到百分位数置信限。百分位数自助法例如：N=1000,分别取第0.025N个和(L0.025N)个数值,即第25和第975个数据，右图为0.133和0.171,故95%的 置信区间为(0.133,0.171)oO竞争条件下臭椿个体基尼系数自助样本频率直方图偏差校正自助法偏差校正的自助法对不以观测统计量为中心的自助法分 布进行调整。其置信区间的上下限是通过确定L来找到的。利用正态分布累积函数校正，对于95%的置信区间，则：下限PF（2Zo-1.96）上限P产 O（2Z0+1.96）定义L为自助重复中小于观测值的比率，查正态分布累积 函数（x）表，Z。为L对应的概率值转换，L 96对应的是 0.025和0.975的概率值转换。得到修正后的匕，%再按百分数自助法找到置信区间上 下限。同样用百分位数法中的例子，56.4%的自助值小于观测值，因此L=0.564,查0(x)表:Z0=0.166,下限百分数P=P(1628)=5.181%,上限百分数P产 0(2.29)=98.91%,即第52个和989个自助法值，查得为(0J36,0.171)o基本自助法利用观测值G=x的假设检验来寻找上下限值x=l和x=u。若原始样本观测值G正好落在平移函数（均值为u）第2,5个百分位数时，G=x在p=0.05下刚好被拒绝；同理原始样本观测值G正好落在平移函数（均值为I）第97.5个百分位数时，G=x在p=0,05下刚好被拒绝。故基本自助法的95%置信区间为:(2G Gq,9755 2G Gq.025)同样以上述为例，第2.5个和97.5个自助法分布的百分位数分布为0.133和0.171,95%的基本自助法置信区间为:(2X 0.155-0.171,2X0.155-0.133),=(0.139,0.177)ot-检验化的自助法当一个估计量g具有正态抽样分布，置信区间 可以表示为：(g-tSg,g+tSg)其中，Sg是g的标准误，t为适当自由度的t分布分位数。则在t检验化自助法中，计算区间方法如下:(g+b|Sg,g+bhSg)其中，Sg是g的标准误，b1和如分别为2.5和 975的百分位数。这个自助法分布，通过从观测值中抽取自助法样本并计算其统计量Gj,和其标准误生 t检验化统计量为：bj=(GrG)/Sj(G=0.155为原样本观测值。)同样以上述为例，即第25个和第975个bj值。对应的值分别为T.80和2.48o观测值G=0.155,标准误为0.0102,故95%置信区间为:(0.155-0.180X0.0102,0.155+2.48X0.0102)=(0.137,0.180)四种自助法的置信区间比较置信区间估计方法总结t检验化自助法脸统 造化量 构验计-本法 基助假验 造检一 构设联校助偏正法位助分自百数法了积布，偏 用态分数正 利正累函校差利助分征接自本特直用样布以上方法不能说谁好谁坏,在实际应用中，可比较 不同方法计算结果差异是否明显。R-project程艳超R简介 R是一个开源（GNU）统计软件，数学计 环境 S语言的一个分支，解释性语言 程序包可以从镜像下载安装:Cran=Comprehensive R Archive的简称（R综合典藏网）Network有UNIX、LINUX、MacOS和WINDOWS版 本，都是可以免费下载和使用的 http:www.rproiect.orgW R简介 k是一个开源（GNU）统计软件，数学计算环境 GNU：GNU计戈ij,又称革奴计划，是由RichardStallman在1983年9月27日发起。其目标是创建一套完全自由的操作系统，“重现当年软件界合作互助的团结精神”。为保证GNU软件可以自由地“使用、复制、修改和发布，所有GNU软件都存在一份在禁止其他人添加任何限制的情况下授权所有权利给任何人 的协议条款，GNU通用公共许可证（GNU General Public License,GPL）o 即“反版权”（或称Copyleft）概念。R简介提供一些集成的统计工具，但更重要的是 它提供各种数学计算、统计计算的函数，从而使使用者能灵活机动的进行数据分析,甚至创造出符合需要的新的统计计算方法。每个人都可以从网上下载到源代码，并可以对它进行修改以符合自己的使用要求。R简介胰 RGui文件编辑查看其他程序包窗口 帮助叵回 国回回 S屋 R Console 亘Tll-IR version 2.14.2(2012-02-29)Copyright(C)2012 The R Foundation for Statistical Confuting ISBN 3-900051-07-0Platform:i386-pc-mingw32/i386(32-bit)R是自由软件，不带任何担保。在某些条件下你可以将其自由散布。用,license()，或licence()，来看散布的详细条件。R是个合作计划，有许多人为之做出了贡献.用,contributors(),来看合作者的详细情况用，ciiaSon()，会告诉你如何在出版物中正确地引用R或R程序包。用，demo()，来看一些示范程序，用，help()，来阅读在线帮助文件，或 用，help.stare()，通过HTML浏览器来看帮助文件。用，q()，退出R.界面1、2、3、4、R基础R对大小写敏感：“A”与“a”是不同的 函数间常常嵌套：如 rm(list=ls()基本算术操作：+、*、/、八乘方允许在“语言上计算”(computing on the language)。R的一些简单函数sum(x)length(x)mean(x)sd(x)var(x)sqrt(x)sort(x)#求和#“长度”#求均值#求标准差#方差#“根号”#排序数据录入X Gini-function(x)+y-sort(x)+n-length(y)+if(n2)return(NA)+i X1 18 18 20 23 25 28 Gini-function(x)+y-sort(x)+n-length(y)+if(n2)return(NA)+i-l:n+Gjackpseudofor(i in 1:length(x)+for(j in 1:length(x)+if(ji)jackji)jackj-1-xj+pseudoi jack pseudo for(i in 1:length(x)+for(j in 1:length(x)+if(ji)jaclcji)jackj-l-xj+pseudoi pseudo1 0.12448166 0.12448166 0.07004870 0.05346122 0.08861512 0.215996501Jackknife在R中的实现#先简单看一下jackknife结果pseudo#看看得到的是什么hist(pseudo)mean(pseudo)var(pseudo)#柱形图#求均值#求方差R R Graphics:Device 2(ACTIVE)Histogram of pseudo0C997 ONInL QL 900.0 A o un bu_pseudoJackknife在R中的实现biassebootfor(i in 1:1000)+booti hi st(boot)#柱形图biasbootcorrectedquantile(boot50.975)quantile(boot,0.025)Bootstrap在R中的实现#置信区间的bootstrap实现 偏差校正的置信区间估计Gini(x)-bias-1.96*sqrt(var(boot)Gini(x)-bias+1.96*sqrt(var(boot)Bootstrap在R中的实现#置信区间的bootstrap实现 基本自助法置信区间估计2*Gini(x)-quantile(boot50.025)2*Gini(x)-quantile(boot50.975)作业计算变异系数CV,bias,se8.26 6.33 10.4 5.27 5.355.61 6.12 6.19 5.2CV函数：CV-function(x)var(x)/mean(x)提示：录入X向量，将Jackkn ife和Bootstrap 命令行中的G i n i函数置换成CV(X)即可。也可以“手动jackknife”或用其他统计软件。