抽屉里有红袜子和黑袜子,随机抽出两只,两只全红的概率为 1/2。
问:
a. 抽屉里的袜子最少有多少只?
b. 如果黑袜子为偶数只,则抽屉里的袜子最少有多少只?
思路参考
a
设抽屉里有 x 只红袜子,y 只黑袜子,总共 N = x + y 只袜子,记随机抽两只均为红色的概率为 P
问题变成解以上不定方程 N > 2 的最小正整数解
解析解
由于对于 y > 0, x > 1,有以下不等式
再结合
我们可以得到不等式
从左半部分不等式可以计算出
从右半部分不等式可以计算出
因此对于 y = 1, x 的范围在,所以 x = 3 是一个候选答案
数值解
暴力解法代码如下, 从 y = 1, x = 2 开始枚举,判断是否满足方程。
y = 1
found = False
while True:
x = 2
while True:
N = x + y
if 2 * x * (x - 1) == N * (N - 1):
print("N: {}, x: {}".format(x, N))
found = True
break
elif 2 * x * (x - 1) > N * (N - 1):
break
x += 1
if found:
break
y += 1
得到 N = 4, x = 3
b
解析解
利用在此前已经推出的以下公式
这里规定了 y 是偶数。依次考虑 y = 2, 4, 6, ...,对每一个 y,我们可以知道 x 的范围,这个范围的长度正好为 1,因此也就是可以知道对应的 x 具体是多少。考察这对 (x, y) 是否满足 P = 1/2 即可。
第一个满足的就是答案。
数值解
依然用暴力法解,只需要把 y 的初始值改为 2, 每轮 y 的增加量改为 2 即可
y = 2
found = False
while True:
x = 2
while True:
N = x + y
if 2 * x * (x - 1) == N * (N - 1):
print("N: {}, x: {}".format(x, N))
found = True
break
elif 2 * x * (x - 1) > N * (N - 1):
break
x += 1
if found:
break
y += 2
得到 N = 21, x = 15
有一个 3 人组成的陪审团,其中两个人独立做决定均有 p 概率做对,另一个人通过抛硬币做决定
还有一个 1 人的陪审团,那个人做决定也有 p 概率做对。
问:哪个陪审团更优可能做出正确决定?
思路参考
第一个陪审团做出正确决定的概率记为 P1, 第二个陪审团做出正确决定的概率记为 P2。
因此 P1 = P2,两个陪审团做出正确决定的可能性相等。
构造数据验证
T = 10
n_trails = int(1e6)
for t in range(T):
p = np.random.random()
# 模拟第一个陪审团 3 个人各自做出的判断,1为正确判断
# 模拟 n_trails 次
c1 = np.random.binomial(1, p, size=n_trails)
c2 = np.random.binomial(1, p, size=n_trails)
c3 = np.random.binomial(1, 0.5, size=n_trails)
# c 中记录了 n_trails 次模拟中每次是否做出了正确判断
c = ((c1 + c2 + c3) >= 2).astype(int)
print("P2 = {:.5f}, P1 = {:.5f}".format(p, np.mean(c)))
实验结果
P2 = 0.34158, P1 = 0.34095
P2 = 0.55780, P1 = 0.55776
P2 = 0.07460, P1 = 0.07439
P2 = 0.01918, P1 = 0.01919
P2 = 0.39175, P1 = 0.39233
P2 = 0.82699, P1 = 0.82643
P2 = 0.84726, P1 = 0.84769
P2 = 0.00396, P1 = 0.00395
P2 = 0.09100, P1 = 0.09023
P2 = 0.45738, P1 = 0.45741
Elmer 如果在三场的系列赛中赢下连续的两场,则可以得到奖励,系列赛的对手安排有两种
这里冠军的水平高于爸爸,问:Elmer 应该选择哪一组,得到奖励的概率更大?
思路参考
一方面,由于冠军水平比爸爸高,应该尽可能少与冠军比赛,另一方面,中间的那一场是更重要的一场,因为如果中间那场输了就不可能连赢两场。
Elmer 与爸爸比赛获胜概率为 P1,与冠军比赛获胜概率为 P2,由于冠军水平高于爸爸,所以 P1 > P2。
对于一组系列赛,Elmer 要得到奖励有两种可能性,一种是前两场都赢了,此时不用看第三场,另一种是第一场输了,但第二第三场都赢了。
下面分别考虑两组系列赛,计算得到奖励的概率
第一组系列赛,前两场都赢的概率为 P1 * P2
,第一场输但是第二、三场都赢的概率为 (1 - P1) * P2 * P1
,得到奖励概率记为 Pa
Pa=P1×P2+(1−P1)×P2×P1
第二组系列赛,前两场都赢的概率为 P2 * P1
,第一场输但是第二、三场都赢的概率为 (1 - P2) * P1 * P2
,得到奖励概率记为 Pb
Pb=P2×P1+(1−P2)×P1×P2
做差比零
Pa - Pb = (P1 \times P2 + (1 - P1) \times P2 \times P1) - (P2 \times P1 + (1 - P2) \times P1 \times P2) \\ = ((1 - P1) - (1 - P2)) \times P1 \times P2 \\ = (P2 - P1) \times P1 \times P2 \\ &< 0 \\
所以 Pa < Pb,Elmer 应该选择第二组系列赛。
构造数据验证
import numpy as np
step = 0.2
for P1 in np.arange(step, 1 + step, step):
for P2 in np.arange(step, P1 - step, step):
print("P1: {:.2f}, P2: {:.2f}".format(P1, P2))
Pa = P1 * P2 + (1 - P1) * P2 * P1
Pb = P2 * P1 + (1 - P2) * P1 * P2
print("Pa: {:.2f}, Pb: {:.2f}".format(Pa, Pb))
print("----")
结果如下:可看到 Pb 的值始终比 Pa 大
P1: 0.60, P2: 0.20
Pa: 0.17, Pb: 0.22
----
P1: 0.60, P2: 0.40
Pa: 0.34, Pb: 0.38
----
P1: 0.80, P2: 0.20
Pa: 0.19, Pb: 0.29
----
P1: 0.80, P2: 0.40
Pa: 0.38, Pb: 0.51
----
P1: 0.80, P2: 0.60
Pa: 0.58, Pb: 0.67
----
P1: 1.00, P2: 0.20
Pa: 0.20, Pb: 0.36
----
P1: 1.00, P2: 0.40
Pa: 0.40, Pb: 0.64
----
P1: 1.00, P2: 0.60
Pa: 0.60, Pb: 0.84
----
P1: 1.00, P2: 0.80
Pa: 0.80, Pb: 0.96
----
一枚骰子掷出第一个 6 平均需要掷多少次。
思路参考 1
记 p(k) := 第一个 6 需要掷 k 次的概率。p := 在一次投掷中得到 6 的概率;q := 在一次投掷中不是 6 的概率,由定义 q = 1 - p。
k = 1: p(1) = p
k = 2: p(2) = pq
k = 3: p(3) = pq^2
...
总结规律后可以得到
p(k)=pq k−1
数学期望为
处理上面的级数的技巧如下,首先两边都乘以 q。
让两式相减
其中几何级数那部分的公式如下
将上面的公式代入到两式相减的公式中
将 p = 1/6 代入,E(k) = 6
思路参考 2
这里我们仍沿用思路参考 1 的数学记号。
记 m 为第一个 6 平均需要的次数(即期望)。
由全期望公式【随机变量 k 的期望就是各种条件下 k 的期望加权求和】。这里条件期望有两种情况:1)当第一次投掷成功时,则需要的平均次数是 1(该情况下期望是 1,对应权重即发生概率为 p);2)当第一次投掷失败时,则需要的平均次数是 1 + m (该情况下期望是 1 + m,权重为 q),可列出如下等式:
蒙特卡洛模拟验证结果
import numpy as np
p = 1/6
def test():
T = int(1e6)
mapping = {}
for t in range(T):
i = 1
while True:
if np.random.rand() < p:
break
i += 1
if i not in mapping:
mapping[i] = 0
mapping[i] += 1
E = 0
for k in mapping:
E += k * (mapping[k] / T)
print("Average times is {:.4f}".format(E))
for i in range(5):
test()
实验结果
Average times is 6.0025
Average times is 5.9918
Average times is 6.0018
Average times is 6.0038
Average times is 6.0059
回顾在之前 试验直到第一次成功 中解决的一个问题,问题和解法如下
一枚骰子掷出第一个 6 平均需要掷多少次。
当时我们用两种方法解决了这个问题:
期望的定义
第一种是直接用期望的定义求期望
由于 i 的所有可能取值和 p(i) 我们都可以算出来,于是通过定义再加上一些公式推导技巧我们求出了 E(X) 的解析解。
全期望公式
第二种是基于 全期望公式 求期望
当时在那篇文章中的思路是这么写的
记 m 为第一个 6 需要的次数
当第一次投掷成功时,则需要的平均次数是 1
当第一次投掷失败时,则需要的平均次数是 1 + m
于是有以下公式,进而可以求出 m
m = p * 1 + q * (1 + m)
从全期望公式的角度看
X 是掷出第一个 6 的次数,m 就是 E(X);
Y 是第一次投掷出的点数是否为 6,是则记 Y=1,否则记 Y=0;
Y=1 的概率为 p,Y=0 的概率为 q;
E(X|Y=1) = 1,E(X|Y=0) = 1 + E(X)
代入公式即可得到解题思路中用的公式 m = p * 1 + q * (1 + m)
E(X)=p(Y=1)E(X∣Y=1)+p(Y=0)E(X∣Y=0)
=p×1+q×(1+E(X))
有向图
全期望公式的应用难点在于找出 E(X|y)。上述问题中 E[X| y=0] = (1 + E(X)), E[X| y=1] = 1 比较简单,我们能快速看出 E[X|y],而在稍微复杂的问题上,我们需要有一种找到 E[X|y] 的方法论 -- 有向图。
用全期望公式求期望的过程可以用有向图的方式表示出来。以上述问题为例,我们可以画出下面的有向图。
其中节点 0 表示未投掷出 6 的状态,节点 1 为表示掷出 6 的状态。0 为起点,1 为终点。
图中的每条有向边有一个概率值,还有另一个权值,节点有一个期望值。
我们要求的是从起始状态走向终点状态要走的步数的期望,每个节点有一个期望值表示从该节点走到终点节点的期望。由于 1 是终点,所以 E[1] = 0。0 是起点,我们要求的就是 E[0]。
从节点 0 会伸出两条边,一条回到 0,表示当前的投掷没有掷出 6,另一条通向 1,表示当前的投掷掷出了 6。边上的概率值表示当前投掷的各个值的概率,而投掷一次步数加一,因此额外的边权都是 1。
定义好这个图之后,我们就可以写出 E[0] 的表达式了
E[0] = p * (1 + E[1]) + q * (1 + E[0])
这个表达式就是全期望公式,只是将全期望公式中 E(X|y) 的部分用图中点权和边权的方式定义了计算方法。
期望 DP
通过以上的分析可以发现,当根据具体问题定义好有向图 D,起点 s,终点 t 后,我们就可以用类似于动态规划的方式求从 s 到 t 的某个指标的期望。
状态定义
dp[u] := 从 u 到 t 的某个指标的期望
初始化
dp[t]
答案
dp[s]
状态转移
dp[u] = g(node[u], sum(p[u][v] * f(w[u][v], dp[v])))
其中 f 是边权 w[u][v] (在上面的题中是 1) 和下一点的期望 dp[v] 的函数,g 是当前点的点权 node[u](在上面的题中没有)和 sum(p[u][v] * f) 的函数, 具体的函数形式取决于问题(在上面的题中 f 就是简单的 w[u][v] + dp[v], g 没有),需要具体分析。
求解 dp[s] 的过程:对反图(将原图中的每条有向边的方向反转之后形成的图)进行拓扑排序,按照拓扑排序的顺序依次求解各个节点 u 的 dp[u],直至求解到 dp[s]。
高斯消元
当建出的有向图中有环的时候,求解 E[s] 的过程如果直接用期望 DP 从 dp[t] 向 dp[s] 推导的话是不行的,因为在推导 DP 状态时会出现类似于下面的情况(就是后面的【用户一次游戏带来的收入】那道题的转移方程,这里提前看一下)
dp[0] = dp[1]
dp[1] = 1 + 0.7 * dp[2] + 0.3 * dp[0]
dp[2] = 2 + 0.6 * dp[3] + 0.4 * dp[1]
dp[3] = 3 + 0.5 * dp[4] + 0.5 * dp[2]
dp[4] = 0
这种情况 dp 状态的转移形成了环,比如要求 dp[1] 要先知道 dp[2],要求 dp[2] 就要先知道 dp[1],没法求了。
如果方程组是线性方程组,还有有办法的,解决办法是利用高斯消元的方式整体解这个线性方程组。
关于高斯消元的推导、代码、题目,可以找相关资料学习。
方法总结
现总结求解这类期望问题的方法论。
1、理论基础:全期望公式;
2、处理具体问题时,我们先分析建图所需的信息:
1)起点 s 和终点 t 状态是什么,还有哪些状态节点;
2)求出各个状态节点一步可以走到哪些状态节点,概率又分别是多少;
3)分析目标期望是否需要额外的点权和边权。
3、建图,初始化 dp[t];
4、分析状态转移方程中的 g,f 并写出实际的转移方程;
5、求解 dp[s],这里注意两点:
题目: 装备升级
下面开始今天要看的题目。题目描述如下:
玩家的装备有从 0 到 4 这 5 个等级,每次升级都需要一个工序,一共有 0->1, 1->2, 2->3, 3->4 这 4 个工序,成功率分别为 0.9, 0.7, 0.6, 0.5。工序成功,玩家的装备就会升一级,但如果失败,玩家的装备就要倒退一级。例如玩家当前的等级为 2,目前执行 2->3 这个工序,如果成功,则装备等级从 2 升为 3,如果失败,装备等级就从 2 降到 1。
问:玩家装备初始等级为 0, 升到 5 平均要经历多少个工序。
思路参考
步骤 1: 建图
按照前面的方法总结。我们首先分析题目并建图,点状态是 0,终点状态是 4,此外还有 1, 2, 3 这三个中间状态。题目描述中明确给出了状态间的转移方向以及概率。我们要求的是从 s=0 走到 t=4 平均需要走的步数。
由于每经过一条边都相当于走了一步,所以边上有额外的权 1。除此之外没有别的权,节点上也没有权。根据这些信息,我们首先把图建出来,如下
步骤 2:期望DP
期望 DP 的方程组如下:
dp[0] = 0.9 * (1 + dp[1]) + 0.1 * (1 + dp[0])
dp[1] = 0.7 * (1 + dp[2]) + 0.3 * (1 + dp[0])
dp[2] = 0.6 * (1 + dp[3]) + 0.4 * (1 + dp[1])
dp[3] = 0.5 * (1 + dp[4]) + 0.5 * (1 + dp[2])
dp[4] = 0
这是一个有环的图,且方程是线性方程组,面试的时候手推就可以。
如果要编程的话,需要整理成标准形式后用高斯消元解决。
步骤3 :高斯消元
标准形式如下
高斯消元求解
#include <vector>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
int main()
{
// 系数矩阵
vector<vector<double>> c{{0.9, -0.9, 0, 0, 0}
,{-0.3, 1, -0.7, 0, 0}
,{0, -0.4, 1, -0.6, 0}
,{0, 0, -0.5, 1, -0.5}
,{0, 0, 0, 0, 1}
};
// 常数列
vector<double> b{1, 1, 1, 1, 0};
int n = 5;
const double EPS = 1e-10;
// 高斯消元, 保证有唯一解
for(int i = 0; i < n; ++i)
{
// 找到 x[i] 的系数不为零的一个方程
for(int j = i; j < n ;++j)
{
if(fabs(c[j][i]) > EPS)
{
for(int k = 0; k < n; ++k)
swap(c[i][k], c[j][k]);
swap(b[i], b[j]);
break;
}
}
// 消去其它方程的 x[i] 的系数
for(int j = 0; j < n; ++j)
{
if(i == j) continue;
double rate = c[j][i] / c[i][i];
for(int k = i; k < n; ++k)
c[j][k] -= c[i][k] * rate;
b[j] -= b[i] * rate;
}
}
cout << std::fixed << std::setprecision(6);
for(int i = 0; i < n; ++i)
{
cout << "dp[" << i << "] = " << b[i] / c[i][i] << endl;
}
}
求解结果
dp[0] = 10.888889
dp[1] = 9.777778
dp[2] = 7.873016
dp[3] = 4.936508
dp[4] = 0.000000
蒙特卡洛模拟验证结果
import numpy as np
import bisect
class Game:
def __init__(self, p):
self.n = p.shape[0]
self.p = p
for i in range(self.n):
for j in range(1, self.p.shape[1]):
self.p[i][j] += self.p[i][j - 1]
def step(self, pos):
return bisect.bisect_left(self.p[pos]
,np.random.rand())
def __call__(self):
pos = 0
n_steps = 0
while pos != self.n:
pos = self.step(pos)
n_steps += 1
return n_steps
p = np.array([[0.1, 0.9, 0.0, 0.0, 0.0]
,[0.3, 0.0, 0.7, 0.0, 0.0]
,[0.0, 0.4, 0.0, 0.6, 0.0]
,[0.0, 0.0, 0.5, 0.0, 0.5]
])
game = Game(p)
def test():
T = int(1e7)
total_n_steps = 0
for t in range(T):
total_n_steps += game()
print("Average Steps: {:.6f}".format(total_n_steps / T))
for i in range(10):
test()
模拟结果
Average Steps: 10.890664
Average Steps: 10.887758
Average Steps: 10.882893
Average Steps: 10.889820
Average Steps: 10.891342
Average Steps: 10.888397
Average Steps: 10.888496
Average Steps: 10.891009
Average Steps: 10.890464
Average Steps: 10.886080
在前一小节中我们通过【有向图+期望DP+高斯消元】这个方法论解决了装备问题,题目描述和解法可以参考前一小节。
下面我们在装备升级问题上做一些变化,我们保持图结构不变,但是我们把期望 DP 转移方程中的边权改为节点的权。
我们考虑一个在游戏相关岗位面试常见的问题: 一次游戏通关带来的收入
一个游戏有四关,通过概率依次为0.9, 0.7, 0.6, 0.5。
第一关不收费,第二到四关每次收费分别为1块, 2块, 3块。
用户每玩一次都会无限打下去直至通关,通关后用户可以提现 10 块钱作为奖励。
问: 公司可以在每次用户游戏中平均挣多少钱。
思路参考
我们首先考虑【公司可以在每次用户游戏中收费多少钱】,然后减去 10 块钱的奖励就是挣的钱。
建图
图与上面的装备升级那道题一样,只是边权没有了,节点有权重表示费用。
期望 DP
期望 DP 的方程组如下:
dp[0] = dp[1]
dp[1] = 1 + 0.7 * dp[2] + 0.3 * dp[0]
dp[2] = 2 + 0.6 * dp[3] + 0.4 * dp[1]
dp[3] = 3 + 0.5 * dp[4] + 0.5 * dp[2]
dp[4] = 0
这是一个有环的图,且方程是线性方程组,面试的时候手推就可以。
如果要编程的话,需要整理成标准形式后用高斯消元解决。
高斯消元
标准形式如下
高斯消元求解
#include <vector>
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;
int main()
{
// 系数矩阵
vector<vector<double>> c{{1, -1, 0, 0, 0}
,{-0.3, 1, -0.7, 0, 0}
,{0, -0.4, 1, -0.6, 0}
,{0, 0, -0.5, 1, -0.5}
,{0, 0, 0, 0, 1}
};
// 常数列
vector<double> b{0, 1, 2, 3, 0};
int n = 5;
const double EPS = 1e-10;
// 高斯消元, 保证有唯一解
for(int i = 0; i < n; ++i)
{
// 找到 x[i] 的系数不为零的一个方程
for(int j = i; j < n ;++j)
{
if(fabs(c[j][i]) > EPS)
{
for(int k = 0; k < n; ++k)
swap(c[i][k], c[j][k]);
swap(b[i], b[j]);
break;
}
}
// 消去其它方程的 x[i] 的系数
for(int j = 0; j < n; ++j)
{
if(i == j) continue;
double rate = c[j][i] / c[i][i];
for(int k = i; k < n; ++k)
c[j][k] -= c[i][k] * rate;
b[j] -= b[i] * rate;
}
}
cout << std::fixed << std::setprecision(6);
for(int i = 0; i < n; ++i)
{
cout << "dp[" << i << "] = " << b[i] / c[i][i] << endl;
}
}
求解结果
dp[0] = 16.000000
dp[1] = 16.000000
dp[2] = 14.571429
dp[3] = 10.285714
dp[4] = 0.000000
因此公司可以在每次用户游戏中收费 dp[0] = 16 块钱,减去用户通关的 10 块钱奖金,公司可以挣 6 块钱。
蒙特卡洛模拟验证结果
import numpy as np
import bisect
class Game:
def __init__(self, transfer, cost):
self.n = transfer.shape[0]
self.cost = cost
self.transfer = transfer
for i in range(self.n):
for j in range(1, self.transfer.shape[1]):
self.transfer[i][j] += self.transfer[i][j - 1]
def step(self, pos):
return bisect.bisect_left(self.transfer[pos]
,np.random.rand())
def __call__(self):
pos = 0
pay = 0
while pos != self.n:
pay += self.cost[pos]
pos = self.step(pos)
pay += self.cost[self.n]
return pay
transfer = np.array([[0.1, 0.9, 0.0, 0.0, 0.0]
,[0.3, 0.0, 0.7, 0.0, 0.0]
,[0.0, 0.4, 0.0, 0.6, 0.0]
,[0.0, 0.0, 0.5, 0.0, 0.5]
])
cost = np.array([0.0, 1.0, 2.0, 3.0, -10.0])
game = Game(transfer, cost)
def test():
T = int(1e5)
total_pay = 0
for t in range(T):
total_pay += game()
print("Average income: {:.4f}".format(total_pay / T))
for i in range(10):
test()
模拟结果
Average income: 5.9768
Average income: 6.0043
Average income: 5.9052
Average income: 6.0664
Average income: 6.0720
Average income: 5.9687
Average income: 6.0101
Average income: 6.0450
Average income: 5.9871
Average income: 6.0704
“祝你好运”一种赌博游戏,经常在嘉年华和赌场玩。
玩家可以在 1, 2, 3, 4, 5, 6 中的某个数上下注。然后投掷 3 个骰子。
如果玩家下注的数字在这三次投掷中出现了 x 次,则玩家获得 x 倍原始投注资金,并且原始投注资金也会返还,也就是总共拿回 (x + 1) 倍原始投注资金。如果下注的数字没有出现,则原始投注资金不会返还。
问:每单位原始投注资金,玩家期望会输多少?
思路参考
设原始投注资金为 s,下注的数字在三次投掷中出现了 x 次,x 的取值为 0, 1, 2, 3,下面我们分别考虑这四种情况。
按照期望的定义
拿回的钱的期望为
每单位投注资金的损失为 17/216 = 0.0787
蒙特卡洛模拟
import numpy as np
def game():
x = 0
for _ in range(3):
if np.random.randint(1, 7) == 1:
x += 1
if x == 0:
return 0
return x + 1
def test():
T = int(1e6)
total_income = 0
for t in range(T):
total_income += game()
print("expected loss: {:.6f}".format((T - total_income) / T))
for i in range(10):
test()
模拟结果
expected loss: 0.078817
expected loss: 0.078182
expected loss: 0.078731
expected loss: 0.078875
expected loss: 0.078526
expected loss: 0.078657
expected loss: 0.078634
expected loss: 0.079033
expected loss: 0.079353
expected loss: 0.079307
阅读量:2029
点赞量:0
收藏量:0