贝叶斯A/B测试简介

from dataclasses import dataclass
from typing import Dict, List, Union

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm

from scipy.stats import bernoulli, expon

RANDOM_SEED = 4000
rng = np.random.default_rng(RANDOM_SEED)

%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")

plotting_defaults = dict(
    bins=50,
    kind="hist",
    textsize=10,
)

本笔记本演示了如何实现A/B测试的贝叶斯分析。我们实现了VWO的贝叶斯A/B测试白皮书中讨论的模型[Stucchio, 2015]，并讨论了这些模型中不同先验选择的影响。本笔记本不讨论其他相关主题，如如何选择先验、早期停止和功效分析。

什么是A/B测试？

来自 https://vwo.com/ab-testing/:

A/B 测试（也称为拆分测试）是一种将同一网页的两个变体同时展示给不同网站访问者群体的过程，并比较哪个变体带来更多的转化。

具体来说，A/B 测试常用于软件行业，以确定是否应向用户发布新功能或对现有功能的更改，以及这些更改对核心产品指标（“转化率”）的影响。此外：

• 我们可以同时测试多个变体。我们还将在这个笔记本中处理如何分析这些测试。

• “转化”的确切含义在不同测试中可能有所不同，但我们关注的两种转化类别是：

  • 伯努利转化 - 一个标志，用于表示访问者是否执行了目标操作（例如，至少完成了一次购买）。

  • 价值转换 - 每位访客的实际价值（例如，美元收入，也可能为0）。

如果你在生物学、心理学和其他科学领域中学习过控制实验，A/B测试听起来会非常像一个控制实验——因为它的确如此！控制组和处理组的概念，以及实验设计原则，是A/B测试的基石。主要区别在于实验运行的环境：A/B测试通常由在线软件公司运行，实验对象是访问网站/应用的用户，感兴趣的结果是可追踪的行为，如注册、购买产品以及返回网站。

A/B测试通常使用传统的假设检验进行分析（参见t检验），但另一种方法是使用贝叶斯统计。这使我们能够结合先验分布，并生成一系列关于“是否存在获胜变体？”和“差异有多大？”的问题的答案。

伯努利转换

让我们首先处理一个简单的两变体A/B测试，其中感兴趣的指标是执行操作的用户比例（例如，至少购买一件商品），这是一个伯努利转化。我们的变体称为A和B，其中A指的是现有的登录页面，B指的是我们想要测试的新页面。我们想要进行统计推断的结果是B是否“优于”A，这取决于每个变体的潜在“真实”转化率。我们可以将其表述如下：

设 \(\theta_A, \theta_B\) 分别为变体 A 和 B 的真实转化率。那么，访问者在变体 A 中是否转化的结果是随机变量 \(\mathrm{Bernoulli}(\theta_A)\)，而变体 B 的结果是 \(\mathrm{Bernoulli}(\theta_B)\)。如果我们假设访问者在落地页上的行为与其他访问者无关（这是一个合理的假设），那么变体的总转化数 \(y\) 服从二项分布：

\[y \sim \sum^N\mathrm{Bernoulli}(\theta) = \mathrm{Binomial}(N, \theta)\]

在贝叶斯框架下，我们假设真实的转化率\(\theta_A, \theta_B\)是未知的，并且它们各自遵循Beta分布。假设潜在的转化率是独立的（我们会随机分配流量到每个变体之间，因此一个变体不会影响另一个变体）：

\[\theta_A \sim \theta_B \sim \mathrm{Beta}(\alpha, \beta)\]

A/B测试期间观察到的数据（似然分布）是：访问页面的访客数量 N，以及至少购买了一件商品的访客数量 y：

\[y_A \sim \mathrm{Binomial}(n = N_A, p = \theta_A), y_B \sim \mathrm{Binomial}(n = N_B, p = \theta_B)\]

通过这个，我们可以从\(\theta_A, \theta_B\)的联合后验分布中采样。

您可能已经注意到，Beta 分布是二项分布的共轭先验，因此我们不需要 MCMC 采样来估计后验（VWO 论文中可以找到精确解）。我们仍然会演示如何使用 PyMC 进行采样，这样做可以更容易地扩展模型，使用不同的先验、依赖假设等。

