Beta

很多时候，在社会调研中会出现一些小数（或百分数），而这些数字背后隐藏的信息也常常被统计人关注。比如 COS 主站上的这篇文章--《从调查报告中的比例数字说统计人如何甄别统计假象》，yihui 生动地为我们展示了一种考量问题的思路。

正如文章中所说的，如果我们对数字足够敏感的话，很容易判断出 0.6667 的分数是 2/3 ，0.625 的分数是 5/8，0.14286 的分数是 1/7。但我们的经验毕竟有限，不可能穷尽所有的数字，通过一个算法来确定分数是十分有必要的。

法里序列（farey sequence）也是考虑这类问题的一个角度。如果给定法里序列的 n 足够大，那么我们必定能够将逼近出一个和小数相等的分数F_i[j]。

法里序列 F_i （i=1 到 n）：

F₁ = {⁰⁄₁, ¹⁄₁}

F₂ = {⁰⁄₁, ¹⁄₂, ¹⁄₁}

F₃ = {⁰⁄₁, ¹⁄₃, ¹⁄₂, ²⁄₃, ¹⁄₁}

F₄ = {⁰⁄₁, ¹⁄₄, ¹⁄₃, ¹⁄₂, ²⁄₃, ³⁄₄, ¹⁄₁}

F₅ = {⁰⁄₁, ¹⁄₅, ¹⁄₄, ¹⁄₃, ²⁄₅, ¹⁄₂, ³⁄₅, ²⁄₃, ³⁄₄, ⁴⁄₅, ¹⁄₁}

F₆ = {⁰⁄₁, ¹⁄₆, ¹⁄₅, ¹⁄₄, ¹⁄₃, ²⁄₅, ¹⁄₂, ³⁄₅, ²⁄₃, ³⁄₄, ⁴⁄₅, ⁵⁄₆, ¹⁄₁}

F₇ = {⁰⁄₁, ¹⁄₇, ¹⁄₆, ¹⁄₅, ¹⁄₄, ²⁄₇, ¹⁄₃, ²⁄₅, ³⁄₇, ¹⁄₂, ⁴⁄₇, ³⁄₅, ²⁄₃, ⁵⁄₇, ³⁄₄, ⁴⁄₅, ⁵⁄₆, ⁶⁄₇, ¹⁄₁}

F₈ = {⁰⁄₁, ¹⁄₈, ¹⁄₇, ¹⁄₆, ¹⁄₅, ¹⁄₄, ²⁄₇, ¹⁄₃, ³⁄₈, ²⁄₅, ³⁄₇, ¹⁄₂, ⁴⁄₇, ³⁄₅, ⁵⁄₈, ²⁄₃, ⁵⁄₇, ³⁄₄, ⁴⁄₅, ⁵⁄₆, ⁶⁄₇, ⁷⁄₈, ¹⁄₁}

但这个过程会比较麻烦，F₁₀₀₀ 已经达到300927 个数字。幸好 R 中的 MASS 包提供了 fractions 函数。这个函数使用有理近似的方式，将小数转化为分数（连分数）形式。比如《从调查报告中的比例数字说统计人如何甄别统计假象》中提到的 29.1667% 这个数值：

> fractions(0.291667)
[1] 7/24

不过，既然是近似算法，这个函数对小数的精确度要求还是蛮高的，而且最好不要用无理数去逗人家。

> fractions(pi)
[1] 4272943/1360120

SPSS 在首页显著位置公布 An important message for our customers and partners，同 IBM 共同宣布 SPSS 被收购的 definitive agreement。这 SPSS 改名还没几天，又有了这么大的动作，BI 界不太平啊！

IBM 的 Press 里有段话很有意思：As companies attempt to control costs and use resources more wisely, IDC estimates that the worldwide market for business analytics software will swell to $25 billion this year, growing 4% over 2008.(1)

深圳福彩3000万巨奖诈骗案发生后，好多朋友开始质疑中国福利彩票的公正性，并向我询问福利彩票是不是真的如传闻那样--福利彩票自产自销？。由于身在中福彩的原因，很多支持福彩公正的内情不太合适发布在网上，但--数据一汇总，上帝也会发笑。这篇博文从一等奖中奖概率分布的角度给各位童鞋聊聊福彩的公正性。

现在每一期的福彩双色球销售额大概是2亿左右（已持续很长时间），每注双色球为2元，就是说每期双色球的销售大概会有1亿人次参与。

在随机选择，且每次抽取都是等概率的的假定条件下，理论的重复彩票注数的分布（0注至41注）如下：

一等奖是完全随机出现的，那么在假设条件下，一等奖同时出现五注的概率最高，六注、四注其次，再次为七注、三注，类推……当然理论上，没有中奖（0注）和中12注以上的概率相比其他情形低了很多。

那么我们再看一看福利彩票双色球每期中一等奖（2008年1月1日至2009年7月23日）的实际分布情况：

