Содержание страницы
Байесовское А/В-тестирование в коде
Здравствуйте и вновь добро пожаловать на занятия по теме «Байесовское машинное обучение на языке Python, часть 1».
В продолжение темы по байесовскому AB-тестированию: в этой статье мы создадим в коде алгоритм байесовского А/В-тестирования, затем немного оглянемся назад что мы изучили, найдем пороговое значение без p-значения, сделаем демонстрацию схождение сэплирования по Томпсону и в завершение пройдем приближённые доверительные интервалы и бета-апостериор.
Если вы не хотите сами писать код, отыщите файл bayesian_bandit.py в репозитарии и просто запустите файл.
Итакначнём с импорта библиотек – Matplotlib, чтобы можнобыло наглядно увидеть наш пример, Numpy и функцию длябета-распределения из библиотеки Scipy.
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats importbeta
Установимколичество испытаний, равное 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):
returnnp.random.random() < self.p
def sample(self):
returnnp.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 inbandits:
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 inbandits:
sample =b.sample()
allsamples.append(“%.4f” % sample)
if sample> maxsample:
maxsample = sample
bestb =b
if i insample_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 importBandit
Напишемфункцию для запуска эксперимента с тремя вероятностями, так что у нас будет трибандита, и количеством попыток.
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 ifnp.random.random() < true_ctr else 0
data[i] = x
Обновляемданные для распределения.
a += x
b += 1 – x
if i inplot_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 –практически полная эквивалентность. Между прочим, это ещё и демонстрацияистинности центральной предельной теоремы. Можно также увеличить график, чтобыпосмотреть, насколько близки наши величины. Как видим, очень, очень близки.
