Python数学编程 第五章 集合与概率 第二节 概率

5.2 概率

我们可以使用集合对概率论中的基本概念进行解释,先从下面几个定义开始:

试验(experiment):简单来说就是我们要做的某种探索行为。我们要做实验是因为对实验的每个可能的结果的概率感兴趣。掷色子、掷硬币、从一副扑克中抽一张牌都可以看作是试验。试验的一次实施可称为一次尝试。

样本空间(sample space):一次实验所有可能的结果构成的集合称为样本空间。在公式中,通常用S表示样本空间。例如,要投掷一个6面的骰子,样本空间为:{1,2,3, 4, 5,6}。

事件(event):就是我们希望计算概率的那些试验结果的集合,即样本空间的子集。例如,我们可能希望知道某个特定结果的概率,如掷出一个3,或多个结果的集合的概率,如掷出一个偶数(2、4或6)。在公式中我们用字母E来表示一个事件。

若有一个均匀分布,即样本空间中的每个结果发生的可能性相同,那么一个事件发生的概率P(E)可由公式计算(稍后将介绍非均匀分布):


其中,n(E)和n(S)分别是事件集合E和样本空间S的基数。P(E)在0到1上取值,值越大表示事件发生的可能性越大。

我们将此公式应用到掷骰子的问题中来计算一个特定结果,比如说掷出3的概率:

S = {1, 2, 3,4, 5, 6}

E = {3}

n(S) = 6 n(E) = 1

P(E) = 1/6

这证明了一个明显的事实:掷骰子掷出一个特定结果的概率是1/6。你可以很容易地在头脑中进行计算,但我们也可以在Python中使用上述公式编写如下函数来计算样本空间(标签为space)中任意事件(标签为event)的概率:


 def probability(space, event):
     return len(event)/len(space)

在这个函数中,不需要使用FiniteSet()函数创建参数space和event,即样本空间和事件,他们也可以是列表或任意支持函数len()的Python对象。

使用此函数,可以编写一个程序来计算掷一个20面骰子时质数出现的概率:


 from sympy import FiniteSet
 def probability(space, event):
     return len(event)/len(space)
 
 def check_prime(number):                    # ①
     if number != 1:
         for factor in range(2 ,number):
             if number % factor == 0:
                 return False
     else:
         return False
     return True
 if __name__ == '__main__':
     space = FiniteSet(*range(1, 21))        # ②
     primes = []
     for num in space:
         if check_prime(num):                # ③
             primes.append(num)
     event = FiniteSet(*primes)              # ④
     p = probability(space, event)
     print('Sample space : {0}'.format(space))
     print('Event : {0}'.format(event))
     print('Probability of rolling a prime : {0}'.format(p))

首先在②处使用range()函数定义了一个样本空间的集合space。为创建事件集,需要从样本空间中找出质数,所以在①处定义了一个check_prime()函数,这个函数取一个整数并判断其是否被2和它自身自建的任意整数除尽(余数为0),如果存在这个整数,立即返回False。因为质数仅能被1和它自身整除,故如果一个整数是质数,则check_prime()函数返回True,否则给返回False。

我们在③处调用check_prime()函数对样本空间中的每个数进行判断,并把质数添加到一个列表primes中。然年后我们使用这个列表来创建一个事件集event。最后调用之前定义的probability()函数,传入样本空间和事件集。运行上述程序得到下述结果:


 Sample space : {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}
 Event : {2, 3, 5, 7, 11, 13, 17, 19}
 Probability of rolling a prime : 0.4

这里n(E) = 8, n(S) = 20,所以事件概率P = 0.4。

在20面骰子程序中,事实上不需要创建集合,只需要调用probability()函数,将样本空间和事件作为列表形式的参数传入即可。


 if __name__ == '__main__':
     space = range(1, 21)
     primes = []
     for num in space:
         if check_prime(num):
             primes.append(num)
     p = probability(space, primes)

在这种情况下,probability()函数仍然可以顺利执行。

5.2.1 事件A或事件B发生的概率

假设我们对两个可能发生的事件感兴趣,并希望计算至少有一个事件发生的概率。例如,回到骰子问题,考虑如下两个事件:

A = 出现质数

B = 出现奇数

如前所述,样本空间S = {1, 2, 3, 4, 5, 6}。事件A表示样本空间中的质数集合{2, 3, 5},事件B表示样本空间中的奇数集合{1,3,5}。为了计算至少有一个事件发生的概率,可以通过计算两个集合的并集的概率来实现。用符号可以表示为:


可以在Python中实现这个计算:


 >>> from sympy import FiniteSet
 >>> s = FiniteSet(1, 2, 3, 4, 5, 6)
 >>> a = FiniteSet(2, 3, 5)
 >>> b = FiniteSet(1, 3, 5)
 >>> e = a.union(b)
 >>> len(e)/len(s)
 0.6666666666666666