有童鞋看出端倪了：实际数据的分布同理论上的分布是不一致的！理论上的一等奖出现概率最大在出现五注的位置，而实际上一等奖出现的最大概率出现在了两注的位置。

为什么会出现这样的情况？主要是因为我们最开始的假设是有问题。

双色球每期销售会有1亿人次的彩民参与？不可能！双色球的覆盖度没有那么大。

一些彩民为了提高中奖概率（或者说迷信一些选号方法），会采用"复式"、"胆拖"、"倍投"等方式投注，当然大部分彩民还是会老老实实的买一注。综合考虑到这些因素以后，凭经验估计样本量应该为现在的一半左右，即5000万。这样看来较为合理的一等奖中奖概率理论上分布为：

此时理论分布同实际分布已经非常相似。

实际双色球一等奖分布的右边尾巴上恰恰显示了"复式"、"胆拖"、"倍投"的投注效果。

再插一句：

深圳福彩3000万大奖诈骗犯身份曝光一文中提到：

警方调查发现，程某先是编写了一个可以自动运行的木马软件，然后利用与福彩中心合作的机会，进入福彩中心机房，植入自动运行的木马程序。一旦摇奖结果出来，这个程序会自动将程某所购买的彩票修改成一等奖的号码。

这里可以推测程某天真地以为满足兑大奖的条件为：

• 数据库里的数据正确；
• 实体彩票存在。

恩，没有问题！但，这两个条件可是通过很多很多很多的手段来监管的。

以前用过几个国内的数据分析（挖掘）软件，每次溜达到商家主页，基本介绍都是诸如：

1. 某某公司依托某某大学……
2. 某某指定专用数据分析软件

这时候我就会想，这些商家怎么不仔细瞅瞅 R 官网的 Members&Donor，看看那几十所 Supporting Institutions ，或者20位 R 的核心团队成员的资料。不过回头一想，这些数据分析软件商家不这么说的话，又能说点啥。

R News 正式更名为 R Journal 后的第一篇文章是 Facets of R，作者 John M. Chambers 可不是一般人物，虽然在 R 的官网上只能看到他是 R 核心团队之一，但仔细查一查就会发现--他可是 R 的老祖宗。

John M. Chambers 自 1966 年供职于贝尔实验室，在 1981-83 年间 head of the Advanced Software Department , and from 1983-88, head of the Data Analysis and Statistics Research Department，J.M. Chambers 在1998 年因为 S 语言（R 语言的前身）获得 ACM（Association for Computing Machinery） 系统奖。这个奖项每年评选一次，一般由 IBM 提供10,000$ 的奖金（wiki 上给的似乎有错误）。

这个奖项有着广泛而深远的影响，下面列举了一些我们熟悉的项目：

• 1983 Unix
• 1986 TeX
• 1989 PostScript
• 1991 TCP/IP
• 1995 World-Wid-Web
• 1997 Tcl/Tk
• 1998 S
• 1999 The Apache Group
• 2002 Java

基本都不用介绍吧，几乎每天我的工作、生活、学习都会涉及到这些软件系统。

最后引用 ACM Software System Award 对 John.M.Chambers 以及 S 语言的评价：

The ACM's citation notes that Dr. Chambers's work will forever alter the way people analyze, visualize, and manipulate data . . . S is an elegant, widely accepted, and enduring software system, with conceptual integrity, thanks to the insight, taste, and effort of John Chambers.

参考：

• John.M.Chambers 的wiki 词条及其链接；
• 谢益辉博士的 《R 语言的 历史-背景、发展历程及现状》

COS 上有人问过如何求1～100的素数。虽说这个问题没准就是计算机系大一新生的一道作业题，但对我这个几乎没有任何 C 编程经验的人来说，似乎还是有些挑战。花了几分钟整理了一下思路，既然素数的定义是只能被1和自身整除，那么：

1. 如果第 n 个数，不能它前面的所有的数字（不包括1）整除，那么即为定义。但需要遍历所有数字，效率肯定不好。
2. 如果 n 不能被 n 前面的所有素数整除的话，那么 n 应该是下一个素数（后来知道这个是算术基本定理）。

根据第二点思路，写出 R 代码：

1
2
3
4
5
6
7

prime2 <- function(m){
    x <- c(2,3,5,7,11)
    for(i in 13:m){
        if(!any(i%%x == 0)) x <- c(x,i)
    }
    return(x)
}

即给出前 5 个素数，而后寻找第 6 个素数；再根据 6 个素数找第 7 个；类推……；直至 n。

很快又有两种解法：

1. cloud_wei 的另外的实现方式：glm包的 isprime 函数（这个包似乎没有 Windows 版本）
2. jo3vul31l3 给出了最优的解法，即埃拉托斯特尼筛法：

1
2
3
4
5
6
7
8
9

prime3<-function(m){
    x<-c(2:m) ; y<-NULL
    repeat{
        z<-x[(x%%x[1])!=0] ; y<-c(y,x[1])
        if(x[1]>sqrt(m))break
        x<-z
    }
    c(y,z)
}