从中心极限定理的模拟到正态分布

昨日翻看朱世武老师的《金融计算与建模》幻灯片(来源,幻灯片“13随机模拟基础”),其中提到了中心极限定理(Central Limit Theorem,下文简称CLT)及其SAS模拟实现。由于我一直觉得我们看到的大多数对CLT的模拟都有共同的误导性,因此在此撰文讲述我的观点,希望能说清楚CLT的真实面目,让读者对它有更深刻的理解。

一、中心极限定理说了什么

从广义的角度来讲,CLT说的是一些随机变量之和(在n趋于无穷的时候)趋于正态分布,条件是这些随机变量之间相关性不要太严重,每个随机变量自身的方差相对来说不要太大。我们平时看到的都是CLT的最简单的形式:随机变量独立同分布,方差有限。实际上独立同分布不是必要的条件。要理解这段话,我们就得去思考一下Lindeberg条件或Liapounov条件的含义,以及非独立情况下的CLT。而CLT为什么偏偏选中了正态分布呢?这主要是泰勒展开的功劳(注1),加上特征函数作为“催化剂”,捣鼓半天,最后一看,特征函数的极限形式正是正态分布。本文不打算走到这样偏僻的数学角落,只是先介绍一下CLT的前世今生,当我们在统计理论上犯迷糊的时候,回溯到事情的最原始状态,往往能让头脑更清醒(注2)。

即使是最简单的CLT,我觉得几乎所有做模拟的人都走入了一个误区,就是把钟形密度曲线解释为正态分布。给一个非正态的总体分布,重复计算若干次样本均值,把密度以直方图或密度曲线的形式画出来,给观众们说,看:钟形曲线!群众们一看,哇,原本是左偏/右偏/水平的密度,现在真的变成了对称的钟形曲线。于是乎不得不信服CLT的魔力。

二、中心极限定理的模拟

我们缺的不是钟形曲线,而是样本均值的分布和正态分布差多少的指标。这一点在下面第3小节中再谈,先看几个常见的模拟例子。

1、对称的钟形曲线不代表正态分布

很多人喜欢用掷骰子的例子来讲中心极限定理,大意是将两次独立掷骰子的结果加起来,看这个和的分布。我犹记当年人大一位老师给我们上抽样就举了这个例子说明CLT的魔力:看,即使样本量为2,得到的分布也是正态分布!恰好昨天在朱老师的幻灯片中又看到这个例子,不禁感叹这样一个糟糕的例子竟然经久不衰。这个例子的骗局在哪儿呢?其实很简单:在于掷骰子得到的结果的“对称性”,换句话说,结果是1、2、3、4、5、6,这6个数字围绕其均值3.5左右对称,因此两次骰子的结果加起来也围绕7对称,再画个图,两边低,中间高,看似正态分布好像出来了。朱老师的SAS代码如下:

data a;
do x1=1 to 6;
do x2=1 to 6;
output;
end;
end;  /*模拟掷骰子两次,生成36行数据*/
data a;
set a;
x=sum(x1,x2);
proc univariate data=a noprint;
var x;
histogram/normal (mu=est sigma=est);
run;

先跑题说一处细节:这段代码中的直方图不合适,因为我们知道结果x只可能是整数,而(SAS)直方图默认的分组也会以整数为边界,每个边界上的整数都被归到右边的条中去了,而这个所谓的直方图其实只是个条形图:展示每种结果的频数。所以这种情况下认为指定非整数的分组边界才能显示数据真实的分布。以下是R代码:

x = rowSums(expand.grid(1:6, 1:6))
hist(x, breaks = seq(min(x) - 0.5, max(x) + 0.5, 1), main = "")
两次掷骰子得到的结果之和的分布
两次掷骰子得到的结果之和的分布

样本量为2的时候真的得到正态分布了么?CLT确实有它的神奇之处,但还没那么神奇。以上结果如我前文所说,仅仅是由于样本空间中的元素的对称性,所以得到了一副对称的图形,看起来像是正态分布,如果我们再用另外6个数字试验一下,马上就看穿这种迷惑了,我把最后一个数字换成15再看直方图:

