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

作者: 生信了 | 来源:发表于2019-12-24 16:47 被阅读0次

R——概率统计与模拟（五）彩票连号、归纳法以及二项分布
R——概率统计与模拟（二）
R统计学(06): 负二项分布
R统计学(05): 泊松分布
R——概率统计与模拟（一）
R——概率统计与模拟（四）拒绝抽样
数据探索之统计分布
统计学学习-1
第一章算法基础——概率论与数理统计基础
R语言学习笔记_03

原创：hxj7

《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 - \cfrac {\binom{n-k+1}{k}} {\binom{n}{k}}$

对于超级大乐透前区， $n=35, k=5$ ，所以连号的概率是：
$1 - \cfrac {\binom{35-5+1}{5}} {\binom{35}{5}} \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{aligned} j_1 & = i_1, \\ j_2 & = i_2 - 1, \\ & \ \ \vdots \\ j_k & = i_k - (k - 1). \end{aligned}$
那么可以证明：

$(i_1, i_2,\dots,i_k)$ 这一组数中有连号的数等价于 $(j_1, j_2,\dots,j_k)$ 这一组数中有重复的数。

对于 $1 \le j_1 \le j_2 \le \cdots \le j_k \le n - k + 1$ ， $(j_1, j_2,\dots,j_k)$ 中没有重复数的全部可能的组合的数量是 $\binom{n-k+1}{k}$ 。

$(i_1, i_2,\dots,i_k)$ 中没有连号数的全部可能的组合的数量是 $\binom{n-k+1}{k}$ 。

$(i_1, i_2,\dots,i_k)$ 中没有连号数的概率是 $\cfrac {\binom{n-k+1}{k}} {\binom{n}{k}}$ 。

$(i_1, i_2,\dots,i_k)$ 中至少有一组连号数的概率是 $1 - \cfrac {\binom{n-k+1}{k}} {\binom{n}{k}}$ 。

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

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

“假设有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{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}$

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

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

（公众号：生信了）

R——概率统计与模拟（五）彩票连号、归纳法以及二项分布
原创：hxj7 《R-概率统计与模拟》的系列文章已经写了四篇。其中有经过整理，内容相对全面，并且有内在联系的文章；...
R——概率统计与模拟（二）
原创：hxj7 本文继续介绍一些和概率统计相关的模拟。前文《R-概率统计与模拟》介绍了一些用 R 进行概率模拟的...
R统计学(06): 负二项分布
前面我们介绍了多种离散型概率分布，大家可以点击下方链接来回顾： R统计学(01): 伯努利分布、二项分布 R统计学...
R统计学(05): 泊松分布
前面我们介绍了多种离散型概率分布，大家可以点击下方链接来回顾： R统计学(01): 伯努利分布、二项分布 R统计学...
R——概率统计与模拟（一）
原创：hxj7 本文记录了三个概率统计相关的小题目，以回顾一些概率统计的知识。正如笔者在前文《公众号一岁啦》中所...
R——概率统计与模拟（四）拒绝抽样
原创：hxj7 本文介绍了拒绝抽样（Reject Sampling）。前文《R-概率统计与模拟（三）变换均匀分布...
数据探索之统计分布
本文用Python统计模拟的方法，介绍四种常用的统计分布，包括离散分布：二项分布和泊松分布，以及连续分布：指数分布...
统计学学习-1
看了几篇统计学资料：恍然间不知道概率分布[二项分布] 与抽样分布[t 分布] 差别。大家知道，统计学分为描...
第一章算法基础——概率论与数理统计基础
1.3 概率论与数理统计基础算法常涉及数据分布情况，而这些分布又与概率紧密相连，常见的分布方式包括二项分布、超几...
R语言学习笔记_03
统计分布和描述性统计分析概率函数通常用来生成特征已知的模拟数据，以及在用户编写的统计函数中计算概率值。 d = ...