R-概率统计与模拟（五）彩票连号、归纳法以及二项分布

《R-概率统计与模拟》的系列文章已经写了四篇。其中有经过整理，内容相对全面，并且有内在联系的文章；也有文内几个小节相对独立，没有内在联系的文章。后一类文章正是笔者在学习过程中当时当地遇到了有意思的问题，即时地去解决并记录下来的结果。本文也仍然是这一类风格。学习是零散的过程，等到积累了一定的素材，才能进行整理归纳。

本文是《R-概率统计与模拟》系列文章的第五篇，包含了三个小节：

• 彩票至少包含一组连续号码的概率
• 用归纳法解决概率问题的一个例题
• 多个独立且符合同一个伯努利分布的变量的和服从二项分布

彩票至少包含一组连续号码的概率

买过超级大乐透的朋友应该熟悉，这种颇受欢迎的彩票分为两区，前区是从1到35这35个数字中无重复地选择5个数字，后区是从1到12这12个数字中无重复地选择2个数字。选定的前（后）区数字与开奖的前（后）区数字重合得越多，奖金越多。

很多彩票分析网站在统计历史开奖数字的各种信息，希望能从中找出些许的规律。这些信息中很多都是在统计“频率”，而这些“频率”其实都是有理论的概率值的。比如，“前区号码中有连号的概率是多少”？也就是：从1到35这35个数字中随机等概率且无重复地抽取5个数字，这5个数字中至少有两个数字是连号的概率是多少？

这个问题可以抽象为：从1到$n$这$n$个数字中随机等概率且无重复地抽取$k(\text{for }n\ge 2k-1)$个数字，这$k$个数字中至少有两个数字是连号的概率是多少？
（可以简单推论，当$n<2k-1$时，抽取出的$k$个数字中有连号数字的概率是1。）

这是一个蛮有意思的概率题，我们可以计算理论值，也可以用模拟的方法去计算其近似值。

其理论值是：
$$1-\frac{\left(\genfrac{}{}{0ex}{}{n-k+1}{k}\right)}{\left(\genfrac{}{}{0ex}{}{n}{k}\right)}$$

对于超级大乐透前区，$n=35,k=5$，所以连号的概率是：
$$1-\frac{\left(\genfrac{}{}{0ex}{}{35-5+1}{5}\right)}{\left(\genfrac{}{}{0ex}{}{35}{5}\right)}\approx 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\ge 2k-1$时，假设$i_1<i_2<\cdots <i_k$是从1到$n$这$n$个数字中随机等概率且无重复抽取的一组数字，按从小到大排序。令$j_s=i_s-(s-1),\text{ for }s=1,2,\dots ,k$，也就是：
$$\begin{array}{rl}{j}_1& =i_1,\\ j_2& =i_2-1,\\ & ⋮\\ j_k& =i_k-(k-1).\end{array}$$
那么可以证明：

1. $(i_1,i_2,\dots ,i_k)$这一组数中有连号的数等价于$(j_1,j_2,\dots ,j_k)$这一组数中有重复的数。
2. 对于$1\le j_1\le j_2\le \cdots \le j_k\le n-k+1$，$(j_1,j_2,\dots ,j_k)$中没有重复数的全部可能的组合的数量是$\left(\genfrac{}{}{0ex}{}{n-k+1}{k}\right)$。
3. $(i_1,i_2,\dots ,i_k)$中没有连号数的全部可能的组合的数量是$\left(\genfrac{}{}{0ex}{}{n-k+1}{k}\right)$。
4. $(i_1,i_2,\dots ,i_k)$中没有连号数的概率是$\frac{\left(\genfrac{}{}{0ex}{}{n-k+1}{k}\right)}{\left(\genfrac{}{}{0ex}{}{n}{k}\right)}$。
5. $(i_1,i_2,\dots ,i_k)$中至少有一组连号数的概率是$1-\frac{\left(\genfrac{}{}{0ex}{}{n-k+1}{k}\right)}{\left(\genfrac{}{}{0ex}{}{n}{k}\right)}$。

用归纳法解决概率问题的一个例题

归纳法是数学中解决问题的重要方法，当然在概率统计里面也会大有用武之地。比如这个题目：

“假设有10个硬币，第$i$个硬币正面朝上的概率是$\frac{1}{2i+1}$，那么10个硬币都掷一次后，正面朝上的硬币数是偶数的概率是多少？”

我们可以把上面的题目抽象为：假设有$n$（$n$是正偶数）个硬币，第$i$个硬币正面朝上的概率是$\frac{1}{2i+1}$，那么$n$个硬币都掷一次后，正面朝上的硬币数是偶数的概率是多少？

我们可以先看看$n=2$时的情形，可以很容易计算出结果是：
$$\frac{1}{3}\times \frac{1}{5}+\frac{2}{3}\times \frac{4}{5}=\frac{3}{5}$$

$n=4$时的计算稍微复杂，结果是：
$$\begin{array}{rl}& \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{array}$$

由$n=2$以及$n=4$时的结果，我们可以推测，更一般的结果是$\frac{n+1}{2n+1}$。这可以用归纳法证明：

当$n=2$时结论成立；
假设当$n=k$（$k$是正偶数）时结论成立，即掷$n$个硬币后正面朝上的硬币数为偶数的概率是$\frac{k+1}{2k+1}$，那么当$n=k+2$时，相应的概率可以这样计算：
$\mathrm{Pr}$($k+2$个硬币正面朝上的硬币数为偶数) = $\mathrm{Pr}$(前$k$个硬币正面朝上的硬币数为偶数且后2个硬币正面朝上的硬币数为偶数) + $\mathrm{Pr}$(前$k$个硬币正面朝上的硬币数为奇数且后2个硬币正面朝上的硬币数为奇数)。
并且，
$\mathrm{Pr}$(前$k$个硬币正面朝上的硬币数为偶数且后2个硬币正面朝上的硬币数为偶数) = $\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)}$；
同时，
$\mathrm{Pr}$(前$k$个硬币正面朝上的硬币数为奇数且后2个硬币正面朝上的硬币数为奇数) = $\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)}$。
所以，
$\mathrm{Pr}$($k+2$个硬币正面朝上的硬币数为偶数) = $\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}$。
所以当$n=k+2$时结论也成立。

回到最开始的题目，当$n=10$时的概率是$\frac{11}{21}$。

我们用代码进行精确计算的值是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.

多个独立且符合同一个伯努利分布的变量的和服从二项分布

这是一个基础的结论。我们可以用模拟其 $\text{p.d.f.}$ 或者 $\text{c.d.f.}$ 来看：
模拟$\text{p.d.f.}$，用R语言中的 hist 或者ggplot2包中的geom_histogram 函数画出模拟的概率直方图。

模拟$\text{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")