【Puzzle】用户一次游戏带来的收入

  |  

摘要: 来自朋友的面试的概率题

【对算法,数学,计算机感兴趣的同学,欢迎关注我哈,阅读更多原创文章】
我的网站:潮汐朝夕的生活实验室
我的公众号:算法题刷刷
我的知乎:潮汐朝夕
我的github:FennelDumplings
我的leetcode:FennelDumplings


写在前面

概率面试题连载已经写了几期了,目前的题目主要是《Fifty challenging problems in probability with solutions》这本书中的题,以及朋友在面试中遇到的题。对于每道题,首先给出思路参考,然后用蒙特卡洛模拟的方式验证结果。

今天依然是看一道朋友面试时遇到的题,我认为本题是概率计算问题中最有代表性的一道题,可以延伸出很多变种,所以准备多交代一些相关内容,以后遇到这类题,90% 都可以不用怕了。

实际上在前面的概率面试题连载文章中,我们解决过类似的问题,只是那道题很简单,因此直接提了思路就直接过了,没有做方法论的总结。

本文会先以此前解决过的这个简单问题为例子,交代一些相关的背景,提炼出一类题的方法论,内容涉及到全期望公式,有向图,期望DP,高斯消元。

然后再引入今天这道题并解决,然后再看一个在今天的题目基础上的变种问题。


$1 回顾: 试验直到第一次成功

回顾在之前的文章中解决的一个问题,问题和解法如下,原文可以看 【Puzzle】试验直到第一次成功

一枚骰子掷出第一个 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|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] 的表达式了

1
E[0] = p * (1 + E[1]) + q * (1 + E[0])

这个表达式就是全期望公式,只是将全期望公式中 E(X|y) 的部分用图中点权和边权的方式定义了计算方法。

期望DP

通过以上的分析可以发现,当根据具体问题定义好有向图 D,起点 s,终点 t 后,我们就可以用类似于动态规划的方式求从 s 到 t 的某个指标的期望。

1
2
3
4
5
6
7
8
9
10
11
状态定义
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 状态时会出现类似于下面的情况(就是后面的【用户一次游戏带来的收入】那道题的转移方程,这里提前看一下)

1
2
3
4
5
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],没法求了。

如果方程组是线性方程组,还有有办法的,解决办法是利用高斯消元的方式整体解这个线性方程组。

关于高斯消元的推导、代码、题目,可以参考这篇文章,高斯消元与线性方程组


$2 方法总结

通过前面的分析,我大致可以总结出求解这类求期望的问题的方法论。

首先理论基础就是全期望公式。

然后针对具体的问题,我们先分析出建图需要的信息

  1. 起点 s 和终点 t 状态是什么,以及还有哪些状态节点

  2. 求出各个状态节点走一步可以走到哪些状态节点,概率分别多少

  3. 分析计算所求的期望是否需要额外的点权和边权

然后建图,初始化 dp[t]

分析状态转移方程中的 g, f 并写出实际的转移方程

最后求解 dp[s]:

  • 如果图没有环,直接按拓扑排序的顺序进行动态规划的方式推导即可;
  • 如果图有环,就看方程组是什么样的,线性方程组可以用高斯消元法求解,其它形式的方程我暂时也不清楚解法。

$3 题目: 装备升级

下面开始今天要看的朋友在面试时遇到的题目。题目描述如下:

玩家的装备有从 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 平均要经历多少个工序。

思路参考

建图

按照前面的方法总结。我们首先分析题目并建图,点状态是 0,终点状态是 4,此外还有 1, 2, 3 这三个中间状态。状态间的转移方向以及概率题目中明确给了。我们要求的是从 s=0 走到 t=4 平均需要走的步数。

由于每经过一条边都相当于走了一步,所以边上有额外的权 1。除此之外没有别的权,节点上也没有权。根据这些信息,我们首先把图建出来,如下

期望DP

期望 DP 的方程组如下:

1
2
3
4
5
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

这是一个有环的图,且方程是线性方程组,面试的时候手推就可以。

如果要编程的话,需要整理成标准形式后用高斯消元解决。

高斯消元

标准形式如下

高斯消元求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#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;
}
}

求解结果

1
2
3
4
5
dp[0] = 10.888889
dp[1] = 9.777778
dp[2] = 7.873016
dp[3] = 4.936508
dp[4] = 0.000000

蒙特卡洛模拟验证结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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()

模拟结果

1
2
3
4
5
6
7
8
9
10
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

$4 变种问题: 一次游戏通关带来的收入

下面我们在装备升级问题上做一些变化,我们保持图结构不变,但是我们把期望 DP 转移方程中的边权改为节点的权。

我们考虑一个在游戏相关岗位面试常见的问题: 一次游戏通关带来的收入

一个游戏有四关,通过概率依次为0.9, 0.7, 0.6, 0.5。
第一关不收费,第二到四关每次收费分别为1块, 2块, 3块。
用户每玩一次都会无限打下去直至通关,通关后用户可以提现 10 块钱作为奖励。
问: 公司可以在每次用户游戏中平均挣多少钱。

思路参考

我们首先考虑【公司可以在每次用户游戏中收费多少钱】,然后减去 10 块钱的奖励就是挣的钱。

建图

图与上面的装备升级那道题一样,只是边权没有了,节点有权重表示费用。

期望 DP

期望 DP 的方程组如下:

1
2
3
4
5
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

这是一个有环的图,且方程是线性方程组,面试的时候手推就可以。

如果要编程的话,需要整理成标准形式后用高斯消元解决。

高斯消元

标准形式如下

高斯消元求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#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;
}
}

求解结果

1
2
3
4
5
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 块钱。

蒙特卡洛模拟验证结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
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()

模拟结果

1
2
3
4
5
6
7
8
9
10
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

Share