前一段时间，John Baez 在自己的主页上更新了一篇文章名为 The Beauty of Roots，这篇文章之后在“科学松鼠会”上被 转载。上面提到了曼德布洛特集（Mandelbrot set），根据其发明者法国数学家 Benoît Mandelbrot 而命名。

曼德布洛特集是一种分形，从一般分形性质来说：

客观自然界中许多事物，具有自相似的“层次”结构，在理想情况下，甚至具有无穷层次。适当的放大或缩小几何尺寸，整个结构并不改变。不少复杂的物理现象，背后就是反映着这类层次结构的分形几何学。

常见的曼德布洛特集是这个样子（分辨率原因，部分细节显示不够）：

假如我们把这个集合的下半部分（最下边的小块）分割出来，就是这个样子（8倍放大）：

由于分辨率的提高，所以显示了第一幅图中并没有显示的细节。

继续放大，上图的左上部分的那个小枝（6倍放大）：

再把上图最靠近左边的那个小枝——放大（50/3倍放大）：

继续放大最左边的小枝，似乎在末端又出现了一个类似的小枝（5倍放大）：

如果继续放大下去可能还是这个样子 ：）

注释：

1. 最后一张图相比第一张图来说相当于局部放大了 4000 倍。
2. 高质量的矢量绘图数据量比较大，R 处理起来有些问题，只好使用局部放大的方式。
3. 更多的使用其他软件绘制图形可以见：http://commons.wikimedia.org/wiki/Mandelbrot_set

最后广告一枚：

COS 统计图形和可视化版 ——用 R 玩分形。

还记得第一次看到水立方时的惊讶么？

是什么这么吸引我们？是有如天空般的颜色？还是那气泡似的形状？

从水立方的外墙结构上看，不但外观美观，而且十分紧凑。水立方外墙为什么会有这样的性质，是因为它上应用了一项最优化的技术，即Voronoi 原理。

Voronoi 图也常常被称为 Dirichlet 格局（Dirichlet tessellation）。通俗讲，其原理是一项从点到面的技术。它的每个多边形只有一个”生成点”，而这个多边形上的每个点到”生成点”的距离总是比到其他”生成点”的距离要小（是不是想到了 K-means 算法？）。

在建筑设计上，有设计人员争论这类方法定义为“参数化设计”。认为这种方法不能同传统的、依靠灵感的设计方式相比，因为这种方法高度依赖计算机，只需要简单改变若干参数就能得到设计方案。但这个论断，恰恰忽略了“参数化设计”背后的数学意义。

既然 Voronoi 是一种最优化的算法，那么除在建筑学上给我们带来的美轮美奂的视觉效果外，它在空间统计上同样也有应用。

下面，我根据各个省会城市（包括香港、澳门）的地理位置，利用 Voronoi 原理，计算每个省最佳控制范围（使用红色的线条标记）：

虽然理论值（最优）和现实值（行政区划、地理）总有差距，但是，比较一下会发现一些值得探讨的现象：

• 内蒙古应该好好的规划一下，从东边到西边实在太远了，把西边的划给宁夏可能好点；东边划给北京、东三省；
• 河北北部，不论是属于北京还是天津都会好些，记得我小的时候，宁可去北京、天津，也不乐意去遥远的省会–石家庄；
• 青海应该把甘肃的北部包括进去；
• 上海、香港、澳门有一部分管辖区也也不错么：）

整体上看，大部分省的行政区划还是符合 Voronoi 原理。也就是说，单纯从空间距离的角度来看，我国的行政区划还是比较不错的。，

很多时候，在社会调研中会出现一些小数（或百分数），而这些数字背后隐藏的信息也常常被统计人关注。比如 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

COS 上有人问过如何求1～100的素数。虽说这个问题没准就是计算机系大一新生的一道作业题，但对我这个几乎没有任何 C 编程经验的人来说，似乎还是有些挑战。花了几分钟整理了一下思路，既然素数的定义是只能被1和自身整除，那么： 1、如果第 n 个数，不能它前面的所有的数字（不包括1）整除，那么即为定义。但需要遍历所有数字，效率肯定不好。 2、如果 n 不能被 n 前面的所有素数整除的话，那么 n 应该是下一个素数（后来知道这个是算术基本定理）。 根据第二点思路，写出 R 代码：

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。 马上 cloud_wei 给出了另外的实现方式：glm包的 isprime 函数（这个包似乎没有 Windows 版本）； 接着 jo3vul31l3 在回帖中给了最优的解法，即埃拉托斯特尼筛法：

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)
}

在这以前我一直以为 素数 越到后面会越”稀疏”，但事实是这样（1：10000区间的素数）：

第一幅图的红色线是顺手做的一个线性回归拟合线；10000以前的素数几乎是平均的出现在各个区间（breaks = 100）。

Hierachial clustering 这类的方法最大的问题就在于对计算机内存的占用，当然对计算量也是一个严重的考验。n×m  的矩阵要进行 n×n 次计算，n^2 这个函数估计大家应该记得很清楚，汗~~~~

Partitioning clustering 好处就是将计算量由x^2变为 x ，即线性函数。而其典型代表就是kmeans：

图一是kmeans 算法中，增加变量（由1到100）的运算时间，图二是增加单维变量长度（x×1到x×100）。如果增加变量情形下，估计HC应该和PC 一致，但如果同时增加case的话，HC必然崩溃……

九连环这个玩具，最早在高中的时候摸过，记得解它的时候因为很多重复步骤，所以恨不得要把它拆掉，还好最后解开了，不然我又得加条”亵渎古人智力”的罪名，呵呵。闲话少扯，下面是 n 连环实际步骤数的求法：

gg <- function(n){
a <- numeric(n)
a[1] <- 1
a[2] <- 1
for(i in 3:n){
a[i] <- a[i-1] + 2*a[i-2] + 1
}
a}

比如九连环的实际步骤：

gg(9)

[1]   1   1   4   7  16  31  64 127 256