通过

2018 年 2 月

第 33 卷，第 2 期

测试运行 - 使用 C# 的汤普森采样

作者 James McCaffrey

James McCaffrey汤普森采样是一种算法，可用于解决多臂老虎机问题。这个词就是从赌博机的绰号“单臂老虎机”引申而来。

假设面前有三台赌博机。如果拉动其中一台赌博机的单臂，它可能会不吐钱，也可能会吐出一美元，而你并不知道赌博机吐钱的概率分布是什么。例如，假设赌博机的平均吐钱概率为 0.5，那么如果拉动赌博机的单臂 100 次，赌博机不吐钱和吐出一美元的次数应分别约为 50 次。

我们的目标是尽可能高效地确定吐钱概率最大的赌博机。刚刚介绍的便是伯努利分布老虎机问题，因为结果不是 0 就是 1。还有其他类型的老虎机问题。例如，如果每台赌博机根据钟形高斯分布吐出的金额介于零和一美元之间（如 55.0 美分），遇到的就是高斯过程老虎机问题。本文仅涉及伯努利分布老虎机问题。

若要尝试确定吐钱概率最大的赌博机，可以使用多种算法。假设总共只能拉动三台赌博机的单臂 100 次。可以尝试拉动每台赌博机的单臂 10 次，同时记录每台赌博机的吐钱情况。然后，可以确定在 30 次拉动探索阶段吐钱最多的赌博机，并将剩余的 70 次拉动机会全都用到这台赌博机上。这种方法的危险之处在于，可能会错误地确定吐钱概率最大的赌博机。如果在探索阶段进行许多次拉动，那么在开发利用阶段所剩的拉动次数就不多了。

汤普森采样是一种非常巧妙的算法，可不断调整每台赌博机吐钱的估计概率。大家很快就会发现，使用汤普森采样解决伯努利分布老虎机问题的关键在于，数学上的 beta 分布。

虽然你的老板不太可能会要求你分析赌博机，但在现实生活中的许多重要场合下都可能会遇到多臂老虎机问题。例如，假设你就职于一家药品生产企业。你刚刚研发出了四种新的实验性药物用于治疗癌症，并且希望确定这四种药物中哪一种最有效，同时最大限度地减少活体测试次数。还有可能是你就职于一家电子商务公司，想要确定几个新广告活动中哪一个最成功。  

了解本文所述观点的一个好方法是查看图 1 中的演示程序。此演示程序设置了三个伯努利分布赌博机，吐钱概率分别为 (0.3, 0.5, 0.7)。当然，在非演示情况下，概率是未知的。总共只能拉动 10 次。目标是确定吐钱概率最大的赌博机（3 号赌博机），并确保它的单臂拉动次数最多。

汤普森采样演示
图 1：汤普森采样演示

在第一次试算中，假设三台赌博机的吐钱概率全都为 0.5。演示程序使用 beta 分布并根据此假设生成三个概率。这些随机概率为 (0.7711, 0.1660, 0.1090)。由于与 1 号赌博机相关的概率为最高，因此拉动的是它的单臂。此时，1 号赌博机并未吐钱。

在第二次试算中，演示程序在后台认为第一台赌博机的吐钱概率现在小于 0.5。使用的是 beta 采样，这一次概率分别为 (0.6696, 0.2250, 0.7654)。因此，选择和拉动的是 3 号赌博机的单臂，并且它的估计吐钱概率会根据赌博机是否吐钱进行调整。

在一段时间后，由于 3 号赌博机的吐钱概率最高，因此它胜出的次数也多于另外两台赌博机；并且只要 3 号赌博机吐钱，它被选择的可能性也会随之增加。

10 次试算之后，1 号赌博机的单臂被拉动了三次，只有一次吐钱，因此模拟猜测的真实吐钱概率大约为 1/3 = 0.33。2 号赌博机的单臂被拉动了三次，并吐钱两次（估计概率为 p = 2/3 = 0.6667）。3 号赌博机的单臂被拉动了四次，并吐钱三次（估计概率为 p = 3/4 = 0.7500）。所以，至少在此示例中，不仅确定了吐钱概率最大的赌博机，而且它的单臂拉动次数也最多。