d = c(1, 2, 3, 4, 5, 15)
x = rowSums(expand.grid(d, d))
hist(x, breaks = seq(min(x) - 0.5, max(x) + 0.5, 1), main = "")
6个数字中2个样本之和的分布
6个数字中2个样本之和的分布

如何?“正态分布”哪里去了?客官是否能看清骰子中的中心极限定理谎言了?

2、直方图上添加正态密度曲线有误导性

在我还会用SPSS的年代(大约是SPSS 15.0),我就发现SPSS的一个荒唐之处,它的直方图有个选项,可以控制是否添加正态密度曲线。我们被“正态”毒害得多深,从这些软件设置就可以看出来。为什么只能加正态密度曲线,而不能加数据自身反映出来的核密度估计曲线?换句话说,数据的分布一定要是正态么?

我们看CLT的模拟,很多人也喜欢在直方图上加上正态分布密度曲线,这有一定道理:可以看直方图跟正态密度的差异有多大。然而,我们却很少见直方图上的核密度曲线(注3)。既然是要作比较,就拿最可比的去比,比如曲线对曲线,直方图对直方图,人眼本来就是不精确的测量工具,那么制图者就应该提供尽量准确而方便的参照系。

3、给出拟合好坏的指标

综上,我在两年多以前已经把这些CLT模拟的想法写成函数收录在R包animation中,参见函数clt.ani()。代码及输出如下:

if (!require(animation)) install.packages("animation")
library(animation)
ani.options(interval = 0.1, nmax = 100)
par(mar = c(4, 4, 1, 0.5))
clt.ani()
中心极限定理模拟:从指数分布到正态分布
中心极限定理模拟:从指数分布到正态分布

这个CLT模拟的过程很简单:给一个总体分布(默认为右偏的指数分布),在给定样本量n时不断重复抽样分别计算样本均值,一直这样计算obs个均值,并画出它们的直方图和相应的核密度估计曲线;然后随着n增大,看相应的样本均值分布如何。此外,我使用了Shapiro正态性检验来检验这些均值的正态性,并把P值取出来画在下半幅图中。这样我们就很清楚地知道,对于每一种样本量(n = 1, …, 100),我们的样本均值究竟离正态分布多远。此处P值就充当了一个拟合好坏的指标。可以看出,上面的动画中,当样本量n超过20之后,P值会普遍偏大,也就是样本均值的分布和正态分布比较接近(严格来说,是“不能拒绝正态分布”),但也不能保证样本量大就一定意味着正态分布,譬如上图中n=40的时候P值就很小。

三、小结

正态分布在统计中的地位如此之重要,以至于人们几乎认为正态是一种自然而然的分布,本文想说明的是,正态有它的自然性(参见注1),但我们不能逮着两边低中间高的东西都叫正态(参见小学课文《小蝌蚪找妈妈》);做模型或分析数据之前,先清空脑子里的这种“正态教义”,用事实说话。

另外,前面胡江堂挖了个大坑,大家跳进去争论了半天SAS和R的问题,我在这里也挖个小坑:就模拟而言,如果有SAS用户愿意研究一下,我想知道朱老师的SAS模拟代码是否有改进的余地。SAS自70年代创立,过了二十多年才引进“函数”这种杀人越货编程必备之工具,一直以来都是“宏”的天下,说它“恐龙”应该也不冤枉吧。朱老师幻灯片中的长篇SAS代码,如果用R改写,应该都不会超过三五行甚至一行(注4),而反过来可能就麻烦了,比如上面的动画(用SAS怎么写?),它用R之所以方便,是因为R大量的统计相关的函数可任意调用,而图形也很灵活,可以把P值动态写在图例的位置,也可以愿意把样本量n作为下标动态写在$\bar{X}$旁边。

脚注

