更新说明：

2021/1/16日更新：修复了一个bug，之后执行了更大规模的仿真（五千亿次模拟抽卡）。文中所有的数据已经更新。

项目源码： GitHub repo - Arknights-Pull-Simulation

如果在 RPubs 页面内无法打开此链接，可以在新的浏览器标签页中打开。

The English version of this article: Estimating the Specific 6-Star Operator Dropping Rate in Arknights - the Monte Carlo Approach

此repo的仿真程序将通过蒙特卡洛模拟的方法，来估计明日方舟内不同类型卡池中在 “ 第 $i$ 次抽卡 ” 与 “ $i$ 次抽卡内 ” 获得特定 6★ 干员的概率值。

在本仿真程序中，保底系统的起点是可以调整的。你可以通过这个功能来探究在不同的数学模型中这些概率值会是多少。可以预计的是，不同的保底系统起点将会对你的抽卡概率有着不同的影响。

在本文中，我将首先介绍如何使用本项目的源码，并会给出一些简单易懂的例子以供尝鲜。随后，我会简单探讨下明日方舟抽卡规则背后的数学模型，以及你很可能会关心的一些随机事件的概率大小。接下来会介绍本程序的核心——随机数发生器。因为我觉得这部分内容应该会对你有所启发（例如 rand() % n 为何不好），所以也放到了这篇文章中。最后一部分是仿真数据与结果分析以及一些将来的工作。

你也可以直接跳到你感兴趣的部分（例如数学模型讨论或结果分析）：

什么是明日方舟
使用本项目
明日方舟中的数学模型
随机数的生成
- 不要再用 rand() % n 了！
- 推荐的方法 - C++11 <random>库
仿真结果分析
将来的工作
结语与声明

什么是明日方舟?

明日方舟是由鹰角网络旗下的蒙塔山工作室开发的一款塔防游戏，登录的平台包括 iOS、安卓与 PC 端。明日方舟于 2019 年 5 月 1 日与中国开始全平台公测，随后与 2020 年 1 月 16 日由悠星Yostar代理在世界范围内上线。（来自 维基百科 ）

明日方舟连续三次登顶 iOS 畅销榜榜首 [1]，并获得了 TapTap2019 - 年度最佳游戏、TapTap2019 - 年度最具影响力国产游戏、TapTap2019 - 年度最受玩家喜爱游戏 [2]，Bilibili 游戏区 2019 年度 - 最受好评游戏、Bilibili 游戏区 2019 年度 - 最多搜索游戏、Bilibili 游戏区 2019 年度 - 用户最喜爱的手游 [3] 等多个奖项。

全球范围内游戏下载量已突破 4 千万 [4]，其中中国大陆地区外下载量已于 2020 年 4 月末突破 7 百万 [5]。

官方网站：

使用本项目

构建代码

git clone 后, cd 进入代码目录并运行 make 。此时会生成一个名为 simulation_sequential 的可执行文件。

运行 make clean 会清除所有的 *.o 文件以及生成的可执行文件。

命令行参数

命令行参数如下：

./simulation_sequential [--help] [-t|--total-pull-time <value>] [--standard|--limited] 
                        [-p|--pity <value>] [-n|--num-rate-up <value>] [-c|--current-pull]

方括号中的参数是可选的，并且它们的顺序并无特别要求。另外，你既可以用短名称参数（-t, -p, -n 与 -c）也可以用对应的长名称参数（--total-pull-time, --pity, --num-rate-up 与 --current-pull），它们是等效的。

每一个参数的含义如下表所示：

参数	解释
`--help`	显示程序的帮助信息，介绍参数含义以及如何使用本程序不应该与其它参数同时提供
`-t` `--total-pull-time`	设置仿真抽卡的总次数应为一个介于 $[1, 18446744073709551615]$ (包含) 之间的整数
`--standard`	估计常规卡池中的概率理论值
`--limited`	估计限定卡池中的概率理论值
`-p` `--pity`	设置保底系统起作用的起点保底系统将在本参数所设置的第 `N` 抽后开始提升 6★ 干员的出率（ `N` 是你提供给本参数的值）应为一个介于 $[0, 4294967295]$ (包含) 之间的整数
`-n` `--num-rate-up`	设置卡池 UP 的干员数量应为 $1$ 或 $2$
`-c` `--current-pull`	设置你目前有多少次抽卡没抽到 6★ 干员应为一个介于 $[0, $ `<-c\|--current-pull value>` $ + 49 )$ (左闭右开) 的整数

默认状态下，如果无任何命令行参数，程序将会在保底系统起点为 $50$ 抽的条件下（与游戏内一样）仿真 $100,000,000$ 次抽卡，并估计双 UP 限定寻访卡池的理论概率值（相当于以参数--limited -t 100000000 -n 2 -p 50 运行程序），除非你通过提供相应的命令行参数来覆盖程序的默认设置。

程序额外说明

仿真抽卡的次数越多，所获得的结果就越精确，但这会消耗更多的时间去完成仿真。在我所使用的平台上，单线程仿真 $100$ 亿次仿真（-t 10000000000）所花费的时间约为 $110$ 秒。我使用的平台的技术参数为：
- 操作系统：Ubuntu 18.04
- CPU：Intel® Xeon® Gold 6140 CPU @ 2.30GHz
- 内存：8GB
你可以通过提供相应的命令行参数来覆盖程序的默认仿真设置。例如，如果你规定了 -t 200000000 （或 --total-pull-time 200000000）参数运行仿真程序，那么默认的仿真出卡次数（ $1$ 亿次）就会被覆盖为 $200000000$ 次，不过仍然会仿真并估计双 UP 限定寻访卡池的理论概率；如果你除了提供 -t 200000000 外还提供了 --standard 与 -n 1 ，那么程序将会通过 $200000000$ 次仿真抽卡去估计一个单 UP 标准寻访卡池的概率理论值。
仿真程序会像一个简易的编译器那样去尽力检查命令行参数中的语法错误（如下图所示），并会显示错误的原因。然而不要完全依靠这个功能，仔细阅读说明文档总是第一位的 ;-)

例子

以下是一些简明易懂的例子，如果你想立刻体验本程序，可以选择几个并运行：

./simulation_sequential --help

仿真并估计在双 UP 限定寻访卡池获得目标 6★ 干员的理论概率

./simulation_sequential --limited -t <大的正整数>

或者

./simulation_sequential -t <大的正整数>
仿真并估计在双 UP 标准寻访卡池获得目标 6★ 干员的理论概率

./simulation_sequential --standard -t <大的正整数>
仿真并估计在单 UP 标准寻访卡池获得目标 6★ 干员的理论概率

./simulation_sequential --standard -n 1 -t <大的正整数>
仿真并估计在单 UP 标准寻访卡池获得目标 6★ 干员的理论概率，但是已经在任一标准寻访卡池中 $42$ 次抽卡没有获得 6★ 干员

./simulaition_sequential --standard -n 1 -c 42 -t <大的正整数>
仿真并估计在单 UP 限定寻访卡池获得目标 6★ 干员的理论概率（还没有在明日方舟中出现！）

./simulation_sequential --limited -n 1 -t <大的正整数>

或者

`./simulation_sequential -n 1 -t <大的正整数>

<大的正整数> 的推荐值大于等于 $100000000$。

因为程序的输出会比较长，故建议把程序输出重定向到一个文本文件：
```
  ./simulation_sequential --limited -t 200000000 >> simu_limited_banner.res