若要更好地理解本文，至少必须拥有中等水平的编程技能，但无需对汤普森采样或 beta 分布有任何了解。虽然演示程序是使用 C# 进行编码，但也可以将代码重构为其他语言（如 Visual Basic 或 Python），应该不会遇到太多麻烦。本文提供了演示程序的全部代码，也可以通过下载随附的文件获取代码。

了解 beta 分布

汤普森采样能否解决伯努利分布老虎机问题取决于是否使用 beta 分布。通常必须了解概率分布，才能理解 beta 分布。概率分布有许多类型，每种分布都有依赖一两个参数的变体。

大家非常熟悉的可能是均匀分布，它依赖两个参数，分别称为 min 和 max（有时也直接称为 a 和 b）。如果是 min = 0.0 和 max = 1.0 的均匀分布，分布将返回介于 0.0 和 1.0 之间的 p 值，其中每个值的返回可能性均相等。因此，如果通过均匀分布采样 1,000 次，应该可以获取约 100 个介于 0.00 和 0.10 之间的 p 值、约 100 个介于 0.10 和 0.20 之间的 p 值（以此类推...直到约 100 个介于 0.90 和 1.00 之间的 p 值）。如果将结果绘制成图，便会看到一个条形图，其中包含 10 个条形且高度全都大致相同。

大家可能还对正态分布（亦称为“高斯分布”）非常熟悉。正态分布的特点也是依赖两个参数，即平均偏差和标准偏差。如果通过平均偏差 = 0.0 且标准偏差 = 1.0 的正态分布采样 1,000 次，应该可以获取约 380 个介于 -0.5 和 +0.5 之间的 z 值、约 240 个介于 +0.5 和 +1.5 之间的 z 值（-0.5 到 -1.5 区间亦是如此）、约 60 个介于 +1.5 和 +2.5 之间的 z 值（-1.5 到 -2.5 区间亦是如此）、约 10 个大于 +2.5 的 z 值（以及不到 10 个小于 -2.5 的 z 值）。如果将结果绘制成图，便会看到一个钟形条形图。

beta 分布的特点是依赖两个参数，通常称为 alpha 和 beta（有时也直接称为 a 和 b）。请注意区分表示整个分布的“beta”和表示两个分布参数中第二个的“beta”，以免发生混淆。

如果通过 a = 1 且 b = 1 的 beta 分布进行采样，结果与平均值 = 0.5 的均匀分布完全相同。如果通过 a 和 b 值不同的 beta 分布进行采样，p 值的平均值为 a / (a+b)。例如，如果 a = 3 且 b = 1，那么在重复采样后将获得介于 0.0 和 1.0 之间的 p 值，但 p 值的平均值为 3 / (3+1) = 0.75。也就是说，介于 0.90 和 1.00 之间的 p 值最常见；介于 0.80 和 0.90 之间的 p 值有些不太常见（以此类推...直到介于 0.00 和 0.10 之间的 p 值几乎没有）。图 2 展示了通过 a = 3 且 b = 1 的 beta 分布进行 10,000 次采样的结果。

通过 Beta(3,1) 分布进行采样
图 2：通过 Beta(3,1) 分布进行采样

beta 分布和伯努利分布老虎机问题之间的联系相当微妙。简而言之，beta 分布是伯努利分布似然函数的共轭先验。尽管这其中蕴含的数学知识非常深奥，但幸运的是，汤普森算法的实现（相对）简单。

演示程序

为了创建演示程序，我启动了 Visual Studio，并新建了名为“汤普森”的控制台应用。由于演示程序没有明显的 .NET Framework 依赖关系，因此可以使用任何一版 Visual Studio。将模板代码加载到编辑器窗口后，我右键单击了文件 Program.cs，并将它重命名为更具描述性一点的 ThompsonProgram.cs，同时允许 Visual Studio 自动为我重命名类 Program。在模板代码的顶部，我删除了对不需要命名空间的所有引用，只留下了对顶级 System 命名空间的引用。