注1:以我的浅见,统计学中绝大多数极限正态分布的来源都是泰勒展开,主要是因为泰勒展开中有一个平方项,这个东西和正态分布的密度函数(或特征函数)看起来形式相似,再加上一些对高次项的假设(使它们在极限中消去),正态分布也就来了。另一个典型例子是极大似然估计的渐近正态性:泰勒展开中一阶导数为零,剩下二阶像正态。此外,正态分布有个独特的特点,就是它的密度函数和特征函数长得很像,冥冥之中也主导了很多统计量的分布(这一点我还没太想清楚)。

注2:我很少见统计学家强调这种想法的重要性,仿佛大家都埋在公式堆里都推导得倍儿高兴。有两个例外:一个是Brian Ripley教授,我在看他的一份幻灯片Selecting amongst large classes of models的时候发现他竟然回顾了AIC的假设条件——这是我见过的唯一一个讲模型选择时会回到数学根源的人;另一个是Ripley的老搭档,Bill Venables,他在他那著名的手稿Exegeses on Linear Models中竟然从泰勒展开来说线性模型的来历,这也是我在所有我读过的线性模型相关的文章和书中看到的唯一一份从泰勒展开角度谈线性模型的材料。我个人并不喜欢数学,但我喜欢看思路。

注3:我的观点是,画直方图尽量加上核密度估计曲线,因为它能刻画数据在任意位置上的密度大小,而直方图的形状则完全受制于分组边界的位置。当然,核密度估计曲线也取决于核函数,但大多数情况下,核函数不会对曲线的形状有太大的影响。不过邱怡轩在前面“有边界区间上的核密度估计”一文中提到的问题很值得注意。

注4:伯努利随机变量之和的分布。SAS代码:

symbol;
goptions ftext= ctext= htext=;

%macro a(n);
data rv;
retain _seed_ 0;
do _i_ = 1 to &n;
binom1 = ranbin(_seed_,&n,0.8); /*ranbin函数返回以seed、n、0.8为参数的二项分布随机变量*/
output;
end;
drop _seed_ _i_;
run;
%mend;

%a(100);

proc univariate data=rv noprint;
var binom1;
histogram/normal( mu=est sigma=est);
inset normal ;
run;

R代码:

a = function(n) hist(rbinom(n, n, 0.8))

什么叫向量化编程?这就是R它爹S语言诞生的原因之一:统计学的编程不应该涉及到那么多繁琐的循环。

关于谢益辉

RStudio码了个工,Iowa State University统计系博了个士。统计之都网站创办者;研究兴趣为统计图形及数据可视化,对统计模型方法的发展感兴趣但不喜欢纯粹抽象的数学理论,以直观、实用为学习标准;偏好以R语言为工具;Email:xie@yihui.name;个人主页:http://yihui.name