首先创建一个集合s,用来表示样本空间,接着创建集合a和b。然后使用union()函数得到事件集e。最后使用前面的公式计算两个集合的并集发生的概率。

5.2.2 事件A和事件B同时发生的概率

假设要计算两个事件同时发生的概率,例如,掷骰子出现的结果既是质数又是奇数的概率。这种情况下需要计算两个事件交集发生的概率:


我们可以通过intersect()函数计算事件A和事件B同时发生的概率,与前一个示例类似:


 >>> from sympy import FiniteSet
 >>> s = FiniteSet(1, 2, 3, 4, 5, 6)
 >>> a = FiniteSet(2, 3, 5)
 >>> b = FiniteSet(1, 3, 5)
 >>> e = a.intersect(b)
 >>> len(e)/len(s)
 0.3333333333333333

5.2.3 生成随机数

概率概念让我们对事情发生的可能性进行推理和推算,为了在实际中能使用计算机程序模拟类似于掷骰子游戏的事件,我们需要一种生成随机数的方法。

掷骰子模拟

为了模拟一个6面的骰子,我们需要一种方法来随机生成1到6之间的整数。Python的标准库中的random模块给我们能提供了各种函数来生成随机数,在这一章中我们将使用到两个函数:一个是randint()函数,它在给定范围内生成一个随机整数;另一个是random()函数,它生成一个0到1之间的浮点数。我们通过一个简单的例子来了解一下randint()是如何工作的:


 >>> import random
 >>> random.randint(1,6)
 4

randint()函数将两个整数作为参数传入,返回一个介于这两个整数之间(包含两端整数)的随机整数。在上面的程序中,我们把范围(1,6)传递给函数,它返回了数字4,但是我们再调用它,很可能得到一个不同的数字:


 >>> random.randint(1,6)
 3

调用randint()函数可以模拟真实的掷骰子游戏。每次调用这个函数,可以得到一个介于1和6之间的整数,就像我们只一个实际的6面骰子一样。注意:randint()函数要求你先提供比较小的数字,所以randint(6,1)是无效的。

能掷出某个数来吗?

下面的程序将模拟一个掷骰子游戏,在这个游戏中,我们将持续掷一个6面的骰子,直到掷出的总和超过20为止:


 '''
 Roll a die until the total score is 20
 '''
 import matplotlib.pyplot as plt
 import random
 target_score = 20
 def roll():
     return random.randint(1, 6)
 if __name__ == '__main__':
     score = 0
     num_rolls  = 0
     while score < target_score:
         die_roll = roll()
         num_rolls += 1
         print('Rolled: {0}'.format(die_roll))
         score += die_roll
     print('Score of {0} reached in {1} rolls'.format(score, num_rolls))

首先我们定义一个与之前创建过的相同函数roll()。然后我们使用while循环调用这个函数,掷出的结果相加,并且记录掷骰子的次数。重复循环直到总和超过20,这时程序输出总数以及掷骰子的次数。


 Rolled: 4
 Rolled: 1
 Rolled: 5
 Rolled: 6
 Rolled: 4
 Score of 20 reached in 5 rolls

运行程序,假如要实现总数超过20,你会发i西安所需投骰子的次数在不断变化。

目标总数是可能的吗?

我们的下一个程序是类似的,但它将告诉我们一个特定的目标总数是否可以在最多掷多少次投资的限定下达到:


 from sympy import FiniteSet
 import random
 def find_prob(target_score, max_rolls):
     die_sides = FiniteSet(1, 2, 3, 4, 5, 6)
     # Sample space
     s = die_sides**max_rolls            # ①
     # Find the event set
     if max_rolls > 1:
         success_rolls = []
         for elem in s:
             if sum(elem) >= target_score:
                 success_rolls.append(elem)
     else:
         if target_score > 6:
             success_rolls = []
         else:
             success_rolls = []
             for roll in die_sides:
                 if roll >= target_score:
                     success_rolls.append(roll)
     e = FiniteSet(*success_rolls)
     # Calculate the probability pf reaching target score
     return len(e)/len(s)
 
 if __name__ == '__main__':
     target_score = int(input('Enter the target score: '))
     max_rolls = int(input('Enter the maximum number of rolls allowec: '))
     p = find_prob(target_score, max_rolls)
     print('Probability : {0:.5f}'.format(p))

运行这个程序,他首先询问目标总数和允许的最大投掷次数病机那个他们作为输入参数,然后输出能完成这一任务的概率。

在①处,对基本样本空间{1,2,3,4,5,6}做笛卡尔积,次幂为max_rolls,得到投掷max_rolls次所有可能的骰子点数组合。如果投掷次数大于1,则在生成的笛卡尔积中对于每项加和与target_score做对比,满足条件,就将该项笛卡尔积加入到success_rolls列表中,作为满足完成任务的事件;如果投掷次数<=1,且target_score>6,显然这是不可能的,故success_rolls为一个空列表;假如投掷次数为且目标分数<=6时,此时die_sides就相当于{1,2,3,4 ,5,6},此时,只需进行一次遍历即可查找到满足条件的事件。