```

明日方舟中的数学模型

在这部分，我将会简要探讨一下明日方舟抽卡规则背后的数学模型，本文中接下来要着重讨论的各种随机事件的定义，以及计算这些随机事件概率理论值的一些方法。

明日方舟抽卡规则与卡池类型

在明日方舟中，获得 6★ 干员（最稀有干员）的基础概率是 $2\%$ 。但游戏中有一种叫做保底系统的机制存在——它会记录你没有获得 6★ 干员的抽卡次数，如果你在连续 $50$ 抽中都没有获得 6★ 干员，那么在你接下来的每一抽中获得 6★ 干员的概率将提升 $2\%$ ，直到你在某一次抽卡中再次获得 6★ 干员。此时这个概率将会被重置为 $2\%$ ，并清空保底系统的抽卡计数。

举例来说，如果某玩家在他的连续 $50$ 抽中均未获得 6★ 干员，那么他在第 $51$ 抽中获得 6★ 干员的概率将会是 $4\%$ ；如果他在第 $51$ 抽里仍然没有获得 6★ 干员，那么他在第 $52$ 抽中获得 6★ 干员的概率将会是 $6\%$ ，以此类推。

这个机制可以确保即使是最不走运的玩家仍然可以在他的第 $99$ 抽里获得一个 6★ 干员。~~如果真是这样估计他就卸载明日方舟了…~~

连续 $98$ 抽都未获得 6★ 干员的概率，理论上是 \[ \begin{aligned} Pr(明日方舟，\ 卸载！) =& (1-2\%)^{50} \times (\prod_{i=51}^{98}(1-(2\% \times (i-50)+2\%))) \times 1\\ =& (0.98)^{50} \times \prod_{i=1}^{48} \frac{49-i}{50}\\ \approx & 1.27248 \times 10^{-21} \end{aligned} \]

这个概率远远远远小于一个（美国）人在一年中被闪电击中的概率 $2 \times 10^{-6}$ （数据来源于美国疾病控制预防中心 [6] ）。

在明日方舟中，干员可以从不同类型的卡池中获得。可以将这些卡池大体上分为两大类：

标准寻访卡池（standard headhunting banner，以下简称 “ 标准寻访 ” ）：不含有限定 6★ 角色（例如年，W 等）的卡池。这种卡池可以是双 UP 卡池（有两位 6★ 干员出率 UP），也可以是单 UP 卡池（有一位 6★ 干员出率 UP）。前者多为常驻卡池，并且会定期轮换，而后者常常随着新活动出现，也往往会有新的 6★ 干员登场。

例如，君影轻灵，瑕光微明，以及常驻标准寻访 #39 都是标准寻访。
限定寻访卡池（limited headhunting-banner，以下简称 “ 限定寻访 ” ）：含有限定 6★ 干员的卡池。目前这种卡池都是双 UP 卡池。

例如，地生五金，遗愿焰火，以及勿忘我 都是限定寻访。

另外，明日方舟中还存一种卡池——联合行动。这种卡池也是一种标准寻访，但是有四位出率 UP 的 6★ 干员。目前本程序还不支持这类卡池的模拟，但我会很快增加这个功能。

在任何一个标准寻访中，没有获得 6★ 干员时，都会累积次数，该次数不会因为标准寻访的结束而清零。因为累积次数而增加的获得概率，也会应用于接下来任意一次标准寻访。（摘自游戏内描述 ）

抽卡次数可以在卡池之间继承的规则只适用于标准寻访，而不适用于限定寻访。限定卡池结束后这个次数将会清零，且限定寻访与标准寻访中的累积次数互不影响。

在标准寻访中，如果你在一次抽卡中获得了 6★ 干员，那么这位 6★ 干员会有 $50\%$ 的（条件）概率是本卡池出率 UP 的干员（之一）；而在限定寻访中，如果你在一次抽卡中获得了 6★ 干员，那么这位 6★ 干员会有 $70\%$ 的（条件）概率是本卡池出率 UP 的干员。

假设你的目标干员（target 6★ operator）是在卡池中被出率 UP 的那一个，那么用条件概率可表达为 \[ \begin{aligned} &Pr(获得目标\ 6★\ 干员\ |\ 在标准寻访卡池中获得\ 6★\ 干员) = 50\% \\ &Pr(获得目标\ 6★\ 干员\ |\ 在限定寻访卡池中获得\ 6★\ 干员) = 70\% \end{aligned} \]

在这里，我们会做一个重要的假设：如果 $n$ 个 6★ 干员在某一卡池中获得了出率 UP ，他们会平分上面所说的条件概率。也就是说，在获得 6★ 干员的条件下，这位 6★ 干员是卡池 UP 干员其中之一的条件概率为 $50\% / n$ （标准寻访中）或 $70\% / n$ （限定寻访中）。个人认为折买目前看起来是十分合理的，即使明日方舟官方并未主动声明这一点。

另外，明日方舟还存在着专门为限定寻访设置的第二个保底机制——玩家一般把它称作 “ 井 ” 或者 “ 硬保底 ” 机制（即外服玩家口中的 Spark System）。每在限定寻访中抽卡一次，你就会获得一个寻访数据契约 代币，本限定池中出率 UP 的两位 6★ 干员都可以通过花费 $300$ 个寻访数据契约 直接购买得到。~~你不会想用这种方法去拿限定 6★ 的…~~

随机事件定义以及所研究的问题

接下来，我们对以下随机事件做出定义，

$A_i$ : 在第 $i$ 抽才获得任意 6★ 干员
$B_i$：$i$ 抽之内都没获得任意 6★ 干员
$S_i$ : 在第 $i$ 抽才获得目标 6★ 干员
$W_i$ : 在 $i$ 抽之内获得目标 6★ 干员
$F_i$ : $i$ 抽之内都没获得目标 6★ 干员

下标 $i$ 用于计数抽卡次数，当前抽卡的结果将决定是上面的哪种随机事件发生。此处，目标 6★ 干员指的是卡池中出率 UP 的 6★ 干员之一。另外根据定义，$A_i$ 与 $S_i$ 都要求前 $i-1$ 次抽卡中未获得任意 6★ 干员 / 目标 6★ 干员。

为了使下标更简洁，

$A_i$ 的下标 $i$ 从 $1$ 开始计数，直到你获得了 6★ 干员使得 $i$ 重新从 $1$ 开始计数。根据保底系统的规则，$i$ 的最大值为 $99$ 。
当你某次抽卡获得 6★ 干员后，$B_i$ 的下标 $i$ 重置为 $0$ ，并且接下来如果每次抽卡未获得 6★ 干员时，$i$ 增加 $1$ 。根据保底系统的规则，$i$ 的最大值是 $98$ 。$i$ 可以理解为未获得 6★ 干员的抽卡次数。
$S_i$ 与 $W_i$ 的下标 $i$ 从 $1$ 开始计数，直到你获得了目标 6★ 干员使得 $i$ 重新从 $1$ 开始计数。因为总有机会抽到的 6★ 干员是非卡池出率 UP 干员，所以 $i$ 可以是$+ \infty$ 。
当你某次抽卡获得目标 6★ 干员后，$F_i$ 的下标重置为 $0$ ，并且如果每次抽卡后未获得目标 6★ 干员时，$i$ 增加 $1$ 。需要注意的是当你获得非目标 6★ 干员时 $i$ 并不会重置为 $0$ 。因为总有机会抽到的 6★ 干员是非卡池出率 UP 干员，所以 $i$ 可以是$+ \infty$ 。$i$ 可以理解为未获得目标 6★ 干员的抽卡次数。

总而言之，当你达成你的目标（出 6★ 或者出目标 6★ 后），下标 $i$ 就会重置计数。

基于上面的定义，这里给几个例子方便理解。准确地理解这些定义将有助于阅读文章后面的各种讨论以及结果分析。

$A_{42}$ 的意思是自从保底系统重置后，在接下来的 $41$ 次抽卡中并未获得 6★ 干员，并在第 $42$ 次抽卡中获得了某一 6★ 干员（只要是 6★ 干员就可）。也就是说 $A_{42}$ 的发生一定意味着 $B_{41}$ 的发生。
$B_{37}$ 的意思是自从保底系统重置后，在接下来的 $37$ 次抽卡中均未获得 6★ 干员。特别地，$B_0$ 与 $A_1$ 的含义相同（即在保底系统重置后的新一系列抽卡中，第一抽就抽中了一位 6★ 干员）。
$S_{54}$ 的意思是自从保底系统重置后，在接下来的 $53$ 次抽卡中并未获得目标 6★ 干员，并在第 $54$ 次抽卡中获得了目标 6★ 干员。也就是说 $S_{54}$ 的发生一定意味着 $F_{53}$ 的发生。
$W_{25}$ 的意思是保底系统重置后，在接下来的 $25$ 次抽卡中获得目标 6★ 干员。注意，在这 $25$ 次抽卡内，只要你获得了目标 6★ 干员，你就停止抽卡。（否则如果继续抽卡的话，仍有可能在后面几次抽卡中再次获得目标 6★ 干员，即 “ 在 $25$ 次抽卡内至少获得一次目标 6★ 干员 ” 。这个事件的概率计算会稍微复杂些。）
$F_{17}$ 的意思是自从保底系统重置后，在接下来的 $17$ 次抽卡中均未获得目标 6★ 干员。特别地，$F_0$ 与 $S_1$ 的含义相同（即在保底系统重置后的新一系列抽卡中，第一抽就抽中了目标 6★ 干员）。

目标 6★ 干员是指卡池出率 UP 的 6★ 干员之一，具体是哪一位是由你决定的。

如果你仔细斟酌这些定义，你会发现 $Pr(A_i)$ ，$Pr(B_i)$ ，$Pr(S_i)$ ，$Pr(W_i)$ 与 $Pr(F_i)$ 本身也是一个条件概率，即在事件 “ 保底系统重置 ” 发生的条件下，在接下来的第 $i$ 次抽卡 / $i$ 次抽卡内获得 / 未获得 6★ 干员 / 目标 6★ 干员的概率。这也是文章后面的仿真结果的那些数据的含义。如果你已经在上一个卡池中抽卡数次且未获得 6★ 干员，在新开的卡池里继续抽卡，在第 $i$ 次抽卡 / $i$ 次抽卡之内获得 / 不获得 6★ 干员 / 目标 6★ 干员概率可能会有些许不同，但我目前认为差别应该不会很大（~~希望不会被打脸~~）。

我们研究的问题是当 $i$ 取不同值时， $S_i$ 发生的概率 $Pr(S_i)$ 以及 $F_i$ 发生的概率 $Pr(F_i)$ 分别是多少。

6★ 干员出率：概率与期望值

计算获得任意 6★ 干员的概率是相对容易的。一些先前的工作 [7] 已经得到了从 $Pr(A_1)$ 到 $Pr(A_{99})$ 的数值。但是 [7] 中并没有给出理论公式，在此我把它们补充上。理论上， \[ \begin{equation} Pr(A_i)=\left\{ \begin{aligned} & (0.98)^{i-1} \times 0.02&1 \le i\le 50 & \\ & (0.98)^{50} \times (\prod_{k=51}^{i-1} \frac{99-k}{50}) \times \frac{i-49}{50} & 50<i\le 98 && \ 其中\ i \in N_+\\ & (0.98)^{50} \times \prod_{k=1}^{48} \frac{49-k}{50} & i = 99 & \\ \end{aligned} \right. \end{equation} \]

基于上面的式子，我们可以得到获得 6★ 干员的期望的抽卡次数，也就是 \[ \sum_{i=1}^{99}i \cdot Pr(A_i) \] 这个值大约是 $34.59$ ，或者说，平均下来在一次抽卡中获得 6★ 干员的概率是 $2.89\%$

这里，期望值 $34.59$ 的含义是，如果我们把所有玩家抽了多少次之后才获得一位 6★ 干员的抽卡次数汇总起来，去计算这些抽卡次数的平均数，这个数值将会是 $34.59$ 。也可以从另一个角度去理解——如果你抽卡 $N$ 次，在这 $N$ 次抽卡中获得了 $m$ 次 6★ 干员。每一次获得 6★ 干员你抽卡抽了 $n_i$ 次，这里 $i \in N_+,\ 1 \le i \le m$ ，则会有下式成立 \[ \begin{aligned} \lim_{m \to +\infty} \frac {\displaystyle \sum_{i=1}^{m} n_i }{m} =& \lim_{N \to +\infty} \frac {\displaystyle \sum_{i=1}^{m} n_i }{m}\\ =& \sum_{i=1}^{99}i \cdot Pr(A_i)\\ \approx& 34.59 \end{aligned} \]

事实上，根据 [8]，在当前次抽卡中获得 6★ 干员的概率仅取决于上一次抽卡有没有获得 6★ 干员，这恰好符合了马尔科夫链的性质——无记忆性。所以，我们可以用一个含有 $100$ 个状态的一阶马尔科夫链去为这个问题（即探究各个 $Pr(A_i)$ 的值）建立模型，其中各状态的含义即为 $B_i$ ，即当前有多少次抽卡未出 6★ 干员，如下图所示：

目标 6★ 干员出率：获得概率值的尝试与方法

当我们想要计算目标 6★ 干员的获得概率时，问题就变得复杂多了。假设你想计算出在你的第 $N^{th}$ 次抽卡获得你的目标 6★ 干员的概率，在你第 $N$ 次抽卡得到这位角色之前，在前面的 $N-1$ 次抽卡中你既有可能获得了数个非目标 6★ 干员、数次重置了保底系统的计数，也有可能一次 6★ 干员也没获得。可见，第 $N^{th}$ 次抽卡获得你的目标 6★ 干员可以对应到非常多种不同且独立的情况，而你需要分别计算出每一种分情况所发生的概率，再把所有这些概率相加。

让我们来具体数一数 “ 第 $N$ 次抽卡获得目标 6★ 干员 ” 到底能对应多少种不同且独立的情况。我们可以根据在前 $N-1$ 次抽卡中是否获得非目标 6★ 干员来去数：可能一位非目标 6★ 干员也没有获得，即 $\binom{N-1}{0}$ 种情形；也可能只获得了一位非目标 6★ 干员，此时可能是第一抽就获得非目标 6★ 干员，也可能是第二抽，也可能是第三抽……故一共有 $\binom{N-1}{1}$ 种可能的情形；也可能获得了两位非目标 6★ 干员，类似地，共有 $\binom{N-1}{2}$ 种可能的情形。所以， “ 第 $N$ 次抽卡获得目标 6★ 干员 ” 所对应的情形一共有 \[ \sum_{i=0}^{N-1}\binom{N-1}{i} = 2^{N-1} \]

我们需要把所有的这些情形发生的概率先计算出来，然后再求和，才能得到事件 “ 第 $N$ 次抽卡才获得目标 6★ 干员 ” 发生的概率。这意味着，如果我们采取这种直接的方法去计算 $Pr(S_i)$ 的理论值，这背后的计算量会相当庞大。

也许可行的获得理论值的方法

一些被认为较为有希望能计算出 $Pr(S_i)$ 精确理论值的方法包括马尔科夫链和 / 或动态规划 [9] 。在 [9] 中，$Pr(S_i),\ i \in N_+,\ 1\le i \le 300$ 的值已被求出，其作者考虑到硬保底机制，并没有计算 $Pr(S_i),\ i \in N_+,\ i \ge 300$ 的值。然而，有其他网友指出在作者的代码中使用了错误的下标 [10] ，故结果也是有误差的，另外，这种方法并不能获得 $Pr(S_i)$ 的代数表达式。

本项目所采取的方法：蒙特卡洛方法

在本项目中，我采用了蒙特卡罗法来估计 $Pr(S_i)$ 的理论值。

我们定义，一次试验是指 ” 进行多次抽卡，直到你获得你的目标 6★ 干员 “ 。模拟 $N$ 次抽卡（$N$ 足够大），设最后一共获得了 $n$ 次目标 6★ 干员。这意味着完成了 $n$ 个试验。设进行了 $m_i$ 次抽卡完成了第 $i$ 个试验，即试验 $1$ ，试验 $2$ ，试验 $3$ ，$...$ ，分别花费了 $m_1,\ m_2,\ ...,\ m_n$ 次抽卡才完成。我们认为其发生的频率近似等于其概率理论值， \[ \frac{m_i}{n} \approx Pr(S_i),\ i \in N_+,\ 1 \le i \le n \]

同时有下式成立：

\[ \begin{aligned} &\lim_{n \to + \infty}\frac{m_i}{n}\\ =& \lim_{N \to + \infty}\frac{m_i}{n}\\ =&Pr(S_i),\ i \in N_+ \end{aligned} \]

这提示我们，为了能获得一个更为精确的结果，我们在程序中所模拟的抽卡次数要尽可能的多，即，在平台性能允许的情况下，我们提供给命令行参数 -t 或 --total-pull-time 的值要尽可能的大。

在程序仿真过程中，如果当前已经执行的模拟抽卡次数已达到了命令行参数 -t 或 --total-pull-time 所规定的值，则会结束本次程序的执行。所以，有可能最后一次试验会被提前终止。

现在来看一下如何计算 $Pr(W_i)$ 。根据随机事件 $W_i$ 与 $S_i$ 的定义，我们可以得知 $W_i$ 意味着事件 $S_1,\ S_2,\ ...,\ S_{i-1},\ S_i$ 的其中之一发生。举例来说，$W_4$ 意味着你在第一次抽卡（根据前文所述，是保底系统重置后的第一次抽卡，以下省略）就获得了目标 6★ 干员（$S_1$），或者第二次抽卡获得了目标 6★ 干员（$S_2$），获得在第三次抽卡时获得（$S_3$），当然也有可能第四抽才获得（$S_4$）。正如前文所述，当你在第 $1$ 到 $4$ 抽中的某一抽获得时，就不在继续抽卡了，这样就可以避免计算一些较为复杂的情况（例如在这四次抽卡里 $S_1$ 发生之后又发生了 $S_3$ ，或者发生了两次 $S_2$ ，或者先发生了 $S_3$ 之后又发生了 $S_1$，又或者欧皇附体发生了四次 $S_1$）。

所以我们有下式：

\[ \begin{aligned} Pr(W_i)=& \sum_{k=1}^{i} Pr(S_k)\\ =& \sum_{k=1}^{i} \frac{m_k}{n},\ n\ 足够大 \end{aligned} \]

随机数的生成

对明日方舟抽卡的仿真涉及到随机数的生成。为了按照游戏内相同的抽卡规则模拟这一过程，并确保仿真结果的精确性，我们需要生成大量且高质量均匀分布的随机数，并把生成得到的数字与一些阈值相比较，来据此判断本次模拟的抽卡是否获得了 6★ 干员或目标 6★ 干员。所以，如何快速地生成高质量分布的随机数便成为了本仿真程序的核心。

不要再用 `rand() % n` 了！

当我们想要获得一个在某区间内呈均匀分布的随机数，例如在 $[0,\ 99]$ 内（闭区间），我们可能往往会这么写：

int seed;  // 种子被初始化为某一具体值
srand(seed);
for(int i = 0; i < 16; ++i) {
    int rand_num = rand() % 100;
    cout << rand_num << endl;
}

顺便一提，上面的代码，以及下面讨论的一些内容取自来自微软^®的大佬 Stephan T. Lavavej 的一场关于用 C++ 进行随机数生成的讲座。如果你对这部分内容感兴趣，个人十分推荐去听一下他的这场讲座。讲座连接 GoingNative2013 - rand() Considered Harmful

而这么写是有问题的！

首先，我们来看一下生成出来的随机数的周期。在 glibc 中，rand() 使用线性同余方法（Linear Congruential Generator，LCG）来生成（伪）随机数。它也有可能使用线性反馈移位寄存器（Linear-Feedback Shift Register，LFSR）去生成。根据 [11] ，如果用的是 LFSR ，rand() 只会记录最近的 $32$ 个输出的结果，并根据这 $32$ 个结果去运算下一个随机数的值。换句话说，它最多能拥有 $2^{32}$ 种不同的状态；如果使用的是 LCG ，则所生成的随机数的周期最大不会超过在 LCG 算法中用于取模的 m 的值 [12] 。事实上， [13] 与 [14] 表明，在 glibc 的实现里，m 的值是 $2^{31}-1$ （也就是十六进制下的 0x7fffffff ）。[15] 与 [16] 的一篇文章指出， rand() 的周期大约在 $2^{32}$ 上下。而在我们这里，为了较为精确地估计 $Pr(S_i)$ 的值，我们至少需要进行十亿次（大于 $2^{29}$ ）模拟的抽卡（每次模拟抽卡都需要生成一个随机数）。如果要提高精确值，我们则需要大约一百亿次（大于 $2^{33}$ ）模拟抽卡。可知，rand() 的周期对于我们的大规模的蒙特卡洛仿真来说是不够的。

[15] 的原文网址已经失效了，我这里提供的是谷歌的网页快照的连接。

其次，通过 ramd() 与 rand() % n 所获得的随机数的分布是低质量的，并不是严格的均匀分布。在 C++ 标准里就没有对 rand() 生成出的随机数的分布做过任何要求 [17] 。事实上，因为采用的生成算法是 LCG 或 LFSR ，生成出的数字就已经不是均匀分布了。即使退一步说，假设 rand() 生成的是严格的均匀分布，rand() % n 或所得的数字也仍然不是均匀分布。在 Stephan 的讲座的例子中，假设 rand() 给出的随机数是在 $[0,\ 32767]$ 上的完美的均匀分布，并且我们想获得在 $[0,\ 99]$ 上的均匀分布，如下面的代码所示

int src = rand();  // 假设在[0, 32767]上均匀分布
int desired_range = 100;
int dist = src % desired_range;  // 这么写是错误的！

这种错误是非常容易犯的。为什么用这种方法所获得的并不是均匀分布呢？原因其实很简单

int src = rand();  // 假设在[0, 32767]上均匀分布
int desired_range = 100;
int dist = src % desired_range;
//          src     →     dist
//        [0, 99]   →   [0, 99]
//     [100, 199]   →   [0, 99]
//     [200, 299]   →   [0, 99]
// ...
// [32700, 32767]   →   [0, 67] -- 问题出在这儿!

我们发现，获得 $[0,\ 67]$ 之间的数字的概率会比获得在 $(67,\ 99]$ 上的概率略微高一点。这个是取模运算导致的。只要 rand() 所能生成的数字的最大值（由头文件中的 RAND_MAX 确定）不能被 desired_range 整除，这个问题就会发生。当 $Pr(S_i)$ 本身就是一个很小的值的时候，这种由取模运算所产生的误差将会对 $Pr(S_i)$ 有着比较显著的影响，会导致相对误差很大。

另外，下面的写法也不是一个完美的均匀分布，因为只有当 src 取到 $99$ 时才会让 dist 取到 $99$ ，而对于小于 $99$ 的 dist 值可以对应于多个 src 的取值。

int src = rand();  // 假设在[0, 32767]上均匀分布
int desired_range = 100;
int dist = static_case<int>(src * 1.0 / RAND_MAX) * (desired_range - 1);

下面的这种写法也是错误的，因为根本存在在一种映射能让 $32768$ 个整数均匀地映射到 $100$ 个整数上去… 说到底问题的本质是抽屉原理。

int src = rand();  // 假设在[0, 32767]上均匀分布
int desired_range = 100;
int dist = static_case<int>(src * 1.0 / (RAND_MAX + 1)) * desired_range;

在 Stephan 的讲座中，他对这些例子做了更详细的解释，如果感兴趣的话可以去查看原讲座内容。

仿真结果分析

由于篇幅所限，在这部分我将只展示当 $i \in N_+, 1 \le i \le 100$ 时的 $Pr(S_i)$ 与 $Pr(W_i)$ 的值，以及根据更大范围内的值所绘制的折线图。你可以在 /res 目录底下找到当 $i \in N_+, 1 \le i < 1000$ 时的全部仿真结果。

如果你不清楚这里的 $Pr(S_i)$ 与 $Pr(W_i)$ 代表什么，可以参考本文的随机事件定义以及所研究的问题部分。

双 UP 限定寻访卡池的概率估计值

下面的结果是由两千亿次模拟抽卡仿真得到，即以下面的命令运行程序

./simulation_sequential --limited -t 500000000000 -n 2 -p 50 -c 0

也可以简单地以这条命令运行
./simulation_sequential -t 500000000000
效果是一样的。

再次推荐将程序输出重定向到一个文本文档，因为输出将会有一千多行。

下面的表格展示了前 $100$ 次抽卡的 $Pr(S_i)$ 与 $Pr(W_i)$ 的数值估计值。

$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$
1	0.700%	0.700%	26	0.587%	16.694%	51	0.748%	30.366%	76	0.457%	53.095%
2	0.695%	1.395%	27	0.583%	17.277%	52	0.980%	31.346%	77	0.453%	53.548%
3	0.690%	2.085%	28	0.579%	17.856%	53	1.180%	32.527%	78	0.448%	53.996%
4	0.685%	2.771%	29	0.575%	18.431%	54	1.338%	33.864%	79	0.444%	54.440%
5	0.681%	3.451%	30	0.571%	19.002%	55	1.445%	35.309%	80	0.440%	54.880%
6	0.676%	4.128%	31	0.567%	19.569%	56	1.501%	36.809%	81	0.436%	55.316%
7	0.671%	4.799%	32	0.563%	20.132%	57	1.506%	38.315%	82	0.432%	55.748%
8	0.666%	5.465%	33	0.559%	20.691%	58	1.467%	39.783%	83	0.428%	56.176%
9	0.662%	6.127%	34	0.555%	21.246%	59	1.394%	41.177%	84	0.425%	56.601%
10	0.657%	6.784%	35	0.551%	21.797%	60	1.296%	42.473%	85	0.421%	57.022%
11	0.653%	7.437%	36	0.548%	22.344%	61	1.184%	43.657%	86	0.417%	57.439%
12	0.648%	8.085%	37	0.543%	22.888%	62	1.068%	44.725%	87	0.413%	57.852%
13	0.643%	8.728%	38	0.540%	23.428%	63	0.955%	45.680%	88	0.410%	58.262%
14	0.639%	9.367%	39	0.536%	23.964%	64	0.851%	46.531%	89	0.406%	58.668%
15	0.634%	10.002%	40	0.532%	24.496%	65	0.761%	47.292%	90	0.402%	59.070%
16	0.630%	10.632%	41	0.529%	25.025%	66	0.686%	47.978%	91	0.399%	59.468%
17	0.626%	11.257%	42	0.525%	25.549%	67	0.626%	48.604%	92	0.395%	59.863%
18	0.621%	11.879%	43	0.521%	26.071%	68	0.579%	49.183%	93	0.391%	60.255%
19	0.617%	12.495%	44	0.518%	26.588%	69	0.544%	49.727%	94	0.388%	60.643%
20	0.613%	13.108%	45	0.514%	27.102%	70	0.518%	50.246%	95	0.384%	61.027%
21	0.608%	13.716%	46	0.510%	27.612%	71	0.500%	50.745%	96	0.381%	61.408%
22	0.604%	14.320%	47	0.507%	28.119%	72	0.486%	51.231%	97	0.378%	61.786%
23	0.600%	14.920%	48	0.503%	28.622%	73	0.476%	51.707%	98	0.374%	62.160%
24	0.595%	15.515%	49	0.500%	29.122%	74	0.468%	52.176%	99	0.371%	62.531%
25	0.591%	16.107%	50	0.496%	29.618%	75	0.462%	52.638%	100	0.368%	62.898%

下面的折线图展示了前 $500$ 次抽卡的 $Pr(S_i)$ 与 $Pr(W_i)$ 的值。

蓝色线条代表 $Pr(S_i)$ 的估计值（左侧纵坐标轴），橙色线条代表 $Pr(W_i)$ 的估计值（右侧纵坐标轴）。

下面是一些较为明显的特征：

在第 $57$ 次抽卡中从双 UP 限定寻访卡池中抽到目标 6★ 干员的概率最大；$Pr(S_{57})=1.506\%$ ；
如果你的游戏资源可以让你抽卡 $41$ 次，那么你在这 $41$ 次抽卡之内抽到目标 6★ 干员的概率是 $25\%$ 。
如果你的游戏资源可以让你抽卡 $70$ 次，那么你在这 $7$0 次抽卡之内抽到目标 6★ 干员的概率是 $50\%$ 。
如果你的游戏资源可以让你抽卡 $132$ 次，那么你在这 $132$ 次抽卡之内抽到目标 6★ 干员的概率是 $75\%$ 。
如果你的游戏资源可以让你抽卡 $212$ 次，那么你在这 $212$ 次抽卡之内抽到目标 6★ 干员的概率是 $90\%$ 。

可以看到，在 $Pr(S_i)$ 的曲线上存在至少 $3$ 个明显的小峰，这是由于保底系统引起的。有趣的是，即使获得 6★ 干员的概率从 $50$ 次抽卡之后就开始增加了，$Pr(S_i)$ 的最大值却在 $i=57$ 时才出现，即在保底系统起作用的 $7$ 抽之后。这里 $Pr(S_{57})=1.506\%$ 。$Pr(S_i)$ 的第二极大值出现在 $i=114$ 处，$Pr(S_{114})=0.466\%$ 。第三极大值出现在 $i=168$ 处，$Pr(S_{168})=0.204\%$ 。

现在我们看一下 $Pr(S_i)$ 的值随着 $i$ 的增大由降转升的点，方便起见，在本文中称它们为转折点（turning point，不是拐点 inflection point）。根据仿真得到的数据，我们发现在 $i = 51,\ 103,\ 164$ 这三个地方 $Pr(S_i)$ 的值都比它的前一个值大，并会继续上升一段距离。

现在来看一下在双 UP 限定池中玩家吃到硬保底的概率（即连续 $300$ 抽都没有抽到目标 6★ 干员的概率）。其实，$Pr(F_{300})$ 的值是很小的，我们可以看一下 $Pr(S_{301})$ 的值。根据定义，$S_{301}$ 意味着 $F_{300}$ 一定发生。这里，$Pr(S_{301})=0.041 \%$ ，它等于 $Pr(F_{300}) \times Pr(在第301抽获得目标\ 6★\ 干员)$ 。注意这里的随机事件 ” 在第 $301$ 次抽卡中获得目标 6★ 干员 “ 并不关心前 $300$ 抽卡的结果。我们有 \[ \begin{aligned} \max Pr(F_{300}) & = \max \frac{Pr(S_{301})}{Pr(在第301抽获得目标\ 6★\ 干员)}\\ & = \frac{Pr(S_{301})}{\min Pr(在第301抽获得目标\ 6★\ 干员)}\\ & = \frac{0.041\%}{2\% \times 70\% \div 2}\\ & \approx 5.587\% \end{aligned} \]

也就是说玩家吃硬保底的概率不会超过$5.587\%$ 。

应该存在一种基于 $Pr(S_i)$ 的值去计算 $Pr(F_i)$ 的方法，可能是 $1-Pr(W_i)$ 。目前我还不太确定我是否是对的😭…

双 UP 标准寻访卡池的概率估计值

下面的结果是由两千亿次模拟抽卡仿真得到，即以下面的命令运行程序

./simulation_sequential --standard -t 500000000000 -n 2 -p 50 -c 0

也可以简单地以这条命令运行
./simulation_sequential --standard -t 500000000000
效果是一样的。

下面的表格展示了前 $100$ 次抽卡的 $Pr(S_i)$ 与 $Pr(W_i)$ 的数值估计值。

$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$
1	0.500%	0.500%	26	0.441%	12.218%	51	0.571%	22.739%	76	0.406%	40.613%
2	0.498%	0.998%	27	0.439%	12.657%	52	0.738%	23.477%	77	0.403%	41.016%
3	0.495%	1.492%	28	0.437%	13.094%	53	0.884%	24.361%	78	0.400%	41.416%
4	0.493%	1.985%	29	0.435%	13.528%	54	0.998%	25.359%	79	0.397%	41.813%
5	0.490%	2.475%	30	0.432%	13.961%	55	1.078%	26.437%	80	0.395%	42.208%
6	0.488%	2.963%	31	0.430%	14.391%	56	1.121%	27.557%	81	0.393%	42.601%
7	0.485%	3.448%	32	0.428%	14.819%	57	1.128%	28.685%	82	0.390%	42.991%
8	0.483%	3.931%	33	0.426%	15.245%	58	1.104%	29.790%	83	0.388%	43.379%
9	0.480%	4.411%	34	0.424%	15.669%	59	1.055%	30.845%	84	0.385%	43.764%
10	0.478%	4.889%	35	0.422%	16.090%	60	0.988%	31.833%	85	0.383%	44.146%
11	0.476%	5.365%	36	0.419%	16.510%	61	0.911%	32.744%	86	0.380%	44.527%
12	0.473%	5.838%	37	0.418%	16.927%	62	0.830%	33.574%	87	0.378%	44.904%
13	0.471%	6.308%	38	0.415%	17.343%	63	0.751%	34.325%	88	0.375%	45.280%
14	0.468%	6.777%	39	0.413%	17.756%	64	0.679%	35.004%	89	0.373%	45.653%
15	0.466%	7.243%	40	0.411%	18.167%	65	0.616%	35.621%	90	0.371%	46.024%
16	0.464%	7.707%	41	0.409%	18.576%	66	0.564%	36.185%	91	0.368%	46.392%
17	0.462%	8.168%	42	0.407%	18.983%	67	0.522%	36.706%	92	0.366%	46.758%
18	0.459%	8.627%	43	0.405%	19.388%	68	0.489%	37.195%	93	0.364%	47.122%
19	0.457%	9.084%	44	0.403%	19.791%	69	0.465%	37.660%	94	0.361%	47.483%
20	0.454%	9.539%	45	0.401%	20.192%	70	0.447%	38.107%	95	0.359%	47.842%
21	0.452%	9.991%	46	0.399%	20.591%	71	0.434%	38.541%	96	0.357%	48.199%
22	0.450%	10.441%	47	0.397%	20.988%	72	0.425%	38.966%	97	0.355%	48.553%
23	0.448%	10.889%	48	0.395%	21.383%	73	0.418%	39.384%	98	0.352%	48.906%
24	0.445%	11.334%	49	0.393%	21.777%	74	0.413%	39.798%	99	0.350%	49.256%
25	0.443%	11.777%	50	0.391%	22.168%	75	0.409%	40.207%	100	0.348%	49.603%

下面的折线图展示了前 $500$ 次抽卡的 $Pr(S_i)$ 与 $Pr(W_i)$ 的值。

line-chart-double-rate-up-standard-banner-zh-CN

蓝色线条代表 $Pr(S_i)$ 的估计值（左侧纵坐标轴），橙色线条代表 $Pr(W_i)$ 的估计值（右侧纵坐标轴）。

下面是一些较为明显的特征：

在第 $57$ 次抽卡中从双 UP 标准寻访卡池中抽到目标 6★ 干员的概率最大；$Pr(S_{57})=1.128\%$ ；
如果你的游戏资源可以让你抽卡 $54$ 次，那么你在这 $54$ 次抽卡之内抽到目标 6★ 干员的概率是 $25\%$ 。
如果你的游戏资源可以让你抽卡 $102$ 次，那么你在这 $102$ 次抽卡之内抽到目标 6★ 干员的概率是 $50\%$ 。
如果你的游戏资源可以让你抽卡 $187$ 次，那么你在这 $187$ 次抽卡之内抽到目标 6★ 干员的概率是 $75\%$ 。
如果你的游戏资源可以让你抽卡 $303$ 次，那么你在这 $303$ 次抽卡之内抽到目标 6★ 干员的概率是 $90\%$ 。

可以看出，在双 UP 标准寻访卡池中，获得目标 6★ 干员比在限定寻访卡池中更为困难，因为目标 6★ 干员所占有的条件概率更低（每个目标 6★ 干员只拥有 $25\%$ 的条件概率，而在限定寻访卡池中这一条件概率是 $35\%$）。

同样由于保底系统，在 $Pr(S_i)$ 的曲线上也存在三个可以被观察到的转折点以及小高峰：

转折点： $i=51,\ 103,\ 163$
局部极大值及其位置：
- $Pr(S_{57})=1.128\%,\ 当\ i=57\ 时$
- $Pr(S_{114})=0.439\%,\ 当\ i=114\ 时$
- $Pr(S_{169})=0.241\%,\ 当\ i=169\ 时$

有趣的是，即使转折点的位置与双 UP 限定池的 $Pr(S_i)$ 曲线不同，它们的局部极大值出现的位置却是相同的。

单 UP 标准寻访卡池的概率估计值

下面的结果是由两千亿次模拟抽卡仿真得到，即以下面的命令运行程序

./simulation_sequential --standard -t 500000000000 -n 1 -p 50 -c 0

也可以简单地以这条命令运行
./simulation_sequential --standard -t 500000000000 -n 1
效果是一样的。

下面的表格展示了前 $100$ 次抽卡的 $Pr(S_i)$ 与 $Pr(W_i)$ 的数值估计值。

$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$	$i$	$Pr(S_i)$	$Pr(W_i)$
1	1.000%	1.000%	26	0.778%	22.995%	51	0.969%	40.468%	76	0.450%	68.597%
2	0.990%	1.990%	27	0.770%	23.765%	52	1.298%	41.766%	77	0.444%	69.040%
3	0.980%	2.970%	28	0.762%	24.527%	53	1.579%	43.345%	78	0.437%	69.478%
4	0.970%	3.940%	29	0.755%	25.282%	54	1.796%	45.141%	79	0.432%	69.909%
5	0.961%	4.901%	30	0.747%	26.029%	55	1.941%	47.082%	80	0.426%	70.335%
6	0.951%	5.852%	31	0.740%	26.769%	56	2.012%	49.094%	81	0.420%	70.755%
7	0.942%	6.793%	32	0.732%	27.501%	57	2.010%	51.104%	82	0.415%	71.170%
8	0.932%	7.725%	33	0.725%	28.226%	58	1.946%	53.050%	83	0.410%	71.580%
9	0.923%	8.648%	34	0.718%	28.944%	59	1.833%	54.882%	84	0.404%	71.984%
10	0.913%	9.562%	35	0.711%	29.655%	60	1.684%	56.567%	85	0.399%	72.383%
11	0.905%	10.466%	36	0.703%	30.358%	61	1.517%	58.084%	86	0.394%	72.777%
12	0.895%	11.362%	37	0.697%	31.055%	62	1.343%	59.427%	87	0.389%	73.166%
13	0.886%	12.248%	38	0.689%	31.744%	63	1.177%	60.604%	88	0.384%	73.549%
14	0.878%	13.126%	39	0.682%	32.426%	64	1.025%	61.629%	89	0.379%	73.928%
15	0.869%	13.994%	40	0.676%	33.102%	65	0.893%	62.522%	90	0.374%	74.302%
16	0.860%	14.854%	41	0.669%	33.771%	66	0.783%	63.306%	91	0.369%	74.671%
17	0.851%	15.705%	42	0.662%	34.433%	67	0.695%	64.001%	92	0.364%	75.035%
18	0.843%	16.548%	43	0.656%	35.089%	68	0.627%	64.628%	93	0.359%	75.394%
19	0.834%	17.383%	44	0.649%	35.738%	69	0.576%	65.204%	94	0.355%	75.749%
20	0.826%	18.209%	45	0.643%	36.381%	70	0.538%	65.743%	95	0.350%	76.099%
21	0.818%	19.027%	46	0.636%	37.017%	71	0.511%	66.254%	96	0.345%	76.444%
22	0.810%	19.836%	47	0.630%	37.647%	72	0.492%	66.745%	97	0.341%	76.785%
23	0.802%	20.638%	48	0.624%	38.270%	73	0.477%	67.223%	98	0.336%	77.122%
24	0.793%	21.431%	49	0.617%	38.887%	74	0.466%	67.689%	99	0.332%	77.454%
25	0.786%	22.217%	50	0.611%	39.499%	75	0.457%	68.147%	100	0.328%	77.781%

下面的折线图展示了前 $500$ 次抽卡的 $Pr(S_i)$ 与 $Pr(W_i)$ 的值。

蓝色线条代表 $Pr(S_i)$ 的估计值（左侧纵坐标轴），橙色线条代表 $Pr(W_i)$ 的估计值（右侧纵坐标轴）。

下面是一些较为明显的特征：

在第 $56$ 次抽卡中从单 UP 标准寻访卡池中抽到目标 6★ 干员的概率最大；$Pr(S_{56})=2.012\%$ ；
如果你的游戏资源可以让你抽卡 $29$ 次，那么你在这 $29$ 次抽卡之内抽到目标 6★ 干员的概率是 $25\%$ 。
如果你的游戏资源可以让你抽卡 $57$ 次，那么你在这 $57$ 次抽卡之内抽到目标 6★ 干员的概率是 $50\%$ 。
如果你的游戏资源可以让你抽卡 $92$ 次，那么你在这 $92$ 次抽卡之内抽到目标 6★ 干员的概率是 $75\%$ 。
如果你的游戏资源可以让你抽卡 $142$ 次，那么你在这 $142$ 次抽卡之内抽到目标 6★ 干员的概率是 $90\%$ 。

单 UP 标准寻访卡池是获得目标 6★ 干员最容易的卡池类型。虽然获得的 6★ 干员是被概率 UP 的 6★ 干员这个条件概率的值比双 UP 限定寻访卡池的低（$50\%$，限定寻访卡池的是 $70\%$），但是这种卡池只有一位概率 UP 的 6★ 干员，所以这位 6★ 干员能独占这 $50\%$ 的条件概率，而在双 UP 限定寻访卡池中，每个被概率 UP 的 6★ 干员所能获得的条件概率是 $70\% \div 2 = 35\%$ 。

然而，在单 UP 标准寻访卡池中，在 $Pr(S_i)$ 的曲线上只存在两个可被观察到的转折点以及小高峰：

转折点： $i=51,\ 104$
局部极大值及其位置：
- $Pr(S_{56})=2.012\%,\ 当\ i=56\ 时$
- $Pr(S_{114})=0.422\%,\ 当\ i=114\ 时$

这里尝试探讨一下背后的原因。在本种卡池中，显然成功获得目标 6★ 干员的概率要比另外两种卡池要高，也就是说失败的概率会更低。同时，在连续抽卡的过程中，应该存在着一对互为拮抗关系的影响概率的因素，我们分别称它们为 $f(i)$ 与 $g(i)$ ，并有 \[ Pr(S_i)=f(i) \cdot g(i) \] 基于我们对随机事件的定义，$S_i$ 的发生一定意味着随机事件 $F_{i-1}$ 的发生与随机事件 ” 在当前第 $i$ 次抽卡中获得目标 6★ 干员 “ 的发生。同样，这里的 ” 在当前第 $i$ 次抽卡中获得目标 6★ 干员 “ 并不关心前面的抽卡结果，它与随机事件 $S_i$ 的含义是不同的。当 $i$ 增加时，随机事件 $F_{i-1}$ 的概率会降低，因为我们需要让更多次数的抽卡也要失败；随机事件 ” 在当前第 $i$ 次抽卡中获得目标 6★ 干员 “ 会上升，因为随着抽卡次数的增多，将会有 “ 更多的机会 ” 让保底系统开始起作用，和 / 或将获得 6★ 干员的概率提升得更高。所以，当 $i$ 属于某一范围时，将会有下面的式子成立 \[ \begin{aligned} f(i) &= Pr(F_i)\\ g(i) &= Pr(在当前第\ i\ 次抽卡中获得目标\ 6★\ 干员)\\ \Delta f(i) &= f(i+1) - f(i) < 0\\ \Delta g(i) &= g(i+1) - g(i) > 0\ ,\ 当\ i\ 属于某一范围时 \end{aligned} \]

根据仿真得到的结果我们知道，在单 UP 标准寻访卡池中，当 $i$ 来到第二个小高峰后面时，$f(i)$ 下降的速度占据了主导地位，它对 $Pr(S_i)$ 的值的影响盖过了 $g(i)$ 的值上升的影响。所以 $Pr(S_i)$ 便持续下降，以至于第三个小高峰就消失了。然而，如果仔细观察 $Pr(S_i)$ 的曲线的话，拮抗因素之一的 $g(i)$ 在 $161 \le i \le 176$ 附近的确减缓了 $Pr(S_i)$ 的下降速度。

对 $f(i)$ 与 $g(i)$ 的分析也许同样也适用于前两种卡池三个小高峰出现的原因，当然也同样可能适用于为何 $Pr(S_i)$ 的第一个局部极大值出现在 $i=52$ 处而不是 $i=51$ 处，要知道保底系统是在 $i=51$ 就开始起作用的。

需要指出的是，在另一次两千亿次抽卡的仿真中，$Pr(S_{167})$ 与 $Pr(S_{168})$ 的值十分接近，二者差值为 $0.000006\%$ 。第三个小高峰消失的原因也有可能是因为仿真误差所导致的，虽然不太可能真的是因为这个。在我们获得 $Pr(S_i)$ 的代数表达式或理论值之前，我们并不能完全确定的说第三个小高峰是不存在的。这也是为什么我在这部分结果分析中使用了 “ 可被观察到的 ” 一词。

单 UP 限定寻访卡池的概率估计值

这一种卡池目前还没有在明日方舟中出现过，所以此处不展示其概率值。不过你可以查阅 res/limited_500000000000_single_up_simu_1.res 这个文件，里面记录了两千亿次模拟抽卡得到的结果。

你可以通过运行下面的命令来执行这种卡池的仿真

./simulaiton_sequential --limited -t <a-big-enough-number> -n 1 -p 50 -c 0

也可以简单地以这条命令运行
./simulation_sequential -t <a-big-enough-number> -n 1
效果是一样的。

不同种类卡池之间的对比

让我们把不通卡池的概率折线图画在同一个坐标系里：

图中实线是单 UP 标准寻访卡池的 $Pr(S_i)$ 与 $Pr(W_i)$ 值，点线是双 UP 限定寻访卡池的概率值，虚线是双 UP 标准寻访卡池的概率值。可以看出，单 UP 标准寻访卡池的 $Pr(S_i)$ 的值上升得最快，下降得也最快，而双 UP 标准寻访卡池的 $Pr(S_i)$ 的值上升、下降得最慢。

最后，我们把这三种卡池的概率值做个汇总，以便玩家参考。

卡池类型	第一局部极大值 $Pr(S_i)$	第二局部极大值 $Pr(S_i)$	第三局部极大值 $Pr(S_i)$	$\ge 25\%$ 的概率获得你的干员所需抽数	$\ge 50\%$ 的概率获得你的干员所需抽数	$\ge 75\%$ 的概率获得你的干员所需抽数	$\ge 90\%$ 的概率获得你的干员所需抽数
双 UP 限定寻访卡池	$Pr(S_{57})=1.506\%$	$Pr(S_{114})=0.466\%$	$Pr(S_{168})=0.204\%$	$41$ 抽	$70$ 抽	$132$ 抽	$212$ 抽
双 UP 标准寻访卡池	$Pr(S_{57})=1.128\%$	$Pr(S_{114})=0.439\%$	$Pr(S_{169})=0.241\%$	$54$ 抽	$102$ 抽	$187$ 抽	$303$ 抽
单 UP 标准寻访卡池	$Pr(S_{56})=2.012\%$	$Pr(S_{114})=0.422\%$	未观测出	$29$ 抽	$57$ 抽	$92$ 抽	$142$ 抽

希望这些结果对你有帮助。祝欧气满满！$￣︶￣*$)

误差分析

应该存在着更为合理的误差分析的方法。由于我的专业并不是数理统计，此处暂且先以下面这种方法估计仿真结果的误差。如果你知道更为合理的衡量误差的方法，欢迎与我交流，我将十分感激 (￣︶￣*)

在我们研究的问题中，三种卡池的 $Pr(S_1)$ 的理论值是十分容求得的。这里，我将把这个理论值与仿真得到的值作比较，得到 $Pr(S_1)$ 仿真结果的绝对误差 $\Delta Pr(S_1)$ 。接下来假定梅森旋转素数引擎 mt19937_64 所产生的随机数的质量如此之高，以至于所有 $Pr(S_i)$ 仿真结果的绝对误差的绝对值与 $Pr(S_1)$ 仿真结果的绝对误差的绝对值 $|\Delta Pr(S_1)|$ 相等。由于当 $i > 50$ 时获得 $Pr(S_i)$ 的真值（也就是理论值）是较为困难的，所以接下来我们将仿真得到的 $Pr(S_i)$ 都与 $|\Delta Pr(S_1)|$ 相比较，如果有 \[ \frac{|\Delta Pr(S_1)|}{Pr(S_i)} \le 5\% \]

，那么我们就认为仿真得到的这个 $Pr(S_i)$ 是 “ 精确的 ” 。此时，$Pr(W_i)$ 有着与 $Pr(S_i)$ 同样的可以被认为是 “精确的” 的范围，因为 $Pr(W_i)$ 是基于 $Pr(S_i)$ 计算得到的。

此处暂时忽略了误差的积累与传播的问题。

需要注意的是，因为我们并不知道绝大多数 $Pr(S_i)$ 的真值，所以它们的相对误差也是无法求出的。

理论上，在双 UP 限定寻访卡池中， \[ \begin{aligned} Pr(S_1) &= 2\% \times 70\% \div 2\\ &= 0.7\% \end{aligned} \]

而我们的仿真得到的 $Pr(S_1)$ 是 $0.699985\%$（数据详见 res/limited_500000000000_double_up_simu_1.res ）。所以此结果的绝对误差的绝对值为 \[ \begin{aligned} \Delta Pr(S_1) &= 0.700026\% - 0.7\% \\ &= 0.000026\% \\ &= 2.6 \times 10^{-7} \\ | \Delta Pr(S_1) | &= 2.6 \times 10^{-7} \end{aligned} \]

于是仿真结果 $Pr(S_i)$ 可被认为是 “ 精确的 ” 的阈值为 \[ \begin{aligned} |\Delta Pr(S_1)| \div 5\% &= 5.2 \times 10^{-6} \\ &= 0.00052\% \end{aligned} \]

也就是说，$Pr(S_i)$ 需要大于或等于这个阈值，我们才有理由认为这个仿真结果是精确的。按照此判据，在我们的双 UP 限定寻访卡池的仿真数据中，当 $1 \le i \le 677$ 时，仿真得到的 $Pr(S_i)$ 与 $Pr(W_i)$ 是足够精确的。

按照同样的方法，在双 UP 标准寻访卡池的仿真结果（详见 res/standard_500000000000_double_up_simu_1.res ）中，当 $1 \le i \le 780$ 时，仿真得到的 $Pr(S_i)$ 与 $Pr(W_i)$ 是足够精确的。在单 UP 标准寻访卡池的仿真结果（详见 res/standard_500000000000_single_up_simu_1.res ）中，当 $1 i $ 时394，仿真得到的 $Pr(S_i)$ 与 $Pr(W_i)$ 是足够精确的。

如果我们把标准稍微放宽一些，从 $5\%$ 放宽到 $10\%$，则单 UP 标准寻访卡池的仿真数据可以说是精确的范围为 $1 \le i \le 433$，较为接近 $500$，所以我最后仍然把直到前 $500$ 次抽卡的数据都包含在折线图里了。如果我们想要达到更高的景精度，则需要进行更大规模的仿真。

需要指出的是，由于这里采用的误差估计方法可能并不十分科学地反映出仿真结果真正的误差大小，也有可能在上述范围之外更为靠后的仿真结果也是同样精确的。

你也许会发现当 $i$ 十分大时（例如接近 $1000$ 时），仿真得到的 $Pr(S_i)$ 开始明显得上下波动。当我们获得 $Pr(S_i)$ 的理论值之前，我们并不可以完全确定的说这些仿真值是不够精确的——即使它们真的很可能的确不够精确。

将来的工作

由于时间所限，我并没有执行更大规模的蒙特卡洛仿真。目前，所有在本文中出现的数据都是基于五千亿次模拟抽卡仿真得到的，以单线程完成这个规模的模拟，所需时间大约是 $95$ 分钟。在将来我可能以多线程并行的方式来执行更大规模的仿真，以便得到更为精确的结果。

另外，我们也可以尝试使用不同的随机数生成算法。明日方舟是基于 Unity 引擎开发的，我们也许可以尝试 Unity 所使用的随机数生成算法，或者尝试不同的随机种子初始化的方法（例如以玩家的 UID 和 / 或时间来初始化随机种子）。

另外，目前我并没有以不同的 -p 的值来执行仿真。我们其实可以探究在不同的保底系统起点下，$Pr(S_i)$ 的曲线将会发生怎样的变化。这应该也是一个十分有意思的问题——为何明日方舟的游戏设计者们选择保底系统在第 $50$ 次抽卡之后开始起作用，而不是第 $20$ 或者是第 $30$ 次？

最后，我们可以有一种自动化的方法基于仿真结果去生成这些折线统计图，以进行数据可视化。我也许会用 Python 实现这一功能。

结语与声明

由于我并不是数理统计专业，本文只能做一些深度有限的讨论。如果你知道计算 $Pr(S_i)$ 与 $Pr(F_i)$ 的理论值的方法，或者进行误差分析的更为合理的方法，或者发现本文中的一些错误，非常欢迎与我讨论，我将十分感激。事实上我对$Pr(S_i)$ 与 $Pr(F_i)$ 的理论值的计算方法十分感兴趣，欢迎发帖交流😀！
感谢我的朋友浮生之夏(@jlgg)，他与我进行了许多理论推导部分的讨论，并指出了许多在我刚接触马克科夫模型时所犯的错误。
目前已经有一些非常棒的代码与文章谈探究过类似的问题，它们给了我很多启发——

💮Frizu - Arknights 6-star Headhunting Rate Analysis：在这个项目中作者对在明日方舟中连续抽卡时获得任意 6★ 干员的概率做了大量研究。Frizu 有一系列文章讲解了他所使用的方法、背后的数学原理以及所得到的结果：Arknights 6★ Headhunting Rate Analysis, Arknights’ Non-6★ Headhunting Streak EDA

💮GalAster - ArkKnightsCapsule：在这个项目中，作者估计了在双 UP 标准寻访卡池与双 UP 限定寻访卡池中，获得概率 UP 的其中一个 6★ 干员以及两个概率 UP 的 6★ 干员都获得的概率，另外，作者也同样估计了让出率 UP 6★ 干员满潜（即连续 $6$ 次获得同一 6★ 干员）的概率。

卡池类型	第一局部极大值 \(Pr(S_i)\)	第二局部极大值 \(Pr(S_i)\)	第三局部极大值 \(Pr(S_i)\)	\(\ge 25\%\) 的概率获得你的干员所需抽数	\(\ge 50\%\) 的概率获得你的干员所需抽数	\(\ge 75\%\) 的概率获得你的干员所需抽数	\(\ge 90\%\) 的概率获得你的干员所需抽数
双 UP 限定寻访卡池	\(Pr(S_{57})=1.506\%\)	\(Pr(S_{114})=0.466\%\)	\(Pr(S_{168})=0.204\%\)	\(41\) 抽	\(70\) 抽	\(132\) 抽	\(212\) 抽
双 UP 标准寻访卡池	\(Pr(S_{57})=1.128\%\)	\(Pr(S_{114})=0.439\%\)	\(Pr(S_{169})=0.241\%\)	\(54\) 抽	\(102\) 抽	\(187\) 抽	\(303\) 抽
单 UP 标准寻访卡池	\(Pr(S_{56})=2.012\%\)	\(Pr(S_{114})=0.422\%\)	未观测出	\(29\) 抽	\(57\) 抽	\(92\) 抽	\(142\) 抽

明日方舟特定六星干员出率理论值的蒙特卡洛仿真估计

zyLiu6707

2021/1/8

什么是明日方舟?

使用本项目

构建代码

命令行参数

程序额外说明

例子

明日方舟中的数学模型

明日方舟抽卡规则与卡池类型

随机事件定义以及所研究的问题

6★ 干员出率：概率与期望值

目标 6★ 干员出率：获得概率值的尝试与方法

也许可行的获得理论值的方法

本项目所采取的方法：蒙特卡洛方法

随机数的生成

不要再用 `rand() % n` 了！

推荐的方法 - C++11 `<random>`库

仿真结果分析

双 UP 限定寻访卡池的概率估计值

双 UP 标准寻访卡池的概率估计值

单 UP 标准寻访卡池的概率估计值

单 UP 限定寻访卡池的概率估计值

不同种类卡池之间的对比

误差分析

将来的工作

结语与声明

明日方舟特定六星干员出率理论值的蒙特卡洛仿真估计

zyLiu6707

2021/1/8

什么是明日方舟?

使用本项目

构建代码

命令行参数

程序额外说明

例子

明日方舟中的数学模型

明日方舟抽卡规则与卡池类型

随机事件定义以及所研究的问题

6★ 干员出率：概率与期望值

目标 6★ 干员出率：获得概率值的尝试与方法

也许可行的获得理论值的方法

本项目所采取的方法：蒙特卡洛方法

随机数的生成

不要再用 rand() % n 了！

推荐的方法 - C++11 <random>库

仿真结果分析

双 UP 限定寻访卡池的概率估计值

双 UP 标准寻访卡池的概率估计值

单 UP 标准寻访卡池的概率估计值

单 UP 限定寻访卡池的概率估计值

不同种类卡池之间的对比

误差分析

将来的工作

结语与声明

不要再用 `rand() % n` 了！

推荐的方法 - C++11 `<random>`库