Содержание страницы
Байесовское А/В-тестирование в коде
Здравствуйте и вновь добро пожаловать на занятия по теме «Байесовское машинное обучение на языке Python, часть 1».
В продолжение темы по байесовскому AB-тестированию: в этой статье мы создадим в коде алгоритм байесовского А/В-тестирования, затем немного оглянемся назад что мы изучили, найдем пороговое значение без p-значения, сделаем демонстрацию схождение сэплирования по Томпсону и в завершение пройдем приближённые доверительные интервалы и бета-апостериор.
Если вы не хотите сами писать код, отыщите файл bayesian_bandit.py в репозитарии и просто запустите файл.
Итак начнём с импорта библиотек – Matplotlib, чтобы можно было наглядно увидеть наш пример, Numpy и функцию для бета-распределения из библиотеки Scipy.
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta
Установим количество испытаний, равное 2000, а также показатели кликабельности для наших бандитов.
NUM_TRIALS = 2000
BANDIT_PROBABILITIES = [0.2, 0.5, 0.75]
Определим класс Bandit, который будет действовать в качестве игрового автомата со своей вероятностью выигрыша и параметрами бета-распределения a и b.
class Bandit(object):
def __init__(self, p):
self.p = p
self.a = 1
self.b = 1
def pull(self):
return np.random.random() < self.p
def sample(self):
return np.random.beta(self.a, self.b)
def update(self, x):
self.a += x
self.b += 1 – x
def plot(bandits, trial):
x = np.linspace(0, 1, 200)
for b in bandits:
y = beta.pdf(x, b.a, b.b)
plt.plot(x, y, label=”real p: %.4f” % b.p)
plt.title(“Bandit distributions after %s trials” % trial)
plt.legend()
plt.show()
Теперь напишем функцию для собственно запуска эксперимента.
def experiment():
bandits = [Bandit(p) for p in BANDIT_PROBABILITIES]
sample_points = [5,10,20,50,100,200,500,1000,1500,1999]
for i in xrange(NUM_TRIALS):
bestb = None
maxsample = -1
for b in bandits:
sample = b.sample()
allsamples.append(“%.4f” % sample)
if sample > maxsample:
maxsample = sample
bestb = b
if i in sample_points:
print “current samples: %s” % allsamples
plot(bandits, i)
x = bestb.pull()
bestb.update(x)
if __name__ == “__main__”:
experiment()
Итак, запустим и посмотрим, что у нас получится.
Итак, после 5 испытаний все дисперсии очень велики.
Далее видим ситуацию после 10 испытаний. «Красный» бандит по-прежнему выше всех. Ещё не совсем ясно, какой у него показатель кликабельности, но он выше других.
Затем видим ситуацию после 20 испытаний. Наша уверенность в показателе кликабельности «красного» крепнет.
После 50 испытаний получаем первую точность в случае «зелёного». Его взвешенная вероятность даже выше, чем у «красного», хотя его настоящий показатель кликабельности равен всего 0,5, но это из-за элемента случайности.
В результате 100 испытаний «зелёный» становится ниже и ближе к своему истинному значению.
Далее идут 200 и 500 испытаний. Заметьте, как ярко выражен «красный», тогда как «зелёный» и «синий» имеют всё ещё довольно «толстые» графики.
1000 испытаний. «Зелёный» и «синий» по-прежнему имеют «толстые» графики, а «красный» становится всё более ярко выраженным.
1500 испытаний. Та же картина.
И последняя контрольная точка – 2000 испытаний. «Зелёный» и «синий» так и продолжают оставаться «толстыми» и располагаются слева.
Интересно отметить, что распределения бандитов с более низкими показателями никогда не становятся точными – они остаются относительно «толстыми» по сравнению с лучшим бандитом, но даже в этом случае основной вес их вероятности меньше, чем у лучшего бандита. Другими словами, мы почти уверены, что их истинный показатель кликабельности ниже и в связи с этим нет необходимости их дальнейшего исследования для уточнения оценки.
Оперативная природа байесовского А/В-тестирования
Давайте сделаем шаг назад и оценим, как далеко мы продвинулись. Мы начали с традиционных А/В-тестов, в которых необходимо было делать всевозможные приближения, сложно решить, сколько примеров собрать для дополнительных приближений, и запускать эксперимент до полного его завершения без какой-либо адаптации к уже полученным результатам. Если же остановиться чересчур рано, то, согласно частотной статистике, не будет достигнут корректный результат.
Теперь же у нас есть способ создать адаптивную систему, которая, если подумать, на самом деле даже не требует каких-либо проверок – мы просто запускаем систему и позволяем ей работать. В конце концов она сходится, и мы получаем наилучший возможный результат с учётом уже собранных данных. Нет никакого порогового значения, при котором результаты становятся корректными – мы просто естественным образом сходимся в точке равновесия.
Любопытный аспект байесовского метода выборки заключается в том, что измерения становятся всё точнее с каждым рассмотренным примером. Как так получается? Если внимательнее присмотреться к коду, то можно заметить, что мы не обновляем a и b после сбора всех данных, вместо этого мы обновляем их после каждой игры. Другими словами, апостериорная вероятность превращается в априорную, когда мы добавляем новые данные.
Сравните это с традиционными моделями машинного обучения, когда мы уточняем одновременно все наши учебные данные. Тогда если мы захотим добавить новые учебные данные, то процесс обучения придётся начинать заново, с нуля. В случае же байесовских методов нам нужно лишь добавить числа к a и b, а бета-распределение, поскольку оно зависит от значений a и b, подстроится автоматически.
Это и есть пример оперативного обучения – модель обучается после каждого нового примера. При этом не возникает никаких проблем с оперативной памятью, поскольку a и b занимают постоянное пространство, и по той же причине нет проблем с быстродействием. В некотором роде это похоже на процесс обучения людей и животных – мы берём некоторый недавний опыт и стараемся сформировать некоторое мнение о нём. Мы не пытаемся учесть весь свой полученный ранее опыт и найти некий глобальный оптимум. На самом деле мы и не можем этого сделать, ведь наша «система памяти» несовершенна.
Таким образом, получается множество замечательных результатов, следующих из использования байесовского метода.
Нахождение порогового значения без p-значения
Вернёмся к нашему обсуждению p-значения и зададимся вопросом, можно ли найти пороговое значение для байесовского А/В-тестирования.
Даже несмотря на то, что байесовская выборка очень эффективна, она всё равно использует ресурсы. Надо полагать, нам всё равно необходимо хранить значения a и b в какой-то базе данных, генерировать случайные числа для всех наших бандитов и так далее. А потому если мы уверены, что одна реклама лучше другой, то почему бы её не использовать постоянно? В этой лекции мы обсудим специфическую меру вероятности, конкретнее говоря, вероятность, что μ1 больше μ2. Не забывайте, что μ – это показатель кликабельность и что μ1 и μ2 имеют бета-распределение.
Прежде, чем перейти к вычислению этой вероятности, давайте ещё раз взглянем на определение p-значения. Вспомните, какое оно путанное, сколько споров идёт по её поводу и как легко ошибиться в формулировках, приводя статистиков в бешенство.
Не забывайте, что вы не можете принять нулевую гипотезу, вы можете только не отвергать её. С такими определениями легко понять, почему огромное множество людей любят машинное обучение, а статистика никогда не привлекала к себе такого внимания, хотя они в значительной степени похожи. Как я уже указывал ранее, людям нетехнических специальностей очень легко неправильно понять смысл p-значения и как легко его утерять при донесении информации. Очень распространённое заблуждение состоит в том, что если, скажем, p-значение равно 1%, то это значит, что вероятность того, что μ1 ≠ μ2 составляет 99%. Используя уже указывавшееся определение, это, конечно же, не так, но многие всё равно хотят использовать подобную вероятность – она проста, легка в интерпретации и понятна любому. К счастью, в связи с тем, что μ1 и μ2 имеют распределения, байесовские методы позволяют найти такую вероятность. Вот только как это сделать?
Начнём с основной предпосылки. Мы хотим узнать вероятность того, что μ2 > μ1, что эквивалентно вероятности того, что μ2 – μ1 > 0. Если мы примем, что X = μ2 – μ1, то получим, что P(X) > 0, а это представляет собой площадь под функцией плотности распределения вероятностей по X от 0 до бесконечности. Но мы не знаем распределения X. Вместо этого мы выберём другой путь: мы рассмотрим объединённую плотность распределения вероятностей μ1 и μ2 и найдём площадь под кривой, где μ2 > μ1. Ну, а объединённая плотность распределения вероятностей просто равна произведению двух бета-распределений, поскольку они независимы. На самом деле можно даже создать график объединённой плотности распределения вероятностей в виде тепловой карты или изображения, а нам нужна область в виде верхнего треугольника (см. слайд).
Доказательство выходит за рамки данного курса, но конечный результат в аналитическом виде можно представить в виде следующих уравнений, которые можно легко решить на компьютере:
Однако обратите внимание, что суммирование имеет сложность O(a2), а потому большая скорость вычислений тут не обязательна.
Другой способ получить пороговое значение – определить функцию потерь и остановиться, когда ожидаемое значение потерь падает ниже некоторого порогового значения. Полагая, что μ2 > μ1, мы можем определить функцию потерь:
Затем мы можем вычислить ожидаемое значение функции потерь под объединённым распределением и остановиться, когда это ожидаемое значение будет меньше некоторого порогового:
Таким образом, мы можем вычислить ожидаемое значение функции потерь под объединённым распределением, используя данные интегралы:
Хотя эти интегралы кажутся сложными, их можно найти, использовав предыдущие, куда более сложные (я рекомендую в качестве упражнения попробовать так и сделать):
Таким образом, возможно определить пороговые значения для байесовского А/В-тестирования, хотя в следующей лекции мы увидим, как использование байесовского А/В-тестирования позволяет автоматически дойти до максимума показателя кликабельности без остановки тестирования.
И имей в виду, что использование байесовских методов по-прежнему не мешает вам использовать традиционное А/В-тестирование. Все данные есть, и ничто не мешает использовать их в обычном А/В-тестировании, нахождении p-значения и так далее.
Схождение сэплирования по Томпсону. Демонстрация
Сейчас устроим демонстрацию того, как показатель кликабельности сходится к лучшему Bandit, когда мы используем байесовский метод.
Итак, начнём с импорта библиотек Matplotlib, Numpy, а из предыдущего нашего примера используем класс Bandit.
import matplotlib.pyplot as plt
import numpy as np
from bayesian_bandit import Bandit
Напишем функцию для запуска эксперимента с тремя вероятностями, так что у нас будет три бандита, и количеством попыток.
def run_experiment(p1, p2, p3, N):
bandits = [Bandit(p1), Bandit(p2), Bandit(p3)]
data = np.empty(N)
for i in xrange(N):
j = np.argmax([b.sample() for b in bandits])
x = bandits[j].pull()
bandits[j].update(x)
Найдём ещё кумулятивное среднее значение. В Numpy есть специальная функция для кумулятивной суммы.
data[i] = x
cumulative_average_ctr = np.cumsum(data) / (np.arange(N) + 1)
И выведем графики на экран.
plt.plot(cumulative_average_ctr)
plt.plot(np.ones(N)*p1)
plt.plot(np.ones(N)*p2)
plt.plot(np.ones(N)*p3)
plt.ylim((0,1))
plt.xscale(‘log’)
plt.show()
Попробуем 100 000 испытаний, но вы можете выставить свои значения.
run_experiment(0.2, 0.25, 0.3, 100000)
Запустим и посмотрим, что у нас получится.
Как видим, в долгосрочной перспективе – примерно после 10000 примеров – мы сошлись на уровне показателя кликабельности в 0,3 – это наш лучший бандит.
Приближённые доверительные интервалы и бета-апостериор
И завершение нашей темы мы выясним, насколько точен доверительный интервал, полученный нами ранее, в частности, мы наглядно отобразим приблизительно гауссово распределение, полученное из центральной предельной теоремы, и сравним его с байесовским бета-апостериором. Тут всё очень похоже на код, написанный нами ранее, поэтому не стесняйтесь попытаться написать его самостоятельно. Если же вы не хотите писать код сами, то соответствующий файл называется ci_comparison.py.
Итак, начнём с импорта библиотек – Matplotlib, Numpy, а из Scipy нам понадобятся функции beta и norm.
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import beta, norm
У нас будет 501 испытание с истинным показателем кликабельности 0,5 и значениями a и b, равными 1, – это равномерное распределение. Переменная plot_indices показывает, когда выводить на график на экран.
T = 501
true_ctr = 0.5
a, b = 1, 1
plot_indices = (10, 20, 30, 50, 100, 200, 500)
data = np.empty(T)
for i in xrange(T):
x = 1 if np.random.random() < true_ctr else 0
data[i] = x
Обновляем данные для распределения.
a += x
b += 1 – x
if i in plot_indices:
p = data[:i].mean()
n = i + 1
std = np.sqrt(p*(1-p)/n)
Далее – отображение гауссианы.
x = np.linspace(0, 1, 200)
g = norm.pdf(x, loc=p, scale=std)
plt.plot(x, g, label=’Gaussian Approximation’)
Теперь то же самое для бета-распределения.
posterior = beta.pdf(x, a=a, b=b)
plt.plot(x, posterior, label=’Beta Posterior’)
plt.legend()
plt.title(“N = %s” % n)
plt.show()
Собственно, это и весь код. Запустим и посмотрим, что получится.
Итак, при N = 11 приближённое гауссово распределение и апостериорная вероятность бета-распределения весьма отличаются. При N = 21 они начинают сходиться, и то же самое при N = 31. При N = 51 они становятся ещё ближе, а при N = 101 они почти одинаковы. При N = 201 – схожесть ещё большая. И последнее, N = 501 – практически полная эквивалентность. Между прочим, это ещё и демонстрация истинности центральной предельной теоремы. Можно также увеличить график, чтобы посмотреть, насколько близки наши величины. Как видим, очень, очень близки.
