【R语言数据科学】（十二）：有趣的概率学（上）
- 🌸个人主页：JoJo的数据分析历险记
- 📝个人介绍：小编大四统计在读，目前保研到统计学top3高校继续攻读统计研究生
- 💌如果文章对你有帮助，欢迎✌
关注、👍点赞、✌收藏、👍订阅专栏- ✨本文收录于【R语言数据科学】本系列主要介绍R语言在数据科学领域的应用包括： R语言编程基础、R语言可视化、R语言进行数据操作、R语言建模、R语言机器学习算法实现、R语言统计理论方法实现。本系列会坚持完成下去，请大家多多关注点赞支持，一起学习~，尽量坚持每周持续更新，欢迎大家订阅交流学习！
文章目录
前言
你真的了解概率吗？在机会游戏中，概率有一个非常直观的定义。例如，我们知道一对骰子出现七的机会是六分之一。但是，现实情况中，概率并不是一个确定的值。在如今，我们常常喜欢用概率来解释现实问题。例如，明天下雨的概率、被鲨鱼咬的概率、癌症的概率等等。
知道如何计算概率会让你在机会游戏中占据优势，纵观历史，许多聪明的学者，包括著名的数学家，如Cardano、Fermat 和 Pascal，都花时间和精力思考这些游戏的数学。结果，概率论诞生了。概率在现代机会游戏中仍然非常有用。例如，在扑克中，我们可以根据桌上的牌来计算赢得一手牌的概率。此外，赌场依靠概率论来开发几乎可以肯定保证盈利的游戏。概率论在许多其他情况下很有用，特别是在依赖于以某种方式受偶然影响的数据的领域。这部分的所有其他章节都建立在概率论之上。因此，概率知识对于数据科学来说是必不可少的。
1.离散型变量概率
我们首先介绍与分类数据相关的一些基本概念。
我们使用P ( A ) P(A) P(A)定义A事件发生的概率。例如，一个抽屉里有五个球，3个蓝色的，2个红色的。现在我们从抽屉里面摸一个球。假设事件A为摸到的是红球，那么 P ( A ) = 0.4 P(A)=0.4 P(A)=0.4。在数据科学领域，我们经常要处理连续变量的情况。例如，工资大于10000的概率，则我们使用 P ( X > 10000 ) P(X>10000) P(X>10000)来表示，后续我们会继续讨论，现在我们先关注分类型变量
1.1 概率分布
如果我们知道不同类别的频率，定义类别结果的分布相对简单。我们只是为每个类别分配一个概率。如果我们从
44%的民主党人、
44%的共和党人、
10%的未决定选民和
2%的绿党人中随机选出可能的选民，这些比例决定了每个群体的概率。概率分布为：
Pr(民主党) = 0.44
Pr(共和党) = 0.44
Pr(未决定) = 0.10
Pr(绿党) = 0.02
1.2 蒙特卡洛模拟
计算机提供了一种实际执行上述简单随机实验的方法：例如，我们从包含三个蓝色珠子和两个红色珠子的袋子中随机挑选一个珠子。随机数生成器允许我们模拟随机挑选的过程。我们可以使用
R中的sample函数。首先，我们通过rep函数生成我们刚刚的珠子情况
beads <- rep(c('red','blue'), times = c(2,3))#两个红色，三个蓝色
beads
- 'red'
- 'red'
- 'blue'
- 'blue'
- 'blue'
下面我们使用sample函数进行抽样
sample(beads,1)
'red'
每次我们运行上述代码都产生一个随机结果。我们想无限次地重复这个实验。在实际情况中，我们往往是将实验重复足够多的次数，以使结果近似等同于无限重复。这是模特卡洛模拟的思想，关于具体的数学证明就不在这里描述了。在本节后面，我们提供了一种实用的方法来确定什么是“足够大”。为了演示模特卡洛模拟，我们使用replicate函数，可以让某一个任务重复n次
B <-1000
events <- replicate(B,sample(beads,1))
我们现在可以看看我们的定义是否真的与这个蒙特卡罗模拟近似一致。我们可以使用table查看分布：
tab <- table(events)
tab
events
blue red
599 401
使用prop.table得到概率
prop.table(tab)
events
blue red
0.599 0.401
可以看出结果基本和我们构造的总体分布一致，当然，B设置的越大，往往越接近真实值，因为我们的样本数越多。
3.独立性
如果一个事件的结果不影响另一个事件，我们说两个事件是独立的。典型的例子是抛硬币。每次我们投掷一枚公平的硬币，无论之前的投掷结果如何，看到正面的概率都是 1/2。当我们有放回的挑选珠子时也是如此。在上面的示例中，无论之前的抽签情况如何，出现红色的概率都是0.40。
然而，有的时候两个事件之间是紧密相关的。例如,有一副扑克，假设第一次抽到了国王，并且不放回，那么再抽到王的概率就是3/53
4.条件概率
当事件不独立时，条件概率很有用。我们已经看到了条件概率的一个例子：我们计算了第二张发牌是国王的概率，假设第一张是国王。在概率上，我们使用以下符号：
P(Card 2 is a king∣Card 1 is a king)=3/51。
让两个事件独立时，有P ( A ∣ B ) = P ( A ) P(A|B)=P(A) P(A∣B)=P(A)
5.加法原则和乘法原则
** 乘法原则： **
如果AB不独立
P
(
A
B
)
=
P
(
A
)
P
(
B
∣
A
)
P(AB) = P(A)P(B|A)
P(AB)=P(A)P(B∣A)
三个事件：
P
(
A
B
C
)
=
P
(
A
)
P
(
B
∣
A
)
P
(
C
∣
A
B
)
P(ABC)=P(A)P(B|A)P(C|AB)
P(ABC)=P(A)P(B∣A)P(C∣AB)
如果AB独立
P
(
A
B
)
=
P
(
A
)
P
(
B
)
P(AB)=P(A)P(B)
P(AB)=P(A)P(B)
** 加法原则: **
P
(
A
o
r
B
)
=
P
(
A
)
+
P
(
B
)
−
P
(
A
B
)
P(AorB)=P(A)+P(B)-P(AB)
P(AorB)=P(A)+P(B)−P(AB)
6.排列组合
在我们的第一个示例中，我们想象了一个带有五颗珠子的盒子。为了计算一次抽样的概率分布，我们简单地列出了所有的可能性。有 5 个，因此，对于每个事件，我们计算了这些可能性中有多少与该事件相关联。选择蓝色珠子的概率是 3/5，因为在五个可能的结果中，三个是蓝色的。
对于更复杂的情况，计算并不那么简单。例如，如果我抽出 5 张牌且没有放回，我得到所有花色相同的牌的概率是多少，这在扑克中称为“同花”，接下来我们介绍一些排列如何，方便我们使用 R 代码来计算答案。
首先，我们要构建一幅扑克牌，需要使用expand.grid和paste函数。paste函数是用来拼接字符串的,具体如下
name <-"JoJo"
year <-"21"
paste(name, year)
'JoJo 21'
函数 expand.grid 为我们提供了两个向量的所有组合。例如，如果我们有蓝色和黑色的裤子和白色、灰色和格子衬衫，那么所有组合是：
expand.grid(pants = c("blue","black"), shirt = c("white","grey","plaid"))
A data.frame: 6 × 2
pantsshirt<fct><fct>blue whiteblackwhiteblue grey blackgrey blue plaidblackplaid
有了这些基础，我们可以定义一副扑克
suits <- c("Diamonds","Clubs","Hearts","Spades")
numbers <- c("Ace","Deuce","Three","Four","Five","Six","Seven","Eight","Nine","Ten","Jack","Queen","King")
deck <- expand.grid(number=numbers, suit=suits)
deck <- paste(deck$number, deck$suit)
构建完这副牌后，我们可以通过计算满足我们条件的可能结果的比例来再次检查第一张牌中出现 King 的概率是 1/13
kings <- paste("King", suits)
mean(deck %in% kings)
0.0769230769230769
现在，假设第一张是
King
，第二张是
King
的条件概率是多少？之前，我们推论过，如果一张K已经出局并且还剩下
51
张，那么这个概率是
3/51
。让我们列出所有可能的结果来确认。这时我们需要使用 gtools 包中的 permutations 函数。对于任何大小为 n 的列表，此函数计算我们在选择 r 个项目时可以获得的所有不同组合。以下是我们可以从由 1、2、3 组成的列表中选择两个数字的所有方法
library(gtools)
permutations(3,2)
A matrix: 6 × 2 of type int
121321233132
# 创建从52张牌中，摸两张牌所有的情况
hands <- permutations(52,2, v = deck)
first_card <- hands[,1]#第一次摸到的牌的情况
second_card <- hands[,2]#第二次摸到的牌牌的情况
kings <- paste("King", suits)
sum(first_card%in%kings & second_card%in%kings)/ sum(first_card%in%kings)#第一次摸到king的情况下，第二次仍然摸到King的概率
0.0588235294117647
正如我们已经推断的那样，这恰好是
3/51
。
注意，上述我们介绍的情况和顺序是有关的，如果和顺序无关呢，我们使用组合数，combinations函数
combinations(3,2)
A matrix: 3 × 2 of type int
121323
7.三门问题
在 1970 年代，有一个名为“让我们做个交易”的游戏节目，Monty Hall 是主持人。在游戏的某个时刻，参赛者被要求从三个门中选择一个。在一扇门后面有一个奖品。其他门后面有一只山羊，以展示他们输掉的参赛者。参赛者选择一扇门后，在透露所选门是否有奖品之前，蒙蒂霍尔会打开剩下的两扇门中的一扇，向参赛者展示那扇门后面没有奖品。然后他会问
“你想换门吗？”你会怎么做？
我们可以用概率来证明，如果你坚持原来的门选择，你中奖的机会仍然是三分之一。但是，如果你切换到另一扇门，你的中奖机会是
三分之二！这似乎违反直觉。许多人错误地认为这两种机会都是
二分之一，因为您要在两个选项之间进行选择。下面我们使用蒙特卡罗模拟来看看哪种策略更好。具体的数学证明不在这里描述，接下来我们来看看如何编写代码。
B <-10000
monty_hall <-function(strategy){
doors <- as.character(1:3)#定义三个门
prize <- sample(c("car","goat","goat"))#三个门后面的东西
prize_door <- doors[prize =="car"]#有礼物的门
my_pick <- sample(doors,1)#第一次随机选择
show <- sample(doors[!doors %in% c(my_pick, prize_door)],1)#展示没有奖品的门
stick <- my_pick#不换的话之前选择的门
switch <- doors[!doors%in%c(my_pick, show)]#交换后选择的门
choice <- ifelse(strategy =="stick", stick, switch)#决定是否进行交换
choice == prize_door#是否得奖}#蒙特卡洛模拟
stick <- replicate(B, monty_hall("stick"))
print(c('不交换获奖的概率:',mean(stick)))#不交换的概率
switch <- replicate(B, monty_hall("switch"))
print(c('交换获奖的概率:',mean(switch)))#交换的概率
[1] "不交换获奖的概率:" "0.3257"
[1] "交换获奖的概率:" "0.6623"
可以看出，交换之后，我们的获奖概率为2/3! 不交换得到的概率是1/3
8.任意两个人同一天生日
假设你在一个有 50 人的教室里。如果我们假设这是一个随机选择的 50 人组，那么至少有两个人生日相同的概率是多少？这里我们使用蒙特卡罗模拟。为简单起见，我们假设没有人在 2 月 29 日出生。
首先，请注意生日可以表示为 1 到 365 之间的数字，因此可以获得 50 个生日的样本，如下所示：
n <-50
bdays <- sample(1:365,n)
上面我们就定义好了50个同学的生日，我们怎么判断是否是同一天呢？可以使用duplicated函数，它的作用是判断某个元素是否出现过，如果出现过返回TRUE。因此，我们只需要使用下列代码来判断是否存在至少有两个人同一天生日
any(duplicated(bdays))
FALSE
可以看出刚刚我们抽取的50个样本，不存在同一天生日的人！（
震惊！！！这种概率真的很小
）
为了估计组中同一天生日的概率，我们通过反复采样 50 个生日的集合来重复这个实验
B <-10000
same_birthday <-function(n){
bdays <- sample(1:365, n, replace=TRUE)
any(duplicated(bdays))}
results <- replicate(B, same_birthday(50))#重复10000遍操作
mean(results)
0.9681
大家可以看出概率
50
个人至少有两个人同一天生日的概率是
0.9681
，刚刚居然不存在。果然，小概率事件时常发生😂。在现实生活中，我们往往倾向于低估这些概率。为了直观地了解为什么它如此之高，想想当组大小接近 365 时会发生什么。在这个阶段，我们用完了所有天，概率是 1。假设我们想利用这些知识与朋友打赌一群人中有两个人生日相同。什么时候机会大于 50%？大于 75%？让我们创建一个查找表。我们可以快速创建一个函数来计算任何组大小：
compute_prob <-function(n,B=1000){
results <- replicate(B,same_birthday(n))
mean(results)}
使用函数 sapply，我们可以对任何函数执行逐元素操作
n <- seq(1,60)
prob <- sapply(n, compute_prob)
我们现在可以绘制一组不同人数情况下，至少存在两个人生日相同的概率估计图
library(tidyverse)
prob <- sapply(n, compute_prob)
qplot(n, prob)
上面是通过蒙特卡洛模拟得到的，其实这个概率我们可以通过概率论的知识直接测算出来。我们来计算不存在生日相同的概率。大家想，第一个人随便可以选一天，第二个人只能在剩下的364天选一天，以此类推。那么不存在相同生日的概率为：
1
×
364
365
×
363
365
.
.
.
365
−
n
+
1
365
1\times\frac{364}{365}\times\frac{363}{365}...\frac{365-n+1}{365}
1×365364×365363...365365−n+1
exact_prob <-function(n){
prob_unique <- seq(365,365-n+1)/3651- prod( prob_unique)#求连乘}
eprob <- sapply(n, exact_prob)
qplot(n, prob)+ geom_line(aes(n, eprob), col ="red")
该图表明，蒙特卡洛模拟提供了对准确概率的非常好的估计。如果不可能计算出准确的概率，我们仍然能够准确地估计概率。
9.如何选择蒙特卡洛重复试验的大小?
这里描述的理论需要不断地重复实验。在实践中，我们无法做到这一点。在上面的例子中，我们使用了B=10000次蒙特卡洛实验，结果证明这提供了准确的估计。这个数字越大，估计就越准确。但在更复杂的计算中，10000 可能还不够。此外，对于某些计算，10000 次实验在计算上可能不可行。在实践中，我们不知道答案是什么，所以我们不知道我们的蒙特卡洛估计是否准确。我们只知道，B越大 ，近似值越好。但是我们需要它有多大？这实际上是一个具有挑战性的问题。我们将在这里描述的一种实用方法是检查估计的稳定性。以下是一组25人的生日问题的示例。
B <-10^seq(1,5, len =100)
compute_prob <-function(B, n=25){
same_day <- replicate(B, same_birthday(n))
mean(same_day)}
prob <- sapply(B, compute_prob)
qplot(log10(B), prob, geom ="line")
在该图中，我们可以看到值在 1000 附近开始稳定。
本章的介绍到此介绍，谢谢阅读。如果文章对你有帮助，请多多点赞、收藏、评论、关注支持！！
版权归原作者 JoJo的数据分析历险记 所有， 如有侵权，请联系我们删除。