最后，请记住我们感兴趣的结果是B是否优于A。在实践中，判断B是否优于A的一个常见指标是转化率的相对提升，即\(\theta_B\)相对于\(\theta_A\)的百分比差异：

\[\mathrm{reluplift}_B = \theta_B / \theta_A - 1\]

我们将在下面使用PyMC实现这个模型设置。

@dataclass
class BetaPrior:
    alpha: float
    beta: float

@dataclass
class BinomialData:
    trials: int
    successes: int

class ConversionModelTwoVariant:
    def __init__(self, priors: BetaPrior):
        self.priors = priors

    def create_model(self, data: List[BinomialData]) -> pm.Model:
        trials = [d.trials for d in data]
        successes = [d.successes for d in data]
        with pm.Model() as model:
            p = pm.Beta("p", alpha=self.priors.alpha, beta=self.priors.beta, shape=2)
            obs = pm.Binomial("y", n=trials, p=p, shape=2, observed=successes)
            reluplift = pm.Deterministic("reluplift_b", p[1] / p[0] - 1)
        return model

现在我们已经定义了一个可以接受先验和合成数据作为输入的类，我们的第一步是选择一个合适的先验。在实践中进行此操作时需要考虑几个因素，但为了本笔记本的目的，我们将重点关注以下内容：

• 我们假设为每个变体设置了相同的Beta先验。

• 当我们将alpha和beta的值设置得较低时，会出现无信息或弱信息先验。例如，alpha = 1, beta = 1会导致均匀分布作为先验。如果我们单独考虑一个分布，设置这个先验意味着我们对该参数的值一无所知，也不了解其周围的置信度。然而，在A/B测试的背景下，我们感兴趣的是比较一个变体相对于另一个变体的相对提升。使用弱信息Beta先验，这种相对提升分布非常宽泛，因此我们隐含地表示变体之间可能存在很大差异。

• 当我们将alpha和beta的值设置得较高时，会出现一个强先验。与上述情况相反，强先验意味着相对提升分布是薄的，即我们的先验信念是变体之间没有太大差异。

我们通过先验预测检查来说明这些点。

先验预测检查

请注意，在这些先验预测检查中，我们可以为观测数据传递任意值。PyMC在从先验预测分布中采样时不会使用这些数据。

weak_prior = ConversionModelTwoVariant(BetaPrior(alpha=100, beta=100))

strong_prior = ConversionModelTwoVariant(BetaPrior(alpha=10000, beta=10000))

with weak_prior.create_model(data=[BinomialData(1, 1), BinomialData(1, 1)]):
    weak_prior_predictive = pm.sample_prior_predictive(samples=10000, return_inferencedata=False)

with strong_prior.create_model(data=[BinomialData(1, 1), BinomialData(1, 1)]):
    strong_prior_predictive = pm.sample_prior_predictive(samples=10000, return_inferencedata=False)

fig, axs = plt.subplots(2, 1, figsize=(7, 7), sharex=True)
az.plot_posterior(weak_prior_predictive["reluplift_b"], ax=axs[0], **plotting_defaults)
axs[0].set_title(f"B vs. A Rel Uplift Prior Predictive, {weak_prior.priors}", fontsize=10)
axs[0].axvline(x=0, color="red")
az.plot_posterior(strong_prior_predictive["reluplift_b"], ax=axs[1], **plotting_defaults)
axs[1].set_title(f"B vs. A Rel Uplift Prior Predictive, {strong_prior.priors}", fontsize=10)
axs[1].axvline(x=0, color="red");

在使用弱先验的情况下，B相对于A的相对提升的94% HDI大约为[-20%, +20%]，而在使用强先验的情况下，它大约为[-2%, +2%]。这实际上是相对提升分布的“起点”，并将影响观察到的转化如何转化为后验分布。

我们在实践中如何选择这些先验取决于运行A/B测试的公司的更广泛背景。一个强先验可以帮助防止错误发现，但当存在获胜变体时，可能需要更多数据来检测它们（更多数据 = 需要更多时间来运行测试）。一个弱先验给予观测数据更多权重，但也可能导致更多错误发现，这是由于早期停止问题造成的。

下面我们将详细介绍两种不同先验选择的推理结果。

数据

我们生成了两个数据集：一个数据集中每个变体的“真实”转化率相同，另一个数据集中变体B的真实转化率更高。