图 3 展示了整体程序结构（为节省空间，进行了少量小幅改动）。所有控制逻辑都位于 Main 方法中。beta 分布采样是使用程序定义的 BetaSampler 类进行实现。为了节省空间，我删除了所有常规错误检查。

图 3：汤普森采样演示程序结构

using System;
namespace Thompson
{
  class ThompsonProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("Begin Thompson sampling demo");
      int N = 3;
      double[] means = new double[] { 0.3, 0.5, 0.7 };
      double[] probs = new double[N];
      int[] S = new int[N];
      int[] F = new int[N];
      Random rnd = new Random(4);
      BetaSampler bs = new BetaSampler(2);
      for (int trial = 0; trial < 10; ++trial)
      {
        Console.WriteLine("Trial " + trial);
        for (int i = 0; i < N; ++i)
          probs[i] = bs.Sample(S[i] + 1.0, F[i] + 1.0);
        Console.Write("sampling probs = ");
        for (int i= 0; i < N; ++i)
          Console.Write(probs[i].ToString("F4") + " ");
        Console.WriteLine("");
        int machine = 0;
        double highProb = 0.0;
        for (int i = 0; i < N; ++i) {
          if (probs[i] > highProb) {
            highProb = probs[i];
            machine = i;
          }
        }
        Console.Write("Playing machine " + machine);
        double p = rnd.NextDouble();
        if (p < means[machine]) {
          Console.WriteLine(" -- win");
          ++S[machine];
        }
        else {
          Console.WriteLine(" -- lose");
          ++F[machine];
        }
      }
      Console.WriteLine("Final estimates of means: ");
      for (int i = 0; i < N; ++i)  {
        double u = (S[i] * 1.0) / (S[i] + F[i]);
        Console.WriteLine(u.ToString("F4") + "  ");
      }
      Console.WriteLine("Number times machine played:");
      for (int i = 0; i < N; ++i) {
        int ct = S[i] + F[i];
        Console.WriteLine(ct + "  ");
      }
      Console.WriteLine("End demo ");
      Console.ReadLine();
    }
  }
  public class BetaSampler
  {
    // ...
  }
}

首先，演示程序设置四个数组：

Console.WriteLine("Begin Thompson sampling demo");
int N = 3;
double[] means = new double[] { 0.3, 0.5, 0.7 };
double[] probs = new double[N];
int[] S = new int[N];
int[] F = new int[N];

虽然演示程序演示的是三台赌博机，但汤普森采样支持任意数量的赌博机。平均吐钱概率采用硬编码。平均概率彼此越接近，问题的难度就越大。probs 数组保留了 beta 分布采样生成的概率，以确定拉动哪台赌博机的单臂。S（“成功”）数组和 F（“失败”）数组保留了每台赌博机在单臂被拉动时吐钱和不吐钱的累计次数。

演示程序使用两个随机生成数字的对象：

Random rnd = new Random(4);
BetaSampler bs = new BetaSampler(2);

rnd 对象用于确定选定赌博机是胜出还是未胜出。请注意，我交替使用了胜出和未胜出、吐钱和不吐钱、成功和失败这三组词语。bs 对象用于执行 beta 分布采样。之所以指定了使用的种子 2 和 4，只是因为它们提供了代表性演示运行。 

主模拟循环开始：

