本文讲述一下蒙特卡洛模拟方法与定积分计算，首先从一个题目开始：设，用蒙特卡洛模拟法求定积分的值。 随机投点法 设服从正方形 上的均匀分布，则可知 分别服从上的均匀分布，且相互独立。记事件 ，则的概率为 即定积分 的值 就是事件出现的频率。同时，由伯努利大数定律，我们可以用重复试验中出现的频率作为 的估计值。即将看成是正方形内的随机投点，用随机点落在区域中的频率作为定积分的近似值。这种方法就叫随机投点法，具体做法如下： 图1 随机投点法示意图 1、首先产生服从 上的均匀分布的 个随机数（ 为随机投点个数，可以取很大，如 ）并将其配对。 2、对这对数据 ，记录满足不等式的个数，这就是事件 发生的频数，由此可得事件 发生的频率 ，则 。 举一实例，譬如要计算，模拟次数时，R代码如下： n=10^4; x=runif(n); y=runif(n); f=function(x) { exp(-x^2/2)/sqrt(2*pi) } mu_n=sum(y<f(x)); J=mu_n/n; J 模拟次数 时，令,其余不变。 定积分的精确值和模拟值如下： 表1 精确值 0.3413447 0.342 0.344 0.34187 0.341539 0.341302 注：精确值用integrate(f,0,1)求得 扩展 如果你很细心，你会发现这个 方法目前只适用于积分区间 ，且积分函数 在区间上的取值也位于 内的情况。那么，对于一般区间 上的定积分 呢？一个很明显的思路，如果我们可以将 与 建立代数关系就可以了。 首先，做线性变换，令 ，此时， ,。 进一步如果在区间上有 ，令 ， 则。此时，可以得到 。 其中,,。 这说明，用随机投点法计算定积分方法具有普遍意义。 举一个实例，求定积分 。 显然，，由于在 上时单调减函数，所以 ，。R中代码为 a=2; b=5; g=function(x) { exp(-x^2/2)/sqrt(2*pi); } f=function(y) { (g(a+(b-a)*y)-c)/(d-c); } ...
阅读全文 »

用过底层语言做计算的人转入R语言的时候一般都会觉得R语言的运算太慢，这是一个常见的对R的误解或者对R的设计的不理解。在二三十年前Chambers等一群人在贝尔实验室设计S语言之前，统计计算主要是通过那些底层语言实现的，典型的如Fortran。当时有一个基于Fortran的统计分析库，用它的麻烦就在于无论做什么样的数据分析，都涉及到一大摞繁琐的底层代码，这让数据分析变得很没劲，统计学家有统计学家的事情，天天陷在计算机程序代码中也不是个办法，要摆脱底层语言，那就只能设计高层语言了。有所得必有所失，众所周知，高层语言一般来说比底层语言低效，但对用户来说更友好。举个简单的例子，用C语言计算均值时，我们得对一个向量（数组）从头循环到尾把每个值累加起来，得到总和，然后除以向量的长度，而均值在统计当中简直是再家常便饭不过了，要是每次都要来这么个循环，大家也都甭干活儿了，天天写循环好了。 前两天COS论坛上有个帖子提到“R语言的执行效率问题”，大意如下： 有120000行数据，用perl写成12万行R命令做t.test，然后执行，大概几分钟就算完了。如果只用R语言，把所有数据先读入，然后用循环，每一行做t.test，花了几个小时，到现在还没有算完。这说明一个问题，在R里执行单行命令要比用循环快，涉及到循环的问题，最好写成单行命令执行。 我不清楚作者是如何写这“12万行”R命令的，本文假设是把t.test(dat), i = 1, 2, ..., 120000写入一个文件，然后用source()执行之。面对这样一个问题，我们有若干种改进计算的可能性。首先我们看“硬”写入程序代码的方法： ## 为了使计算可重复，设定随机数种子 set.seed(123) ## 生成数据，随机数，120000行 x 100列矩阵 dat = matrix(rnorm(120000 * 100), 120000) nr = nrow(dat) nc = ncol(dat) ## 六种方法的P值向量 p1 = p2 = p3 = p4 = p5 = p6 = numeric(nr) ## via source() f = file("test.t") writeLines(sprintf("p1 = t.test(dat)$p.value", seq(nr), seq(nr)), f) system.time({ source(f) }) # user...
阅读全文 »

公元1777年，法国科学家D·布丰(D.Buffon 1707～1788)设计了一个巧夺天工的实验：往间距为a的平行线族之间投掷长为L 的针，可以计算出针和平行线相交的概率。
阅读全文 »

介结了在R中如何实现极大似然估计的一种实现方法。当然其在R中还有其它的实现方法。
阅读全文 »

即日起，统计之都网站成立“统计学博文导读”栏目，归属于“网站导读”栏目。我们号召广大读者和作者将喜爱的统计学博客文章推荐给我们，以方便更多读者在这个信息爆炸的时代能够快速阅读到优秀的文章；本文是统计之都“统计学博文导读”第一篇，权当示范本栏目的作用。这次我们重点推荐两篇博文，分别来自于刘思喆和左辰，向大家展示统计学理论的生活和思维魅力。
阅读全文 »