def generate_binomial_data(
    variants: List[str], true_rates: List[str], samples_per_variant: int = 100000
) -> pd.DataFrame:
    data = {}
    for variant, p in zip(variants, true_rates):
        data[variant] = bernoulli.rvs(p, size=samples_per_variant)
    agg = (
        pd.DataFrame(data)
        .aggregate(["count", "sum"])
        .rename(index={"count": "trials", "sum": "successes"})
    )
    return agg

# Example generated data
generate_binomial_data(["A", "B"], [0.23, 0.23])

	A	B
trials	100000	100000
successes	22979	22970

我们还将编写一个函数来封装数据生成、采样和后验图，以便我们可以轻松比较在两种情况下（相同真实速率与不同真实速率）两种模型（强先验和弱先验）的结果。

def run_scenario_twovariant(
    variants: List[str],
    true_rates: List[float],
    samples_per_variant: int,
    weak_prior: BetaPrior,
    strong_prior: BetaPrior,
) -> None:
    generated = generate_binomial_data(variants, true_rates, samples_per_variant)
    data = [BinomialData(**generated[v].to_dict()) for v in variants]
    with ConversionModelTwoVariant(priors=weak_prior).create_model(data):
        trace_weak = pm.sample(draws=5000)
    with ConversionModelTwoVariant(priors=strong_prior).create_model(data):
        trace_strong = pm.sample(draws=5000)

    true_rel_uplift = true_rates[1] / true_rates[0] - 1

    fig, axs = plt.subplots(2, 1, figsize=(7, 7), sharex=True)
    az.plot_posterior(trace_weak.posterior["reluplift_b"], ax=axs[0], **plotting_defaults)
    axs[0].set_title(f"True Rel Uplift = {true_rel_uplift:.1%}, {weak_prior}", fontsize=10)
    axs[0].axvline(x=0, color="red")
    az.plot_posterior(trace_strong.posterior["reluplift_b"], ax=axs[1], **plotting_defaults)
    axs[1].set_title(f"True Rel Uplift = {true_rel_uplift:.1%}, {strong_prior}", fontsize=10)
    axs[1].axvline(x=0, color="red")
    fig.suptitle("B vs. A Rel Uplift")
    return trace_weak, trace_strong

场景1 - 相同的潜在转化率

trace_weak, trace_strong = run_scenario_twovariant(
    variants=["A", "B"],
    true_rates=[0.23, 0.23],
    samples_per_variant=100000,
    weak_prior=BetaPrior(alpha=100, beta=100),
    strong_prior=BetaPrior(alpha=10000, beta=10000),
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]

100.00% [24000/24000 00:09<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 21 seconds.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]