从中心极限定理的模拟到正态分布》有50个想法

  1. 终于看到有业内人士进行一些正本清源的工作了。统计中揣着糊涂当明白的人太多了,当然也有很多真明白的人, 但我有个疑问,为什么这些明白人不去指出一些基本概念中的误导呢?比如说方差的含义和表示方差的公式,两者之间好像是有点小问题。

  2. 极限理论有两个关键的定理:central limit theorem 和 law of rare events。
    前一个渐近分布是正态,后一个渐近分布是poisson。
    并不是所有的“和”都能渐近正态。

    统计里面大部分的渐近正态性由taylor展开得到,但是小样本时候也有所谓的Bartlett Correction,也就是taylor展开到稍微高阶的项。
    投色子,两项的和怎么可能渐近正态呢?看画出的直方图尾巴那么厚,目测都知道不是正态,即使他是对称的。《金融计算与建模》的作者估计还未打通任督二脉。

    1. 中心极限定理没有错,但是样本量为2的时候话中心极限定理,这是没有真正理解“极限”的含义。

  3. 谢谢益辉的工作,记得第一次R语言会议席间,还给我更正过这方面的错误呢。今天读了这篇文章,加深了理解。

  4. 谢先生在这里还是以代码的长度论输赢, 恐怕还是书生气了一点。 做industry的东东,关键是能不能get something done with reasonable resources and cost. 我的R 和SAS 都不怎么样,都用最平常的写法, 现在分析一个50G的数据, 用SAS 在台式单机(2G内存, 200G硬盘)15小时搞定。到LINUX SERVER (8 CORES, 32G内存)跑R, 怎么都是死翘翘。 不要高诉我在R里用call C PERL什么的也一样, C和PERL能做到的不等于R也能做到, 也不是每个分析师都有时间学C

    另外关于sas作图,可以看下面的连接里的例子。
    http://www.robslink.com/SAS/Home.htm
    当然这些图肯定还是比不了R的功能, 我只是说SAS如果学好了, 一样可以作出漂亮的图, 与软件无关, 与人有关。

    我其实对R 和SAS都感兴趣,他们各有长处。 我所反对的是根据一些不全面的信息而作出误导性的结论。

    1. 如果我真的以代码长度论“输赢”的话,那我确实太书生气了——我没有推广到输赢的意思。不过如果有简单的办法能实现目的,为什么要选择臃肿的办法呢?如果代码长短没关系,那何不统一退回到汇编语言呢?用你自己的评判标准,在本文的模拟中,R的“resource and cost”都免费,而且实现起来更方便,要啥有啥,SAS写起来像恐龙语言,而且输出往往是有用没用的一大堆,你用SAS也知道,给3个数字都能输出好几页报表。

      你若是统计学家,论证也应该有点统计学的风格:一个学统计的人有多大概率会去分析50G的数据呢?举例子是个很好的论证方式,但抽样要有代表性。另外,你50G数据硬要装到32G内存中,这不是故意调戏你的服务器么……这问题有两种可能的解决方式:一、拿昂贵的SAS购买费用去多买一些便宜得像白菜的内存;二、看看R处理大数据的包,比如bigmemory包。无论如何,不能用“不知道”代替“不可能”。

      我也不想总在这个问题上争来争去,有人愿意花钱愿意多写代码,那也没办法。我这篇文章写出来就是避免空泛的争论,大家如果有争议,可以针对具体的工作来讨论。我知道我的文字看起来很偏激,这非我本意,我只能说就我的工作而言,R可以包打天下,我没有任何用SAS的必要,我一再强调,读者要切记,这是我的个人观点,至于我是不是个有代表性的统计工作者,只能留给旁观者判断了。

      从始至终,我也没说“SAS不能做什么”,只是在说“R能做什么”,这里面如果有什么误导的话,我愿意洗耳恭听。

      1. 其实我真的不愿发这些帖子, 因为我其实对R非常感兴趣, 也在积极努力的提高。因为R虽然有些缺点(依赖内存,package良莠不齐,文档差), 但是又有那个软件没有缺点呢? 也是R是缺点比较少的语言了把。

        但是我现在还很少看到一个SAS programmer去攻击R(也许SAS PROGRAMMER 大多是老帮菜, 很少看BLOG, 哈哈), 相反R论坛上充斥了对SAS的以偏概全和道听途说的攻击, 已经产生了很多误导效应,所以不得不厚着老脸出来提供一点 不同之声。

        谢先生说50G的数据对统计师来说是OUTLIER, 我不得不说谢先生对当前DATA MINING界发展的把握是略欠火候了。 我证券交易所的朋友随便一个dataset轻轻松松就20G以上, 如果是double click 的online 数据, 平均体积是35-60 GB, 我是搞生物的, 下一代基因测序技术的出现使每个病人产生的数据量达到TB。 这些都是要靠统计师在挥汗如雨的操作。 你在感叹内存的廉价, 殊不知真实世界的数据增长大大超过了内存的增长。

        总之R快速增长已经使其成为现代统计师不可不学的工具, R的重要性和优点根本不需要通过攻击SAS来得以彰显, 因为象牙塔以外的世界是如此的复杂和挑战, 我们太需要各种不同的工具并存去发挥他们的长处。

      2. 多谢指点!我知识范围有限那简直是一定的(a.s.),不过呢,仍然是我的概率观点:搞统计不等于做Data Mining,换句话说,DM的人占统计工作者多大比例呢?我知道相对来说,做DM的人手上的数据都会比较大,我也并非没有耳闻,前年我去新浪溜达的时候了解过一些情况(比如你说的TB单位的数据),但并没有参与他们的工作,我自己仅仅在某政府部门用R亲手做过5G左右的数据(在普通的PC机上,而且是Windows XP系统,大数据并不代表所有数据都有用或一定要整体读入)。

        要说误导,我真不知道“误”在何处,有谁说过“SAS不能做什么”没有?本文中我都给了SAS代码,说它可以做,但问题是孰繁孰易。抱歉,我似乎有点不见棺材不落泪……

        看!又跑题了……本文的例子中,哪里有50G的数据啊,两个骰子36条数据,手指头都能数过来,若要谈大数据的处理问题,各位作者们谁有兴趣可以写下一篇文章,再专门讨论。

    2. R语言的精炼你不要不承认,你要是非抬杠,那岂不是大家都去学机器语言了,还要现在的高级语言干什么?

      “与软件无关, 与人有关。 ”我喜欢你这句话,俗话说没有金刚钻,不揽瓷器活。人笨,不要嫌刀钝。你那50G的例子是你自己搞不定,不代表别人也搞不定。个案,不具有普遍性。要是SAS的分析师都成了门外汉,我看SAS也该近黄昏了。

    1. K-L信息量这个东西确实是衡量两个分布差异的工具,但我一直觉得它有个小问题,就是它没有上界,就如同方差衡量变异程度一样,没有界限的东西就很难知道到底差多大算大、多小算小(当然0意味着两个分布a.s.一样)。

      1. 这种准则性的东西还是两两比较有优势,所谓不怕不识货,就怕货比货。

  5. 我来打打擦边球(这好像是我的一贯习惯)

    — 我顶评注1. 这一点是做asymptotical analysis的必备道具,基本上我见过的证明一个估计量是渐进正太的论文,都是按照这个思路来的:1.泰勒展开,2.证明二次项主导其他所有的项。。。这也是我厌倦理论分析的原因。

    但是还有另外一种分析,走得完全是另外一条路子,那就是不仅假设n->infinity,同时,模型的变量也->无穷,这个就很难证了。而且也更好玩。

  6. 谢谢益辉,貌似明白了一些我以前的有疑惑的地方,但还是有些不太明白的地方。
    我再仔细看看,喜欢你这种 正本清源 的精神 #_#

  7. 鄙人觉得R和SAS的起因本身是不同的. SAS的设计初衷是极度应用性的软件包,为了用户能写几行代码就可迅速得到答案, 而不是为了开发. R的初衷是做开发环境. 你可以自己开发像SAS里的那种几行代码出一大摞结果的包, 也可以当它是一般编程语言来使用做任何想做的事情. 我觉得可以把SAS比喻成Financial Calculator, 把R比喻成Scientific Calculator. 比如你在算金融中的特定问题时, Financial Calculator预设的功能可以让你几个按键便出结果. 但一般科学函数计算器则要从first principle开始一步一步地算, 有些高级一些的函数计算器带macro储存, 那也可以把这一步一步的案件存成macro,以后直接调用.
    总的来说我觉得”分析员没时间学C语言”这种说法算一种不太好的借口. C语言影响了几乎全部当今的主流编程语言, 这是任何程序员的基本功中的基本功, 如果连C语言都不愿意学的话, 那实在是不能算具有专业水准的.

  8. 我有一点还在惶惑,中心极限定律和大数定律一直都是放在一起讲的,那么中心极限定律和大数定律又是什么关系呢,期待有人讲一讲。

    1. 大数定律只是一个点的性质,而中心极限定理是分布的性质,后者比前者更深刻。

      1. 我觉得粗糙地理解的话可以认为CLT蕴含LLN. CLT给的是sample mean的极限分布, 并且方差趋向于0. LLN给的是sample mean的收敛极限. 那既然方差收敛于0了, 那粗略的看应该就是随机变量按概率收敛到一个点了.

    2. 中心极限定理是 大数定理往后的延伸,先有大数定理,才能有中心极限,不然中心极限定理的最核心的就是样本量的扩大化后,呈正态分布 ,这个就不能成立了……

  9. 益辉,那段模拟二项分布的SAS代码,最后缺了结束宏的%mend;。说代码的简洁,其实,这段宏编译之后,%a(n)就可以当作函数用,n是唯一的参数,这就可以相当于附上的R代码a = function(n) hist(rbinom(n, n, 0.8))——这段R代码里面的几个函数,也是背后编译好的。正文的那段比较,可能有失公平。

    SAS的data steps脱胎于以前的大机语言PL/I,有明显的过程式编程风格(C、非OOP的C++等等),而R是向量式的编程语言,个人认为,对传统的程序员来讲,SAS会更亲切,R则是让统计学家感到舒服的语言。

    1. 江堂兄指出了要害:
      SAS中的宏=R中已有的函数
      其实从简洁性的角度来说是差不多的。

      大家喜欢SAS就用SAS,喜欢R就用R.

      我个人的观点还是:找工作用SAS,搞研究用R

      1. 这个等号要换成约等号。如江堂所说,SAS的编程风格是过程式的,把一个工作的流程一起写完然后一起提交代码,顺序执行,这是宏(或批处理)的特征,表面看起来宏确实挺像函数,但有个本质的区别就是函数可以“批处理”一些事情,它还能返回值,而宏最多只能创造一些全局变量供后面的过程使用。

        在注4的SAS代码中,画直方图只能根据全局的数据rv来做,宏只能这样:

        %a(100);
        proc univariate data=rv noprint;

        而不能直接把数据传给需要它的地方,如

        proc univariate data=%a(100) noprint;

        用宏写一段小程序没问题,如果是大型开发,而程序语言不支持函数,程序员肯定会疯掉。当然,SAS在问世三十多年后,现在终于可以支持函数了。

      2. 益辉说SAS”不能”把“把数据传给需要它的地方”, 看来你对这段宏的理解还欠火候, 让小的斗胆来评论一下:
        proc univariate data=rv noprint; 里面的rv, 表面看来是一个固定的数据, 但是他的内容是由%a(n)里的n (也就是这个SAS 函数/宏的唯一参数)来决定的, n=100, rv 就是一个长度为100 的vector, n=1000, rv就是一个长度为1000的vector, 这个rv就是%a(n)的唯一输出变量, 也就是益辉所说的返还值。

        那么宏和函数的差别到底在那里呢? 宏的本质是文本替代, 换句话说, 如果用一个%do 1=1 %to 100的 循环来运行这个%a(n), cpu在编译的时候其实是实际上产生了100X的代码, 然后执行。 但是函数是抽象的, 编译好了之后也只有一段code, 但是CPU要把这一段code动态运行100次。 所以对CPU来说宏和函数的差别是本质的。 但是对于终用户来说, 都有输入参数和返回值, 宏和函数的功能是很相似的。
        由于宏的历史比较悠久了, 现在嘲笑宏好像也是很时髦的事情。但是应该承认宏在现代计算语言中还是有一席之地的。 C就是支持宏的, JAVA不支持宏也时常为人诟病。 有兴趣的还可以看一看ROSS IHAKA的一个讲座:
        http://www.stat.auckland.ac.nz/~ihaka/downloads/Singapore-2008.pdf
        IHAKA最近在研究超越R的下一代统计平台, 在里面宏的重要性几乎是和函数并立的。

    2. %a(n)确实可以当做rbinom(n,n,0.8)来用,但生成n个随机数是统计模拟中最常见的功能,如果在庞大的SAS中还需要一个个生成,那未免也太悲剧了。这个比较并非针对功能问题——长短代码都可以实现类似的功能,而是(统计)程序员的工作量。

      %mend;已经加上,谢谢。

      1. 我的R只是”hello world”的水平,先说说看。文中SAS的随机数生成函数ranbin(),与R函数rbinom()类似——SAS也不需要“一个个生成”随机数,而那段宏大致与hist(rbinom())相当。

        区别在于,R用hist()生成直方图,而SAS需要用一个proc univariate来实现。

      2. 其实SAS “一个一个”生成随机数, 这个说法是没错的。
        binom1 = ranbin(_seed_,&n,0.8); CPU产生随机数, 下面的OUTPUT是把该记录写入硬盘。sas的工作模式是从硬盘读一条记录, 在CPU/RAM操作之, 然后将其写回硬盘, R则是讲整个object读入内存, SAS这个特点是他在操作大数据方面很有优点, 但是代价是整个过程是 data centric, 逻辑比较不直观。 这里又不能不提到 Ross Ihaka的另一片雄文了。
        http://www.stat.auckland.ac.nz/~ihaka/downloads/Compstat-2008.pdf
        在里面他提到SAS这种“one record at a time” 也要被整合到他的所谓下一代计算语言中了…….

      3. to xin,

        你是对的,ranbin()产生一个随机数,n个循环则产生n个随机数。我提到SAS不需要“一个个生成”随机数,是认为“一个个”这个字眼显得太手工,循环很方便完成这个任务,而循环如果不是明显写出来的,就是暗含在里面的。

      4. 所以说向量化编程可以减轻统计工作者的编程任务啊,R把这些常见的功能都封装到底层语言中,代码写起来既轻巧又不失效率

    1. 谢谢,这个建议很好,固定坐标轴范围之后可以明显观察到:样本均值的分布越来越集中于一个点,即:x轴越来越窄,y轴越来越高。这是自然而然的结果,因为最终样本均值会以概率收敛到一个点(总体均值)。

      这里面有个问题就是,y轴的范围要画完图之后才能确定,事前没法确定,虽然也不是太麻烦(用hist()事先计算一下就可以),但我不太想增加代码了……

  10. 你好,想要做CLT的一个仿真gif图,通过google找到这来了。刚刚开始学2天R,问的问题可能比较菜,见谅!
    clt.ani里面似乎只能做独立同分布随机变量的演示,如果随机变量本身属于不同的分布,或者分布类型相同,但方差等不同该怎么做?还有是不是只能做一些常用分布的演示,如果我已知了分布的表达式想要做演示该怎么做?

    我再把我的问题说具体些:有一系列相位随机的正弦运动,每个运动的幅值当作一个随机变量,想演示这样的n个随机变量的和的分布随n增加的演变过程

    1. 你说得对,clt.ani()只能做独立同分布的演示。如果你已知了分布的表达式,那么下一步需要知道的是如何从那个分布中随机抽样,就像rexp()从指数分布中抽样、runif()从均匀分布中抽样一样。

      我不太了解你说的正弦运动,你说的随机变量是振幅?(我不知道幅值是什么)随机变量的分布是什么?

      1. 问问题果然不那么容易,还是没问清

        那么如何从一个已知表达的分布中抽样呢?我搜索不到相关的资料,对rexp()这类函数的介绍只有如何用它,而没有函数原理的说明。

        y =A*sin(t) t是一个均匀分布的随机变量,A是系数,我要演示的就是N个y的和的分布情况。运到的问题是:

        1.这样表达式所表示的y取值的分布不是常用的分布,不能直接用rexp() runif()这样的函数进行抽样
        2.每个分布的类型相同,但是方差是不同的

      2. 我不知道均匀分布的随机变量在取了sin()之后是什么分布。这里如果A是常数,那么CLT的演示是没有问题的,如:

        library(animation)
        A = 1
        ani.options(interval = 0.2, nmax = 100)
        par(mar = c(4, 4, 1, 0.5), 
            mgp = c(1.5, 0.5, 0))
        clt.ani(FUN = function(n) A*sin(runif(n)), 
                mean = NULL)
      3. 感谢你耐心解答我的问题^_^ 现在大致明白了,这里有几个参数的意义你看我理解的对不对:
        1. nmax表示动画的帧数,同时也决定了样本量的数目,因为动画是从样本量为1,逐渐增加到nmax为止;
        2. obs(默认值300),是样本均值的数目,或者说是对每个样本抽样了obs次,决定绘制直方图时的精细程度(但是试了下图片上没有变化)
        3. 自定义函数中的n呢?和obs是什么关系?和nmax有什么关系?感觉好像n就是去了obs的值。如果我吧常熟A写成用n索引的一个数组,是不是能实现对不同分布的演示,即clt.ani(FUN = function(n) A[n]*sin(runif(n)), mean = NULL),试了下好像不行,估计这个n的意义理解的不对,或者clt.ani只能做同分布的演示

      4. 1、完全正确。
        2、完全正确。300已经够大了,所以你要是再增加,效果应该也差不多,但是要是减少的话,应该就能看出来直方图没那么精细了。
        3、这个动画的过程是这样:对于特定样本量n,重复随机生成obs批样本,每一批样本计算一个均值,这样就出来了obs个均值,也就可以画一个直方图了。如此,逐渐增加n,也就出来一幅幅直方图。n便是样本量。如果你写A[n],那么意思是从向量A中取出第n个元素(只是一个数字),也就是随着样本量增加,sin前面的系数在向量A中向后移动。

  11. SAS的macro和R的函数在是不是文本替代的问题上其实区别不大。R的函数里面其实很大程度上也是经过parse的所谓“文本”。我觉得更大的区别在于variable scoping。SAS的macro variables有一定的规则来决定是%local的还是%global的,但是它的macro variable本质上不是data set,而它的data set是没有这种默认的scoping——如不特殊指定,一切都在work库下面。所以,写大程序的时候,数据名冲突很容易发生。R的函数则不用去分什么是数据,什么是变量,一切皆为object;除非特殊指定,几乎所有scoping都是局部的,因此,用R写函数要比用SAS安全许多。至于更细节的(enclosing) environment与evaluation frame的问题对大多数使用R的人来说,都不必了解很深就可以写出非常reliable的代码。Being functional is highly desirable to most statisticians.

  12. 你不应该死抠字眼,通常我们所说的正态分布随机数或者说我们需要的正态分布随机数都不是完全正态分布的。
    真正的正态分布的随机数只会被用在理论研究和一些消耗随机数很大的项目中,大部分情况下,随机模拟需要的正态分布随机数不能是完全正态分布的,因为其中有些数据在实际生活中不会出现(实际上很少统计结果真的是正态分布的),这些数据必须被筛选掉(比如在±2σ之外的数),使用中央极限定理是一种优化,他能很好的控制随机数的范围和分布,一来产生随机数所需的计算量显著降低了,而来产生的随机数符合实际需要以至于根本不需要筛选,我们甚至可以省略变换的需要。

  13. 为什么有时候用hist()画出来的密度直方图的纵坐标刻度范围不是0到1???

      1. 画密度的时候,要保证的是直方图的矩形面积加起来为1,这并不代表单个矩形的高度要小于1。如果数据的范围很小,那么相应的密度值也就要大一些。例如:

        hist(runif(1000,0,.1), freq=FALSE)

        换句话说,密度不是概率,没有[0, 1]的要求。

  14. 话说就没一个POSITION要求用到R和SAS的么?

    截然不同的2条道路啊

邱怡轩进行回复 取消回复

邮箱地址不会被公开。 必填项已用*标注