执行程序:


 Enter the target score: 7
 Enter the maximum number of rolls allowec: 1
 Probability : 0.00000

投掷一次不可能得到7点,故概率为0。


 Enter the target score: 11
 Enter the maximum number of rolls allowec: 2
 {(5, 6), (6, 5), (6, 6)}
 Probability : 0.08333

投掷两次且目标分数为11,满足该条件的事件为{(5, 6), (6, 5), (6, 6)},可计算出概率为:


其中6^2为笛卡尔积的基数,也就是投掷两次骰子1所有可能出现的组合个数。

5.2.4 非均匀随机数

目前为止,我们对概率的讨论都假设了样本空间中的每个结果发生的可能性相同。例如,random.randint()函数返回指定范围内的一个整数,并假设每个整数返回的可能性相同,我们称这种概率为均匀概率(uniform probability),并将randint()函数生成的随机数称为均匀随机数。然而,我们要模拟投掷一个不均匀的硬币,这个硬币正面出现的概率是反面出现概率的两倍,这时我们就需要一种生成非均匀随机数的方法了。

在编写程序实现该方法之前,我们先来看一下生成非均匀随机数的思路。

考虑数轴上的0到1区间,等分该区间,如下图所示。

我们视这条线为概率数轴,每一部分表示一个等可能性的结果,好比硬币的正面或者反面。现在,下图的数轴表示了一个不同的版本。

这里,正面对应的部分长度为2/3,反面对应的部分长度为1/3,这代表正面出现的概率是反面的两倍,即正反面出现的概率不相等。下面的Python函数模拟量这样掷硬币:


 import random
 def toss():
     # 0 -> Heads, 1-> Tails
     if random.random() > 2/3
         return 0
     else:
         return 1

假设函数返回0表示反面,返回1表示正面,那么我们可以使用random.randint()函数来返回0到1之间的随机数。如果生成的随机数小于等于2/3,就返回1,否则返回0。

现在我们来考虑如何推断之前的函数来模拟具有多个可能结果的非均匀事件。设想我们有一个xunideATM机,当按下按钮时,可以分发出一张5美元、10美元、20美元或者50美元的钞票。不同的面额有着不同的分发概率,如下图所示。

这里,面额为5美元或10美元的钞票被分发的概率都是1/6,面额为20美元或50美元被分发的概率为1/3。

我们创建一个列表来存储概率的波动和,然后生成一个介于0到1之间的随机数。从存储概率的列表的左端开始,返回该列表的第一个索引,其对应的和小于或等于生成的随机数。get_index()函数实现了这一想法。


 '''
 Simulate a fictional ATM that dispensescdollar bills
 of various denomination with varying probability
 '''
 
 import random
 def get_index(probability):
     c_probability = 0
     sum_probability = []                # ①
     for p in probability:
         c_probability += p
         sum_probability.append(c_probability)
     r = random.random()                 # ②
     for index, sp in enumerate(sum_probability):
         if r <= sp:                     # ③
             return index
     return len(probability) - 1         # ④
 
 def dispense():
     dollar_bills = [5, 10, 20, 50]
     probability = [1/6, 1/6, 1/3, 1/3]
     bill_index = get_index(probability)
     return dollar_bills[bill_index]

我们将一个列表作为输入参数来调用get_index()函数,该列表包含了相应位置的事件发生的概率。在①处,我们构造了列表sum_probability,其中第i个元素是列表probability的第i个元素之和。也就是说,sum_probability的第一个元素等于probability的第一个元素,第二个元素等于probability的前两个元素之和,以此类推。在②处,用标签r指代生成的一个介于0到1之间的随机数。接下来,我们遍历sum_probability并返回第一个超过r的元素的索引。

函数的最后一行(④处),处理了一种情形,最好通过一个例子来说明其原理。考虑有三个事件的列表,每个事件发生的概率百分比为0.33.这种情形下,列表sum_probability可表达成[0.33,0.66,0.99]。现在假设生成的随机数r是0.99314.对这个r值,我们希望选择事件列表中的最后一个元素,你可能会说,这并不正确,因为这样一来最后一个是事件被选中的可能性大于0.33。按照③处的条件,sum_probability中没有大于r的元素,因此函数根本不会返回任何索引。④处语句处理了这种情形,并返回最后一个索引。

如果调用dispenses()函数模拟ATM机分发的大量钞票,你进啊感到每种钞票出现的次数的比例严格遵守其指定的概率。下一章我们将看到在创建分分形时这一技巧非常有用。

发表评论
留言与评论(共有 0 条评论) “”
   
验证码:

相关文章

推荐文章