100.00% [24000/24000 00:09<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 20 seconds.

• 在这两种情况下，真正的提升值0%都位于94%的高密度区间内。

• 然后，我们可以使用这个相对提升分布来决定是否将变体B中的新登录页面/功能作为默认设置。例如，我们可以决定如果94%的高密度区间（HDI）高于0，我们将推出变体B。在这种情况下，0在HDI中，因此决定是不推出变体B。

场景2 - 不同的基础利率

run_scenario_twovariant(
    variants=["A", "B"],
    true_rates=[0.21, 0.23],
    samples_per_variant=100000,
    weak_prior=BetaPrior(alpha=100, beta=100),
    strong_prior=BetaPrior(alpha=10000, beta=10000),
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]

100.00% [24000/24000 00:09<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 21 seconds.
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]

100.00% [24000/24000 00:09<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 20 seconds.

(Inference data with groups:
 	> posterior
 	> log_likelihood
 	> sample_stats
 	> observed_data,
 Inference data with groups:
 	> posterior
 	> log_likelihood
 	> sample_stats
 	> observed_data)

• 在这两种情况下，后验相对提升分布表明B的转化率高于A，因为94% HDI远高于0。在这种情况下，决策是将变体B推广给所有用户，这个结果是“真正的发现”。

• 尽管如此，在实践中我们通常也会对变体B有多好感兴趣。对于具有强先验的模型，先验实际上将相对提升分布拉近到0，因此我们对相对提升的中心估计是保守的（即低估的）。我们需要更多的数据才能使我们的推断更接近9.5%的真实相对提升。

上述示例展示了如何使用简单的Beta-Binomial模型对两变量测试进行A/B测试分析，以及选择弱先验与强先验的优缺点。在下一节中，我们将提供处理多变量（“A/B/n”）测试的指南。

推广到多变量测试

我们将继续在本节中使用伯努利转换和Beta-二项式模型，以简化问题。重点是如何分析具有3个或更多变体的测试——例如，我们不仅仅有一个不同的登录页面进行测试，而是有多个想法想要同时测试。我们如何判断其中是否有赢家？

这里我们有两个主要的方法可以采用：

1. 将A作为‘对照组’。逐一将其他变体（B、C等）与A进行比较。

2. 对于每个变体，与其他变体的max()进行比较。

方法1对大多数人来说很直观，也容易解释。但如果有两个变体都优于对照组，而我们想知道哪一个更好呢？我们无法通过单独的提升分布来做出这个推断。方法2确实处理了这种情况——它有效地尝试找出所有变体中是否存在一个明显的赢家或明显的输家。

我们将在下面为这两种方法实现模型设置，清理我们之前的代码，使其能够推广到n变体的情况。请注意，我们也可以为2变体情况重新使用此模型。

class ConversionModel:
    def __init__(self, priors: BetaPrior):
        self.priors = priors

    def create_model(self, data: List[BinomialData], comparison_method) -> pm.Model:
        num_variants = len(data)
        trials = [d.trials for d in data]
        successes = [d.successes for d in data]
        with pm.Model() as model:
            p = pm.Beta("p", alpha=self.priors.alpha, beta=self.priors.beta, shape=num_variants)
            y = pm.Binomial("y", n=trials, p=p, observed=successes, shape=num_variants)
            reluplift = []
            for i in range(num_variants):
                if comparison_method == "compare_to_control":
                    comparison = p[0]
                elif comparison_method == "best_of_rest":
                    others = [p[j] for j in range(num_variants) if j != i]
                    if len(others) > 1:
                        comparison = pm.math.maximum(*others)
                    else:
                        comparison = others[0]
                else:
                    raise ValueError(f"comparison method {comparison_method} not recognised.")
                reluplift.append(pm.Deterministic(f"reluplift_{i}", p[i] / comparison - 1))
        return model

def run_scenario_bernoulli(
    variants: List[str],
    true_rates: List[float],
    samples_per_variant: int,
    priors: BetaPrior,
    comparison_method: str,
) -> az.InferenceData:
    generated = generate_binomial_data(variants, true_rates, samples_per_variant)
    data = [BinomialData(**generated[v].to_dict()) for v in variants]
    with ConversionModel(priors).create_model(data=data, comparison_method=comparison_method):
        trace = pm.sample(draws=5000)

    n_plots = len(variants)
    fig, axs = plt.subplots(nrows=n_plots, ncols=1, figsize=(3 * n_plots, 7), sharex=True)
    for i, variant in enumerate(variants):
        if i == 0 and comparison_method == "compare_to_control":
            axs[i].set_yticks([])
        else:
            az.plot_posterior(trace.posterior[f"reluplift_{i}"], ax=axs[i], **plotting_defaults)
        axs[i].set_title(f"Rel Uplift {variant}, True Rate = {true_rates[i]:.2%}", fontsize=10)
        axs[i].axvline(x=0, color="red")
    fig.suptitle(f"Method {comparison_method}, {priors}")

    return trace

我们生成数据，其中变体B和C远高于A，但彼此非常接近：

_ = run_scenario_bernoulli(
    variants=["A", "B", "C"],
    true_rates=[0.21, 0.23, 0.228],
    samples_per_variant=100000,
    priors=BetaPrior(alpha=5000, beta=5000),
    comparison_method="compare_to_control",
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]

100.00% [24000/24000 00:09<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 22 seconds.

• B和C的相对提升后验概率表明，它们明显优于A（94% HDI远高于0），相对提升约为7-8%。

• 然而，我们无法推断B和C之间是否有胜者。

_ = run_scenario_bernoulli(
    variants=["A", "B", "C"],
    true_rates=[0.21, 0.23, 0.228],
    samples_per_variant=100000,
    priors=BetaPrior(alpha=5000, beta=5000),
    comparison_method="best_of_rest",
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [p]

100.00% [24000/24000 00:10<00:00 Sampling 4 chains, 0 divergences]

Sampling 4 chains for 1_000 tune and 5_000 draw iterations (4_000 + 20_000 draws total) took 22 seconds.

• A的提升图告诉我们，与变体B和C相比，它明显处于劣势（A的相对提升的94% HDI远低于0）。

• 请注意，B和C的相对提升计算实际上忽略了变体A。这是因为，例如，当我们计算B的reluplift时，其他变体的最大值很可能是变体C。同样，当我们计算C的reluplift时，它很可能与B进行比较。

• B和C的提升图告诉我们，我们还不能在这两个变体之间明确地选出一个赢家，因为94%的高密度区间（HDI）仍然与0重叠。我们需要更大的样本量来检测23%与22.8%的转化率差异。

• 这种方法的一个缺点是我们不能直接说明这些变体相对于变体A（对照组）的提升是多少。这个数字在实践中通常很重要，因为它允许我们估计如果A/B测试的更改推广到所有访问者时的总体影响。不过，我们可以通过重新定义问题为“A与其他两个变体相比有多糟糕”（这在变体A的相对提升分布中显示）来近似得到这个数字。

值转换

现在，如果我们想要比较A/B测试变体在产生多少收入方面的差异，以及/或者估计一个获胜变体带来的额外收入是多少？我们不能使用Beta-Binomial模型来解决这个问题，因为每个访问者的可能值现在处于范围[0, Inf)内。VWO论文中提出的模型如下：

单个访问者产生的收入是 revenue = probability of paying at all * mean amount spent when paying：

\[\mathrm{Revenue}_i = \mathrm{Bernoulli}(\theta)_i * \mathrm{Exponential}(\lambda)_i I(\mathrm{Bernoulli}(\theta)_i = 1)\]

我们假设在支付时，支付的概率与平均支付金额无关。这是实践中典型的假设，除非我们有理由相信这两个参数有依赖关系。基于此，我们可以为支付的游客总数和购买游客的总支付金额分别创建模型（假设每个游客的行为之间相互独立）：

\[c \sim \sum^N\mathrm{Bernoulli}(\theta) = \mathrm{Binomial}(N, \theta)\]

\[r \sim \sum^K\mathrm{Exponential}(\lambda) = \mathrm{Gamma}(K, \lambda)\]

其中 \(N\) 是访问者的总数，\(K\) 是至少有一次购买的访问者的总数。

我们可以重用之前的Beta-二项式模型来建模伯努利转化率。对于平均购买金额，我们使用Gamma先验（这也是Gamma似然的共轭先验）。所以在两变体测试中，设置如下：

\[\theta_A \sim \theta_B \sim \mathrm{Beta}(\alpha_1, \beta_1)\]

\[\lambda_A \sim \lambda_B \sim \mathrm{Gamma}(\alpha_2, \beta_2)\]

\[c_A \sim \mathrm{Binomial}(N_A, \theta_A), c_B \sim \mathrm{Binomial}(N_B, \theta_B)\]

\[r_A \sim \mathrm{Gamma}(c_A, \lambda_A), r_B \sim \mathrm{Gamma}(c_B, \lambda_B)\]

\[\mu_A = \theta_A * \dfrac{1}{\lambda_A}, \mu_B = \theta_B * \dfrac{1}{\lambda_B}\]

\[\mathrm{reluplift}_B = \mu_B / \mu_A - 1\]

\(\mu\) 这里表示每位访客的平均收入，包括那些没有进行购买的访客。这是捕捉整体收入效应的最佳方式——某些变体可能会提高平均销售价值，但会减少实际支付的访客比例（例如，如果我们在着陆页上推广更昂贵的商品）。

下面我们将模型设置转化为代码，并执行先验预测检查。

@dataclass
class GammaPrior:
    alpha: float
    beta: float

@dataclass
class RevenueData:
    visitors: int
    purchased: int
    total_revenue: float

class RevenueModel:
    def __init__(self, conversion_rate_prior: BetaPrior, mean_purchase_prior: GammaPrior):
        self.conversion_rate_prior = conversion_rate_prior
        self.mean_purchase_prior = mean_purchase_prior

    def create_model(self, data: List[RevenueData], comparison_method: str) -> pm.Model:
        num_variants = len(data)
        visitors = [d.visitors for d in data]
        purchased = [d.purchased for d in data]
        total_revenue = [d.total_revenue for d in data]

        with pm.Model() as model:
            theta = pm.Beta(
                "theta",
                alpha=self.conversion_rate_prior.alpha,
                beta=self.conversion_rate_prior.beta,
                shape=num_variants,
            )
            lam = pm.Gamma(
                "lam",
                alpha=self.mean_purchase_prior.alpha,
                beta=self.mean_purchase_prior.beta,
                shape=num_variants,
            )
            converted = pm.Binomial(
                "converted", n=visitors, p=theta, observed=purchased, shape=num_variants
            )
            revenue = pm.Gamma(
                "revenue", alpha=purchased, beta=lam, observed=total_revenue, shape=num_variants
            )
            revenue_per_visitor = pm.Deterministic("revenue_per_visitor", theta * (1 / lam))
            theta_reluplift = []
            reciprocal_lam_reluplift = []
            reluplift = []
            for i in range(num_variants):
                if comparison_method == "compare_to_control":
                    comparison_theta = theta[0]
                    comparison_lam = 1 / lam[0]
                    comparison_rpv = revenue_per_visitor[0]
                elif comparison_method == "best_of_rest":
                    others_theta = [theta[j] for j in range(num_variants) if j != i]
                    others_lam = [1 / lam[j] for j in range(num_variants) if j != i]
                    others_rpv = [revenue_per_visitor[j] for j in range(num_variants) if j != i]
                    if len(others_rpv) > 1:
                        comparison_theta = pm.math.maximum(*others_theta)
                        comparison_lam = pm.math.maximum(*others_lam)
                        comparison_rpv = pm.math.maximum(*others_rpv)
                    else:
                        comparison_theta = others_theta[0]
                        comparison_lam = others_lam[0]
                        comparison_rpv = others_rpv[0]
                else:
                    raise ValueError(f"comparison method {comparison_method} not recognised.")
                theta_reluplift.append(
                    pm.Deterministic(f"theta_reluplift_{i}", theta[i] / comparison_theta - 1)
                )
                reciprocal_lam_reluplift.append(
                    pm.Deterministic(
                        f"reciprocal_lam_reluplift_{i}", (1 / lam[i]) / comparison_lam - 1
                    )
                )
                reluplift.append(
                    pm.Deterministic(f"reluplift_{i}", revenue_per_visitor[i] / comparison_rpv - 1)
                )
        return model

对于Beta先验，我们可以设置一个与之前类似的先验——以0.5为中心，alpha 和 beta 的大小决定了分布的“薄度”。

我们需要对Gamma先验更加小心。Gamma先验的均值是 \(\dfrac{\alpha_G}{\beta_G}\)，并且需要根据现有的平均购买值设置为一个合理的值。例如，如果 alpha 和 beta 被设置为使得均值为1美元，但网站的每位访客的平均收入要高得多，为100美元，这可能会影响我们的推断。

c_prior = BetaPrior(alpha=5000, beta=5000)
mp_prior = GammaPrior(alpha=9000, beta=900)

data = [
    RevenueData(visitors=1, purchased=1, total_revenue=1),
    RevenueData(visitors=1, purchased=1, total_revenue=1),
]

with RevenueModel(c_prior, mp_prior).create_model(data, "best_of_rest"):
    revenue_prior_predictive = pm.sample_prior_predictive(samples=10000, return_inferencedata=False)

fig, ax = plt.subplots()
az.plot_posterior(revenue_prior_predictive["reluplift_1"], ax=ax, **plotting_defaults)
ax.set_title(f"Revenue Rel Uplift Prior Predictive, {c_prior}, {mp_prior}", fontsize=10)
ax.axvline(x=0, color="red");

与伯努利转化模型类似，先验预测提升分布的宽度将取决于我们先验的强度。有关使用弱先验与强先验的优缺点的讨论，请参见伯努利转化部分。

接下来我们为模型生成合成数据。我们将生成以下场景：

• 相同的购买倾向和相同的平均购买价值。

• 购买倾向较低，平均购买价值较高，但每位访客的总收入相同。

• 更高的购买倾向和更高的平均购买价值，以及每位访客的总体更高收入。

def generate_revenue_data(
    variants: List[str],
    true_conversion_rates: List[float],
    true_mean_purchase: List[float],
    samples_per_variant: int,
) -> pd.DataFrame:
    converted = {}
    mean_purchase = {}
    for variant, p, mp in zip(variants, true_conversion_rates, true_mean_purchase):
        converted[variant] = bernoulli.rvs(p, size=samples_per_variant)
        mean_purchase[variant] = expon.rvs(scale=mp, size=samples_per_variant)
    converted = pd.DataFrame(converted)
    mean_purchase = pd.DataFrame(mean_purchase)
    revenue = converted * mean_purchase
    agg = pd.concat(
        [
            converted.aggregate(["count", "sum"]).rename(
                index={"count": "visitors", "sum": "purchased"}
            ),
            revenue.aggregate(["sum"]).rename(index={"sum": "total_revenue"}),
        ]
    )
    return agg

def run_scenario_value(
    variants: List[str],
    true_conversion_rates: List[float],
    true_mean_purchase: List[float],
    samples_per_variant: int,
    conversion_rate_prior: BetaPrior,
    mean_purchase_prior: GammaPrior,
    comparison_method: str,
) -> az.InferenceData:
    generated = generate_revenue_data(
        variants, true_conversion_rates, true_mean_purchase, samples_per_variant
    )
    data = [RevenueData(**generated[v].to_dict()) for v in variants]
    with RevenueModel(conversion_rate_prior, mean_purchase_prior).create_model(
        data, comparison_method
    ):
        trace = pm.sample(draws=5000, chains=2, cores=1)

    n_plots = len(variants)
    fig, axs = plt.subplots(nrows=n_plots, ncols=1, figsize=(3 * n_plots, 7), sharex=True)
    for i, variant in enumerate(variants):
        if i == 0 and comparison_method == "compare_to_control":
            axs[i].set_yticks([])
        else:
            az.plot_posterior(trace.posterior[f"reluplift_{i}"], ax=axs[i], **plotting_defaults)
        true_rpv = true_conversion_rates[i] * true_mean_purchase[i]
        axs[i].set_title(f"Rel Uplift {variant}, True RPV = {true_rpv:.2f}", fontsize=10)
        axs[i].axvline(x=0, color="red")
    fig.suptitle(f"Method {comparison_method}, {conversion_rate_prior}, {mean_purchase_prior}")

    return trace

场景1 - 相同的潜在购买率和平均购买价值

_ = run_scenario_value(
    variants=["A", "B"],
    true_conversion_rates=[0.1, 0.1],
    true_mean_purchase=[10, 10],
    samples_per_variant=100000,
    conversion_rate_prior=BetaPrior(alpha=5000, beta=5000),
    mean_purchase_prior=GammaPrior(alpha=9000, beta=900),
    comparison_method="best_of_rest",
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [theta, lam]

100.00% [6000/6000 00:14<00:00 Sampling chain 0, 0 divergences]

100.00% [6000/6000 00:13<00:00 Sampling chain 1, 0 divergences]

Sampling 2 chains for 1_000 tune and 5_000 draw iterations (2_000 + 10_000 draws total) took 28 seconds.

• 94% HDI 包含 0，符合预期。

场景2 - 较低的购买率，较高的平均购买额，相同的每位访客总收入

scenario_value_2 = run_scenario_value(
    variants=["A", "B"],
    true_conversion_rates=[0.1, 0.08],
    true_mean_purchase=[10, 12.5],
    samples_per_variant=100000,
    conversion_rate_prior=BetaPrior(alpha=5000, beta=5000),
    mean_purchase_prior=GammaPrior(alpha=9000, beta=900),
    comparison_method="best_of_rest",
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [theta, lam]

100.00% [6000/6000 00:13<00:00 Sampling chain 0, 0 divergences]

100.00% [6000/6000 00:14<00:00 Sampling chain 1, 0 divergences]

Sampling 2 chains for 1_000 tune and 5_000 draw iterations (2_000 + 10_000 draws total) took 28 seconds.

• 每位访客的平均收入（RPV）的94% HDI包含0，如预期所示。

• 在这些情况下，绘制theta（购买任何东西的比率）和1 / lam（平均购买价值）的相对提升分布图也很有用，以了解A/B测试如何影响访问者行为。我们在下面展示了这一点：

axs = az.plot_posterior(
    scenario_value_2,
    var_names=["theta_reluplift_1", "reciprocal_lam_reluplift_1"],
    **plotting_defaults,
)
axs[0].set_title(f"Conversion Rate Uplift B, True Uplift = {(0.04 / 0.05 - 1):.2%}", fontsize=10)
axs[0].axvline(x=0, color="red")
axs[1].set_title(
    f"Revenue per Converting Visitor Uplift B, True Uplift = {(25 / 20 - 1):.2%}", fontsize=10
)
axs[1].axvline(x=0, color="red");

• 变体B的转化率提升的HDI远低于0，而每位转化访客的收入HDI远高于0。因此，模型能够捕捉到购买访客的减少以及平均购买金额的增加。

场景3 - 更高的购买倾向和平均购买价值

_ = run_scenario_value(
    variants=["A", "B"],
    true_conversion_rates=[0.1, 0.11],
    true_mean_purchase=[10, 10.5],
    samples_per_variant=100000,
    conversion_rate_prior=BetaPrior(alpha=5000, beta=5000),
    mean_purchase_prior=GammaPrior(alpha=9000, beta=900),
    comparison_method="best_of_rest",
)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [theta, lam]

100.00% [6000/6000 00:15<00:00 Sampling chain 0, 0 divergences]

100.00% [6000/6000 00:13<00:00 Sampling chain 1, 0 divergences]

Sampling 2 chains for 1_000 tune and 5_000 draw iterations (2_000 + 10_000 draws total) took 30 seconds.

• 如预期，94% HDI 对于变体 B 高于 0。

请注意，在实际使用值转换时（在我们只是模拟合成数据时不会出现这种情况），一个值得关注的问题是异常值的存在。例如，一个变体中的访问者可能会花费数千美元，观察到的收入数据不再遵循像Gamma这样的“良好”分布。在进行统计分析之前，通常需要对这些异常值进行插补（我们必须小心地完全移除它们，因为这可能会导致推断偏差），或者在决策时回退到伯努利转换。

进一步阅读

在实践中实施贝叶斯框架来分析A/B测试时，还有许多其他需要考虑的因素。一些包括：

• 我们如何选择我们的先验分布？

• 在实践中，人们每天都会查看A/B测试结果，而不仅仅是在测试结束时查看一次。我们如何平衡更快地发现真实差异与最小化错误发现（即“早期停止”问题）？

• 如果我们使用贝叶斯模型来分析结果，我们如何使用功效分析来规划A/B测试的长度和大小？

• 除了转化率（每个访客的伯努利随机变量）之外，在线软件中的许多价值分布无法用像正态分布、伽马分布等良好的密度函数来拟合。我们如何对这些进行建模？

各种教科书和在线资源更详细地探讨了这些领域。《Doing Bayesian Data Analysis》[Kruschke, 2014] 由 John Kruschke 撰写，是一个很好的资源，并且已被翻译为 PyMC 这里。

我们还计划在这些主题上创建更多的PyMC教程，敬请期待！

作者

• 由 Cuong Duong 于2021年5月撰写（pymc-examples#164）

• 由 percevalve 于2022年5月重新执行（pymc-examples#351）

参考资料

[1]

克里斯·斯图奇奥。VWO的贝叶斯A/B测试。2015年。URL: https://vwo.com/downloads/VWO_SmartStats_technical_whitepaper.pdf.

[2]

约翰·克鲁施克。贝叶斯数据分析：R、JAGS和Stan的教程。学术出版社，2014年。

水印

%load_ext watermark
%watermark -n -u -v -iv -w -p pytensor,xarray

Last updated: Thu Sep 22 2022

Python implementation: CPython
Python version       : 3.8.10
IPython version      : 8.4.0

pytensor: 2.7.3
xarray: 2022.3.0

matplotlib: 3.5.2
sys       : 3.8.10 (default, Mar 25 2022, 22:18:25) 
[Clang 12.0.5 (clang-1205.0.22.11)]
numpy     : 1.22.1
pymc      : 4.0.1
arviz     : 0.12.1
pandas    : 1.4.3

Watermark: 2.3.1