for (int trial = 0; trial < 10; ++trial) {
  Console.WriteLine("Trial " + trial);
  for (int i = 0; i < N; ++i)
    probs[i] = bs.Sample(S[i] + 1.0, F[i] + 1.0);
...

建议将试算次数设置为较大的数字（如 1,000），可以看到汤普森采样对效果最佳的赌博机快速归零。整个演示程序的关键在于，调用 Sample 方法。成功和失败次数都会传递到此方法。由于 beta 分布要求 a 和 b 参数均大于 0，因此添加了常数 1.0 作为“入侵”内容。如果以前就了解过一点赌博机的吐钱概率，可以将此常数调整为反映相应信息。

显示更新后的采样概率后，演示程序选择采样概率最高的赌博机：

int machine = 0;
double highProb = 0.0;
for (int i = 0; i < N; ++i) {
  if (probs[i] > highProb) {
    highProb = probs[i];
    machine = i;
  }
}

由于是概率采样，因此即使赌博机一次也没有胜出，且失败过许多次，也仍可能会选择拉动它的单臂。这种机制增强了算法的探索功能。

最后，主迭代循环拉动选定赌博机的单臂：

...
  Console.Write("Playing machine " + machine);
  double p = rnd.NextDouble();
  if (p < means[machine]) {
    Console.WriteLine("win"); ++S[machine];
  }
  else {
    Console.WriteLine("lose"); ++F[machine];
  }
}

NextDouble 方法返回介于 0.0 和 1.0 之间的均匀随机值，用于实现伯努利分布过程。最后，演示程序会显示每台赌博机的估计吐钱概率（不用费心检查是否有除数为零），以及每台赌博机的单臂拉动次数。

实现 beta 分布

有些令人惊奇的是，据我所知，.NET Framework 没有库方法可以执行 beta 分布采样。我决定从头开始实现一个，而不是依赖外部库。用于生成 beta 示例的算法有很多。我使用了 R. C. H. Cheng 于 1978 年开发的“BA”算法。图 4 展示了完整代码。

图 4：程序定义的 beta 分布取样器

public class BetaSampler
{
  public Random rnd;
  public BetaSampler(int seed)
  {
    this.rnd = new Random(seed);
  }
  public double Sample(double a, double b)
  {
    double alpha = a + b;
    double beta = 0.0;
    double u1, u2, w, v = 0.0;
    if (Math.Min(a, b) <= 1.0)
      beta = Math.Max(1 / a, 1 / b);
    else
      beta = Math.Sqrt((alpha - 2.0) /
(2 * a * b - alpha));
    double gamma = a + 1 / beta;
    while (true) {
      u1 = this.rnd.NextDouble();
      u2 = this.rnd.NextDouble();
      v = beta * Math.Log(u1 / (1 - u1));
      w = a * Math.Exp(v);
      double tmp = Math.Log(alpha / (b + w));
      if (alpha * tmp + (gamma * v) - 1.3862944 >=
 Math.Log(u1 * u1 * u2))
        break;
    }
    double x = w / (b + w);
    return x;
  }
}

beta 分布采样本身就是一个非常有吸引力的主题。快速扫一眼代码，应该就会相信，确实有某种十分巧妙的非实验性数学算法。此实现非常遵循源代码研究论文，包括使用相同的变量名称等等。请注意可能的无限循环，这对于研究伪代码很常见，但对于生产代码，这是绝对禁止的。建议添加循环计数器变量，用于在值超过设定的阈值时抛出异常。

总结

若要尝试使用汤普森采样解决伯努利分布多臂老虎机问题，请参阅本文及其代码，其中应该包含了所需的全部信息。借助其中的信息，还应该可以探索使用其他类型的吐钱函数来解决老虎机问题。例如，如果有赌博机是根据 Poisson 似然函数吐钱，将会通过 gamma 分布（即 Poisson 的共轭先验）进行采样，而不是使用 beta 分布。

多臂老虎机问题便是所谓的加强学习 (RL)。在 RL 机器学习中，不是使用已知正确答案的标记数据创建预测系统，而是使用某种回报函数快速生成预测模型。RL 是机器学习研究的重心。

Dr.James McCaffrey 供职于华盛顿地区雷蒙德市沃什湾的 Microsoft Research。他参与过多个 Microsoft 产品的工作，包括 Internet Explorer 和必应。Scripto可通过 jamccaff@microsoft.com 与 McCaffrey 取得联系。

衷心感谢以下 Microsoft 技术专家对本文的审阅：Chris Lee、Ricky Loynd、Adith Swaminathan

在 MSDN 杂志论坛讨论这篇文章