《R-概率统计与模拟》的系列文章已经写了四篇。其中有经过整理，内容相对全面，并且有内在联系的文章；也有文内几个小节相对独立，没有内在联系的文章。后一类文章正是笔者在学习过程中当时当地遇到了有意思的问题，即时地去解决并记录下来的结果。本文也仍然是这一类风格。学习是零散的过程，等到积累了一定的素材，才能进行整理归纳。
本文是《R-概率统计与模拟》系列文章的第五篇，包含了三个小节：
彩票至少包含一组连续号码的概率用归纳法解决概率问题的一个例题多个独立且符合同一个伯努利分布的变量的和服从二项分布买过超级大乐透的朋友应该熟悉，这种颇受欢迎的彩票分为两区，前区是从1到35这35个数字中无重复地选择5个数字，后区是从1到12这12个数字中无重复地选择2个数字。选定的前（后）区数字与开奖的前（后）区数字重合得越多，奖金越多。
很多彩票分析网站在统计历史开奖数字的各种信息，希望能从中找出些许的规律。这些信息中很多都是在统计“频率”，而这些“频率”其实都是有理论的概率值的。比如，“前区号码中有连号的概率是多少”？也就是：从1到35这35个数字中随机等概率且无重复地抽取5个数字，这5个数字中至少有两个数字是连号的概率是多少？
这个问题可以抽象为：从1到 n n n这 n n n个数字中随机等概率且无重复地抽取 k ( for n ≥ 2 k − 1 ) k \ (\text{for } n \ge 2k-1 ) k (for n≥2k−1)个数字，这 k k k个数字中至少有两个数字是连号的概率是多少？ （可以简单推论，当 n < 2 k − 1 n < 2k-1 n<2k−1时，抽取出的 k k k个数字中有连号数字的概率是1。）
这是一个蛮有意思的概率题，我们可以计算理论值，也可以用模拟的方法去计算其近似值。
其理论值是： 1 − ( n − k + 1 k ) ( n k ) 1 - \cfrac {\binom{n-k+1}{k}} {\binom{n}{k}} 1−(kn)(kn−k+1)
对于超级大乐透前区， n = 35 , k = 5 n=35, k=5 n=35,k=5，所以连号的概率是： 1 − ( 35 − 5 + 1 5 ) ( 35 5 ) ≈ 0.4766043 1 - \cfrac {\binom{35-5+1}{5}} {\binom{35}{5}} \approx 0.4766043 1−(535)(535−5+1)≈0.4766043
用代码模拟的近似值是0.4766973：所用的代码如下：
## Q1 oneTry <- function(n, k) { res <- sort(sample(1:n, size=k)) ifelse(1 %in% (res[-1] - res[-length(res)]), 1, 0) } ntry <- 10000000 nn <- 35 nk <- 5 set.seed(123) sum(replicate(ntry, oneTry(nn, nk))) / ntry # 0.4766973 when ntry = 10000000 and seed = 123.理论值的推导：（参考《Probability and Statistics (4ed)》 DEGROOT等著）
当 n ≥ 2 k − 1 n \ge 2k-1 n≥2k−1时，假设 i 1 < i 2 < ⋯ < i k i_1 < i_2 < \cdots < i_k i1<i2<⋯<ik是从1到 n n n这 n n n个数字中随机等概率且无重复抽取的一组数字，按从小到大排序。令 j s = i s − ( s − 1 ) , for s = 1 , 2 , … , k j_s = i_s - (s - 1), \text{ for } s = 1,2,\dots,k js=is−(s−1), for s=1,2,…,k，也就是： j 1 = i 1 , j 2 = i 2 − 1 , ⋮ j k = i k − ( k − 1 ) . \begin{aligned} j_1 & = i_1, \\ j_2 & = i_2 - 1, \\ & \ \ \vdots \\ j_k & = i_k - (k - 1). \end{aligned} j1j2jk=i1,=i2−1, ⋮=ik−(k−1). 那么可以证明：
( i 1 , i 2 , … , i k ) (i_1, i_2,\dots,i_k) (i1,i2,…,ik)这一组数中有连号的数等价于 ( j 1 , j 2 , … , j k ) (j_1, j_2,\dots,j_k) (j1,j2,…,jk)这一组数中有重复的数。对于 1 ≤ j 1 ≤ j 2 ≤ ⋯ ≤ j k ≤ n − k + 1 1 \le j_1 \le j_2 \le \cdots \le j_k \le n - k + 1 1≤j1≤j2≤⋯≤jk≤n−k+1， ( j 1 , j 2 , … , j k ) (j_1, j_2,\dots,j_k) (j1,j2,…,jk)中没有重复数的全部可能的组合的数量是 ( n − k + 1 k ) \binom{n-k+1}{k} (kn−k+1)。 ( i 1 , i 2 , … , i k ) (i_1, i_2,\dots,i_k) (i1,i2,…,ik)中没有连号数的全部可能的组合的数量是 ( n − k + 1 k ) \binom{n-k+1}{k} (kn−k+1)。 ( i 1 , i 2 , … , i k ) (i_1, i_2,\dots,i_k) (i1,i2,…,ik)中没有连号数的概率是 ( n − k + 1 k ) ( n k ) \cfrac {\binom{n-k+1}{k}} {\binom{n}{k}} (kn)(kn−k+1)。 ( i 1 , i 2 , … , i k ) (i_1, i_2,\dots,i_k) (i1,i2,…,ik)中至少有一组连号数的概率是 1 − ( n − k + 1 k ) ( n k ) 1 - \cfrac {\binom{n-k+1}{k}} {\binom{n}{k}} 1−(kn)(kn−k+1)。归纳法是数学中解决问题的重要方法，当然在概率统计里面也会大有用武之地。比如这个题目：
“假设有10个硬币，第 i i i个硬币正面朝上的概率是 1 2 i + 1 \frac{1}{2i+1} 2i+11，那么10个硬币都掷一次后，正面朝上的硬币数是偶数的概率是多少？”
我们可以把上面的题目抽象为：假设有 n n n（ n n n是正偶数）个硬币，第 i i i个硬币正面朝上的概率是 1 2 i + 1 \frac{1}{2i+1} 2i+11，那么 n n n个硬币都掷一次后，正面朝上的硬币数是偶数的概率是多少？
我们可以先看看 n = 2 n=2 n=2时的情形，可以很容易计算出结果是： 1 3 × 1 5 + 2 3 × 4 5 = 3 5 \frac{1}{3} \times \frac{1}{5} + \frac{2}{3} \times \frac{4}{5} = \frac{3}{5} 31×51+32×54=53
n = 4 n=4 n=4时的计算稍微复杂，结果是： 1 3 × 1 5 × 1 7 × 1 9 + 2 3 × 4 5 × 6 7 × 8 9 + ∑ 1 ≤ i < j ≤ 4 1 2 i + 1 × 1 2 j + 1 × ∏ k = 1 4 ( 1 − 1 2 k + 1 ) ( 1 − 1 2 i + 1 ) ( 1 − 1 2 i + 1 ) = 5 9 \begin{aligned} & \frac{1}{3} \times \frac{1}{5} \times \frac{1}{7} \times \frac{1}{9} \\ & + \frac{2}{3} \times \frac{4}{5} \times \frac{6}{7} \times \frac{8}{9} \\ & + \sum_{1 \le i < j \le 4} \frac{1}{2i+1} \times \frac{1}{2j+1} \times \frac {\prod_{k=1}^4 \left(1-\frac{1}{2k+1}\right) } {\left(1-\frac{1}{2i+1}\right) \left(1-\frac{1}{2i+1}\right)} \\ &= \frac{5}{9} \end{aligned} 31×51×71×91+32×54×76×98+1≤i<j≤4∑2i+11×2j+11×(1−2i+11)(1−2i+11)∏k=14(1−2k+11)=95
由 n = 2 n=2 n=2以及 n = 4 n=4 n=4时的结果，我们可以推测，更一般的结果是 n + 1 2 n + 1 \frac{n+1}{2n+1} 2n+1n+1。这可以用归纳法证明：
当 n = 2 n=2 n=2时结论成立； 假设当 n = k n=k n=k（ k k k是正偶数）时结论成立，即掷 n n n个硬币后正面朝上的硬币数为偶数的概率是 k + 1 2 k + 1 \frac{k+1}{2k+1} 2k+1k+1，那么当 n = k + 2 n=k+2 n=k+2时，相应的概率可以这样计算： Pr \Pr Pr( k + 2 k+2 k+2个硬币正面朝上的硬币数为偶数) = Pr \Pr Pr(前 k k k个硬币正面朝上的硬币数为偶数且后2个硬币正面朝上的硬币数为偶数) + Pr \Pr Pr(前 k k k个硬币正面朝上的硬币数为奇数且后2个硬币正面朝上的硬币数为奇数)。 并且， Pr \Pr Pr(前 k k k个硬币正面朝上的硬币数为偶数且后2个硬币正面朝上的硬币数为偶数) = k + 1 2 k + 1 × [ 1 2 ( k + 1 ) + 1 × 1 2 ( k + 2 ) + 1 + 2 ( k + 1 ) 2 ( k + 1 ) + 1 × 2 ( k + 2 ) 2 ( k + 2 ) + 1 ] = 2 k 2 + 5 k + 3 ( 2 k + 1 ) ( 2 k + 5 ) \frac{k+1}{2k+1} \times \left[\frac{1}{2(k+1)+1} \times \frac{1}{2(k+2)+1} + \frac{2(k+1)}{2(k+1)+1} \times \frac{2(k+2)}{2(k+2)+1} \right] = \frac{2k^2+5k+3}{(2k+1)(2k+5)} 2k+1k+1×[2(k+1)+11×2(k+2)+11+2(k+1)+12(k+1)×2(k+2)+12(k+2)]=(2k+1)(2k+5)2k2+5k+3； 同时， Pr \Pr Pr(前 k k k个硬币正面朝上的硬币数为奇数且后2个硬币正面朝上的硬币数为奇数) = [ 1 − k + 1 2 k + 1 ] × [ 1 2 ( k + 1 ) + 1 × 2 ( k + 2 ) 2 ( k + 2 ) + 1 + 2 ( k + 1 ) 2 ( k + 1 ) + 1 × 1 2 ( k + 2 ) + 1 ] = 2 k ( 2 k + 1 ) ( 2 k + 5 ) \left[1 - \frac{k+1}{2k+1} \right] \times \left[\frac{1}{2(k+1)+1} \times \frac{2(k+2)}{2(k+2)+1} + \frac{2(k+1)}{2(k+1)+1} \times \frac{1}{2(k+2)+1} \right] = \frac{2k}{(2k+1)(2k+5)} [1−2k+1k+1]×[2(k+1)+11×2(k+2)+12(k+2)+2(k+1)+12(k+1)×2(k+2)+11]=(2k+1)(2k+5)2k。 所以， Pr \Pr Pr( k + 2 k+2 k+2个硬币正面朝上的硬币数为偶数) = 2 k 2 + 5 k + 3 ( 2 k + 1 ) ( 2 k + 5 ) + 2 k ( 2 k + 1 ) ( 2 k + 5 ) = k + 3 2 k + 5 = ( k + 2 ) + 1 2 ( k + 2 ) + 1 = n + 1 2 n + 1 \frac{2k^2+5k+3}{(2k+1)(2k+5)} + \frac{2k}{(2k+1)(2k+5)} = \frac{k+3}{2k+5} = \frac{(k+2)+1}{2(k+2)+1} = \frac{n+1}{2n+1} (2k+1)(2k+5)2k2+5k+3+(2k+1)(2k+5)2k=2k+5k+3=2(k+2)+1(k+2)+1=2n+1n+1。 所以当 n = k + 2 n=k+2 n=k+2时结论也成立。
回到最开始的题目，当 n = 10 n=10 n=10时的概率是 11 21 \frac{11}{21} 2111。
我们用代码进行精确计算的值是0.5238095，模拟计算的值是0.523596。所用代码如下：
######### Q2 ############ ## Q2 - exact answer # get the probability that there are n coins facing up. getProbNUp <- function(nu) { # nu: num of coins facing up. if (nu == 0) { return(prod(down_prob)) } else if (nu == dice_num) { return(prod(up_prob)) } else if (nu < 0 | nu > dice_num) { print("Error: wrong param.") return(-1) } all_comb <- combn(dice_num, nu) vprob <- apply(all_comb, 2, function(x) prod(up_prob[x]) * prod(down_prob[-x])) sum(vprob) } dice_num <- 10 up_prob <- 1 / (2 * 1:dice_num + 1) down_prob <- 1 - up_prob sum(sapply(seq(0, dice_num, 2), getProbNUp)) # 0.5238095 ## Q2 - simulation oneTry <- function() { vres <- sapply(1:dice_num, function(i) # 0 - up; 1 - down; sample(0:1, size=1, replace=T, prob=c(up_prob[i], down_prob[i]))) ifelse(sum(vres == 0) %% 2 == 0, 1, 0) } ntry <- 1000000 set.seed(123) sum(replicate(ntry, oneTry())) / ntry # 0.523596 when ntry = 1000000 and seed = 123.这是一个基础的结论。我们可以用模拟其 p.d.f. \text{p.d.f.} p.d.f. 或者 c.d.f. \text{c.d.f.} c.d.f. 来看： 模拟 p.d.f. \text{p.d.f.} p.d.f.，用R语言中的 hist 或者ggplot2包中的geom_histogram 函数画出模拟的概率直方图。
模拟 c.d.f. \text{c.d.f.} c.d.f.，用R语言中的 ecdf 函数画出模拟的累积分布曲线。
具体代码如下：
obinom <- function(n, p) { sum(sample(0:1, size=n, replace = T, prob=c(1-p, p))) } sbinom <- function(n, p, N) { replicate(N, obinom(n, p)) } SEED <- 123 ber.p <- 0.3 ber.n <- 100 ber.N <- 100000 set.seed(SEED) ber.simu <- sbinom(ber.n, ber.p, ber.N) set.seed(SEED) ber.theo <- rbinom(ber.N, ber.n, ber.p) # p.d.f. (by plotting histogram and density) library(ggplot2) library(dplyr) library(tidyr) data.pdf <- tibble(Simulative=ber.simu, Theoretical=ber.theo) %>% gather(type, value) data.pdf %>% ggplot(aes(x=value, y=..density..)) + geom_histogram(aes(fill=type), position="identity", alpha=.3) + geom_density(aes(color=type), alpha=.3) # c.d.f. simu.cdf <- ecdf(ber.simu) theo.cdf <- ecdf(ber.theo) plot(simu.cdf, do.points=F, verticals=T, col="red", main="Simulation for sum of n i.i.d. variables\n(Bernoulli distribution)") lines(theo.cdf, do.points=F, verticals=T, col="blue") legend("bottomright", legend=c("Simulative", "Theoretical"), col=c("red", "blue"), lty=1, bty="n")