Введение

В статье Ретеншен — основная метрика F2P игры, вероятностный подход я рассматривал ретеншен как биномиально распределенную случайную величину и показал, что точность оценки его отдельных измерений невысока.

Однако, можно получить более точную оценку, если использовать данные за несколько дней наблюдений и применить их, используя известные закономерности поведения метрики. В частности, можно отметить следующие доступные средства:

1. Если рассматривать все измеренные данные как временной ряд с некоторыми известными свойствами, можно построить аппроксимирующую кривую, которая будет являться точечной оценкой матожидания метрики.

2. Часто можно показать, что мы имеем дело не с биномальной случайной величиной, а с суммой нескольких случайных величин, что сужает ее дисперсию.

3. Как правило, известны не только данные о возвратах пользователей в один конкретный день после инсталла, но другие данных об их поведении, используя которые, можно гипотетически повысить точность оценки. Например, можно использовать измеренные данные о ретеншене глубины i, для уточнения оценки ретеншена глубины k.

4. Если мы располагаем большим объемом данных о значениях метрики от большого количества проектов, это, потенциально, позволяет более точно подобрать гиперпараметры аппроксимирующей модели и также повысить точность оценки метрики.

В данной статье рассматривается первый пункт из этого списка.

Постановка задачи

Рассмотрим классический ретеншн N-го (например, 1-го дня) дня как случайный временной ряд, значения которого мы определяем экспериментально. На графике этот процесс выглядит примерно так.

Задача - оценить параметры этого случайного ряда, а именно:

1. Получить точечную оценку математического ожидания в каждой точке ряда.

2. Получить интервальную оценку математического ожидания в каждой точке.

Задача по определению является задачей анализа временных рядов. Существуют специализированные фреймворки для решения задач подобного класса, например Prophet, Neural Prophet, Nixtla/statsforecast, о которых также можно найти информацию на хабре: 

Но сегодня я не буду использовать подобный фреймворк, а решаю задачу в общем виде, как задачу нелинейной аппроксимации процесса с некоторыми известными априори свойствами, подробно обосновывая принципы создания аппроксимирующей функции.

К слову, попытка использовать готовый фреймворк как черный ящик для данной задачи, не вникая в суть его работы и не вникая в суть аппроксимируемого процесса, может привести к провалу из-за некоторых особенностей задачи, в частности:

1. Процесс не стационарный. Метрика убывает со временем, на что накладывается выпуск патчей, дисперсия тоже меняется вследствие изменения количества установок приложения.

2. Исходных данных мало. Например, если проект находится в релизе целый год, у нас есть не более 365 точек наблюдения. Зачастую же надо измерить показатели патча, выпущенного только что или даже проекта только что выпущенного в релиз.

3. Данные очень сильно зашумлены.  Например тут автор рассматривает вопрос анализа временного ряда метрики CCU при количестве активных пользователей ~20к-60к в час. Я же говорю о необходимости исследования метрики проекта с тысячей инсталлов в день.

4. Метрика меняется не только вследствие изменения качества игры (что мы хотим измерить), но и-за изменения трафика.

В общем, сегодня и здесь будет деконструкция метрики и последующая нелинейная аппроксимация, а Prophet и его друзья будут потом, в одной из следующих публикаций.

Рисунок 2 иллюстрирует ожидаемый результат аппроксимации.

Обозначения на графике:

• Исходные измеренные данные изображены тонкой линией синего цвета

• Исходная интервальная оценка каждой точки наблюдения изображена плавающим коридором темно-синего цвета

• Аппроксимирующая кривая изображена как сплошная линия красного цвета.

• Новая интервальная оценка, полученная при помощи аппроксимации - коридор серого цвета.

Примечание. Здесь и далее я буду показывать не настоящие данные, а смоделированные. К сожалению, у меня нет возможности публиковать данные с реальных проектов.

Анализ задачи

В задаче анализа временных радов часто используется так называемая аддитивная модель, которая представляет аппроксимируемый случайный процесс в виде суммы:

r_i = g(t_i) + s(t_i) + a(t_i) + e(t_i) (1)

где

r_i , t_i - исходный временной ряд случайной функции

g(t) - линия общего тренда, одномерная функция

s(t) - сезонные явления, одномерная функция

h(t) - аномальные выбросы, случайная функция

e(t) - ошибка, нормально-распределенная случайная величина с нулевым матожиданием

В некоторых случаях корректнее соединять линию тренда и сезонные явления не аддитивно, а мультипликативно.

r_i = g(t_i) * s(t_i) + a(t_i) + e(t_i) (2)

Обе эти математические модели являются традиционными в задаче анализа временных рядов и я буду их использовать.

Покажем, что подобные модели могут быть применены к рассматриваемой здесь метрике "ретеншен".

В конце первой статьи я привел список некоторых хорошо известных закономерностей поведения ретеншна. Большинство из них (но не все!) действительно можно разложить по четырем перечисленным выше составляющим. Рассмотрим это подробнее.

Линия тренда g(t)

Линия тренда - это и есть то, что мы хотим измерить. Это и есть среднее значение метрики, отражающее качество проекта.

Из практики известно, что ретеншен всегда плавно снижается снижается от момента релиза, но также может резко изменить значение в момент выпуска патча. Таким образом линию тренда удобно представить в виде двух слагаемых

g(t) = d(t) + p(t) (3)

где

d(t) - линия общего тренда, показывающая общую тенденцию поведения метрики - постоянное монотонное снижение.

p(t) - влияние патчей. Эта функция - кусочная, потому что она состоит из нескольких участков, каждый из которых соответствует одному конкретному патчу, и патчи не пересекаются.

Функции соединены аддитивно, поскольку это удобнее в работе, и, кроме того, в пределах действия одного патча значения функции d(t) существенно не меняются, следовательно мы сможем подобрать такой вид p(x), чтобы применять его аддитивно.

На рисунке 3 показана линия общего тренда (зеленая линия), линия патчей (черные линии) и их сумма - красная линия.

Ошибка e(t)

В первой статье показано, что ретеншен можно считать биномиально-распределенной случайной величиной, а также, одновременно, нормально-распределенной величиной, потому что при большом количестве установок, согласно центральной предельной теореме, биномиальное распределение становится нормальным. По правде сказать, утверждение о биномиальности взято несколько с потолка, скорее мы имеем дело с суммой нескольких распределений, вероятно, биномиальных. Но, в любом случае, количество установок большое, сотни как минимум, поэтому центральная предельная теорема работает и мы получим распределение, близкое к нормальному.

Ошибка e(t) является случайной функцией от времени, поскольку дисперсия ошибки, в общем случае, не является постоянной, а должна быть вычислена как дисперсия биномиального распределения с известными матожиданием и количеством экспериментов Бернулли (то есть установок приложения). То есть

e(t) = N(0,σ(t)) (4)

где

(5)

σ(t) - дисперсия биномиального распределения ретенешна с матожиданием M(r(t)) и количеством установок (количеством экспериментов Бернулли) I(t)

Аномальные выбросы a(t)

Помимо нормально распределенной ошибки формулах (1) и (2) присутствует еще одна случайная величина - a(t), аномальные выбросы.

Дело в том, что существует ряд аномальных факторов, существенно влияющих на метрику и не подчиняющихся закону биномиального распределения, основанному на количестве установок. Например, технические сбои, выпуск патчей - они приводят к существенным выбросам метрики, это уже не просто нормальный шум! И это надо как-то учитывать!

Оценка метрики в момент действия аномалии - сложная задача, выходящая за рамки настоящей статьи. Но и игнорировать факт существования аномалий нельзя, ведь они происходят действительно очень часто и сильно искажают измерения.

Буду действовать максимально просто:

1. Данные с большими аномалиями - выбрасывать из выборки.

2. Мелкие аномалии, незаметные на фоне шума, можно попробовать игнорировать - делать вид, что ничего не было.

В литературе по машинному обучению пишут, что задача обнаружения аномалий довольно сложна. К счастью, у нас не тот случай! Обнаружить аномалии можно, и не очень и сложно. Да аномалии (как правило) непредсказуемы - мы ведь не знаем (как правило), когда произойдет сбой и как количественно он повлияет на метрики. Но мы, к счастью, анализируем метрику постфактум и обладаем данными не только об этой метрике, но и огромным объемом других данных. Есть, например, метрика CCU, показывающая количество пользователей в игре в текущий момент с точностью до секунд, есть системные логи с ошибками - это отличные индикаторы аномалий. Современные средства машинного обучения отлично справятся с задачей их обнаружения. Кстати, аномальными днями зачастую следует считать и дни выпуска патчей, ведь патч тоже имеет признаки технического сбоя, о чем я писал в статье 1.

Сезонность s(t)

Функция s(t) показывает зависимость метрики от сезона. Сезонные изменения могут отражать недельные, месячные и годовые периодические процессы.

Недельная сезонность наиболее ярко выражена и хорошо объясняется поведением пользователей в разные дни недели, о чем я писал в первой статье.

Месячная сезонность как правило, очень слабо выражена, хотя иногда, в некоторых регионах, заметна. Например, в странах, где законодательно жестко зафиксирован день выплаты зарплаты, существует явная сезонность метрик, связанных с продажами, что влияет и на игровые метрики.

Годовая сезонность выражается факторами двух типов: краткосрочные события, соответствующие коротким праздникам (причем, в каждой стране - свои праздники) и долгосрочные сезонные изменения, связанные со сменой времен года. Годовая сезонность также может изменять паттерн недельной и месячной сезонности. Например, во время длинных каникул, все дни ведут себя как выходные дни.

Я не буду рассматривать здесь годовую и месячную сезонность, ограничившись недельной, как наиболее ярко выраженной.

Изменение метрики из-за изменения трафика

В формулах (1) и (2) нет очень важной составляющей  - метрика меняется из-за изменения качественного состава пользователей, то есть из-за изменения траффика.

Некорректно включать это в состав функции тренда g(x), поскольку мы измеряем ретеншен как показатель качества продукта, а не влияние состава пользователей на эту метрику. Считать изменения трафика аномалией и относить их в слагаемое a(t) тоже некорректно -  это не какое-то разовое событие, нельзя просто отбросить данные за 1-2 дня и успокоиться. 

Хороший и простой способ исключить фактор траффика - производить все измерения на некотором эталонном трафике - из эталонного источника, эталонного региона, эталонного объема. Это позволяет, с определенной натяжкой, утверждать, что мы будем измерять качество продукта на некоторых унифицированных пользователях, следовательно состав этих пользователей уже не влияет на метрику, и, следовательно, мы измеряем именно качество продукта. Такой подход требует определенной дисциплины в маркетинге и уменьшает размер ежедневной выборки, но это едва ли не единственное, что действительно можно сделать. Именно так мы и поступим.

Если вы плохо знаете о траффике и маркетинге F2P игры, и не понимаете, в чем проблема, под катом я поместил чуть более подробную информацию о ней.

Чуть подробнее о проблеме трафика

Предположим, у вас был один источник трафика, скажем, Unity Ads, но потом почему-то вы вдруг решили, что он плохой (например, дорогой какой-то, и плохо масштабируемый), и от него избавились. Совсем. И вместо него включили Google Admob, например.

И о чудо, ретеншен изменился! Был 40%, например, а стал 35%.

Как же так? Может быть вы сломали игру? Ее качество упало? Да нет же, в игре точно ничего не меняли, вы не выпустили патч, не стартовали игровой ивент, ничего не изменилось! Это просто другой трафик, такова жизнь. Ретеншн является, строго говоря, не метрикой качества игры, а метрикой совместимости игры с некоторыми пользователями. Меняется источник трафика, значит меняются пользователи, значит меняется и метрика.

И проблема в том, что это происходит вообще постоянно. Пример выше с резким переходом с одного источника данных на другой немного утрирован (хотя и такое бывает), но мелкие изменения в маркетинге происходят вообще всегда и постепенно они могут накапливаться в очень крупные - вы действительно можете поменять полностью всю маркетинговую стратегию проекта за несколько месяцев.

И как быть? Как с этим жить?

Одна из неплохих рекомендаций - надо всегда иметь некоторый эталонный источник траффика из эталонного региона для контроля эталонного значения метрик. Вот, предположим, у вас игра ориентирована на Россию, у вас здесь много пользователей, это важный для вас регион. ОК, вот тогда пусть какая-то довольно существенная доля этих российских пользователей будет всегда из одного и того же источника с одной и той же рекламной кампанией. Всегда! А если вы хотите поменять эталонный источник, надо некоторое время, довольно продолжительное, держать включенными оба эти источника - старый и новый, чтобы подстроить метрики одного источника под другой. Проблема этого подхода (как и любой другой сегментации) - снижение статистической значимости данных. Если у вас всего 1000 пользователей в день, и из них 100 - из эталонного источника трафика, следовательно, доверительный интервал измерений (смотрим статью 1) расширяется примерно в 3.5 раза! Не так и просто, увы.

Кстати, часто эталонным трафиком считают органический трафик, что далеко не всегда корректно, ведь органика ведет себя очень непредсказуемо. Например, фичеринги попадают в органику, но это другая органика, это не ваша основная аудитория, причем состав этой органики зависит от типа фичеринга, а это полный рандом. Сегодня будет один фичеринг - трафик будет одной структуры, через неделю другой - трафик другой структуры.

Глобальные факторы несезонного характера

Некоторые факторы, очень сильно влияющие на метрику, не носят какой-то сезонный характер. Например, ковидный локдаун. Каждый такой фактор надо рассматривать индивидуально. 

В настоящей статье эти факторы не рассматриваются.

Метод наименьших квадратов

Традиционно для оценки математического ожидания случайного процесса используют аппроксимацию по методу наименьших квадратов (или МНК). Я тоже буду использовать МНК, но прежде хотелось бы обосновать возможность применения этого метода, потому что он вообще-то может быть использован не всегда.

МНК используется и дает максимально точный результат в том случае, если:

1. Аппроксимируемый случайный процесс является серией нормально-распределенных испытаний.

2. Изматывания независимы.

3. Дисперсия испытаний постоянна.

4. Кроме того также делается предположение, что исследователь знает аналитически точно функцию математического ожидания этой случайной функции как обычную одномерную функцию с некоторыми коэффициентами, однако не знает значения некоторых коэффициентов, входящих в эту функцию. Задача аппроксимации (также именуемой, кстати задачей регрессии) - определить эти коэффициенты.

Рассмотрим каждое исходное требование. 

Условие 1. Нормальная форма распределения

Выше мы уже договорились исключить из выборки аномальные выбросы и объяснили, что функция ошибки, присутствующая в формулах (1) и (2) - нормально распределенная величина. Условие выполнено.

Уловия 2. Независимость испытаний

В литературе по регрессии условие независимости испытаний часто называют самым натянутым и спорным. Однако, на практике все вынуждены допускать независимость испытаний, чтобы не переусложнять задачу. Что касается именно ретеншена, то требование независимости интуитивно кажется довольно правильным. Метрика складывается из поведения совершенно различных людей, слабо связанных между собой. Теоретически можно предположить, что в течении некоторого промежутка времени в игру приходят тысячи людей, которые активно общаются между собой. Первые пришедшие люди из этой группы передают свое впечатление об игре другим и таким образом влияют на их поведение. Но такая организованная группа совершенно маловероятна.

Таким образом, да, будем считать условие выполненным.

Условие 3. Дисперсия во всех точках постоянна

В общем виде требования нарушено, ведь количество инсталлов ежедневно меняется, и, следовательно, согласно формуле дисперсии биномиальной случайной величины, дисперсия тоже будет постоянно меняться. К счастью, это легко решается введением весовых коэффициентов в метод наименьших квадратов.

Итак, это условие немного нарушено, но существует простое решение.

Условие 4. Формула математического ожидания известна.

Требование звучит довольно жестко. Необходимо заранее знать формулу функции, которую мы хотим найти! Как же так? Мы хотим заранее знать то, что хотим найти.

Ну да, выбор аппроксимирующей функции - этот процесс не простой, и, отчасти искусство. Нельзя взять одну аппроксимирующую функцию, скажем линейную, и применять ее во всех случаях жизни. Здесь необходимо использовать техническую эрудицию, обобщить все априорные данные об исследуемом процессе, а также учитывать объем обучающей выборки, чтобы аппроксимируемая функция не оказалась слишком сложной.

Таким образом, необходимо выбрать хороший вид аппроксимирующей функции, которая будет качественно отражать природу рассматриваемого процесса. Конечно, нельзя гарантировать, что функция выбрана идеально точно, и это неверный выбор будет вносить некоторую погрешность в измерения.

Методика аппроксимации

Итак, подытожу все написанное выше и сформулирую методику аппроксимации.

Исходные данные задачи

• фактические измеренные значения ретеншена - временной ряд {r_i}, {t_i}, i = 1...n

• фактическое количество установок приложения - временной ряд {I_i}, {t_i}, i = 1...n

Значения измерены в течении максимально доступного нам количества дней после релиза на некотором эталонном трафике, не меняющемся в течении всего периода наблюдений

Математическая модель

Будем рассматривать r_i как функцию вида:

r_i = d(t_i) + p(t_i) + s(t_i) + a(t_i) + N(0,σ(t_i)) (6)

или

r(t_i) = ( d(t_i) + p(t_i) ) * s(t_i) + a(t_i) + N(0,σ(t_i)) (7)

где:

d(t_i) - линия общего тренда,

p(t_i) - влияние патчей,

s(t_i) - недельная сезонность,

a(t_i) - аномальные выбросы,

N(0,σ(t_i)) - нормально-распределенная ошибка с нулевым матожиданием и дисперсией σ(t)

- дисперсия биномиального распределения ретенешна

с матожиданием M(r(t_i)) и количеством установок i(t_i)

Задачи

1. Определить функцию общей линии тренда d(t_i), как при условии известной функции сезонности s(t), так и при неизвестной.

2. Определить линию тренда с учетом воздействия патчей, то есть d(t) + p(t), как при условии известной сезонности s(t), так и при неизвестной.

3. Определить функцию недельной сезонности s(t)

4. Определить общий вид аппроксимирующей функции, тое есть функцию:

   R(t) = d(t) + p(t) + s(t) (аддитивная модель) (8)

   R(t) = ( d(t) + p(t) ) * s(t) (мультипликативная модель) (9)

Методика аппроксимации

1. Выбрасываем из исходной выборки аномальные дни, для обнаружения которых используем какие-то другие, дополнительные, данные.  

2. Создаем аппроксимирующую функцию, состоящую из трех составляющих

   d(t) - линия общего тренда,

   p(t) - влияние патчей,

   s(t) - недельная сезонность,

3. Для построения аппроксимации применяем метод взвешенных наименьших квадратов с учетом меняющейся дисперсии σ(t)

4. Качество работы модели будем оценивать грубо аналитически, а также статистически на основе предыдущего опыта ее применения, поскольку наша обучающая выборка слишком мала, чтобы выделять из нее тестовую и валидационную.

Предложенная методика реализована в виде Jupyter-ноутбука retention-rate-approximator, выложенного на github в открытый доступ

Аппроксимирующая функция

Ключевой момент рассматриваемого метода аппроксимации - это, конечно, вид аппроксимирующей функции.

Аппроксимирующая функция состоит из трех функций

1. Линия общего тренда d(t)

2. Линии патчей p(t)

3. Сезонность s(t)

Как я уже писал выше, линия тренда и линия патчей объединяются аддитивно, а сезонность может быть применена аддитивно или мультипликативно. Мне кажется более корректно применять сезонность мультипликативно, но, на всякий случай, в retention-rate-approximator я реализовал оба варианта, поскольку это было не сложно.

Линия общего тренда - убывающая функция

Из практики известно, что в момент релиза игры метрика максимальна, далее снижается с убывающей скоростью, но не доходит до нуля никогда. Из этих простых соображений сформулируем требования к линии общего тренда d(t).

1. Функция определена на значениях аргумента >= 0.

2. Функция непрерывная, монотонно убывающая, первая производная всегда отрицательна.

3. На бесконечных значениях аргумента функция медленно приближается к ассимптоте.

4. Функция вогнутая, то есть вторая производная всегда положительна.

Множество функций удовлетворяет этим требованиям. Например, можно использовать варианты предложенные в главе "Как моделировать удержание" книги Василия Сабирова. "Игра в цифры" (правда, автор их применяет для аппроксимации несколько другого процесса - он рассматривает ретеншен как функцию от количества дней, прошедших с момента установки, а не как функцию времени от момента релиза проекта). Можно попробовать использовать экспоненту, гиперболический тангенс и сигмоиду - все эти функции формально удовлетворяют заявленным требованиям в диапазонах значений x > 0 и правильно выбранных коэффициентах. Также при небольшом количестве наблюдений можно использовать линейную функцию и просто константу.

Список некоторых возможных кандидатов приведен в таблице 1

Недостатком всех этих функций (за исключением первых двух) является их нелинейность как по отношению к аргументу так и весам. Следовательно, регрессия будет нелинейная, свойства ее документированы плохо, сходимость может быть плохой.

Есть способ перейти от нелинейной регрессии к линейной, если вместо одной простой функции использовать ряд из нескольких базисных функций без весов в знаменателе. Однако, в этом случае количество весов будет большое, и, следовательно, на нашем небольшом объеме данных мы не сможем правильно обучить модель. Целесообразнее использовать одну нелинейную функцию но с малым количеством весов.

В Jupyter-ноутбуке retention-rate-approximator я реализовал несколько функций из предложенного выше набора, и вы можете их попробовать на практике. Также можете реализовать свою функцию по аналогии с уже реализованными, создав класс из 3-х простых методов, это совсем нетрудно.

Ниже по тексту я напишу немного о своих собственных экспериментах с аппроксимирующей функцией и сделанных мною выводах.

Линия патча

Мы хотим проанализировать результаты внедрения патча как можно скорее, поэтому аппроксимирующая функция должна быть максимально проста.

В большинстве случаев единственно возможным вариантом аппроксимации патча будет просто константа - у каждого патча один весовой коэффициент. То есть, патч добавляет просто постоянное смещение к линии тренда. p(t) = w₀

Если данных много и они слабо зашумлены вместо константы можно использовать линейную функцию. Два весовых коэффициента на каждый патч. p(t) = w₀ + w_{1 *} t

В ряде случаев, наверное, можно использовать и функцию с 3-мя и 4-мя коэффициентами из числа рассмотренных выше. Но для этого необходимо гигантское количество инсталлов.

Аппроксимация недельной сезонности

Для аппроксимации недельной сезонности часто использую для каждого дня недели свой весовой коэффициент, задающий смещение этого дня недели. Мне это подход кажется разумным, так и поступлю без дополнительных обоснований.

Нетрудно заметить, что, если каждому дню недели будет соответствовать свой вес, то есть отдельное отдельное слагаемое в аппроксимирующей функции, и мы применяем функцию недельной сезонности аддитивно (формула (6)), то один из этих 7-ми весовых коэффициентов будет являться линейной комбинации 6-ти других значений. Метод оптимизации никогда никогда не сможет найти оптимальное значение весов такой модели, ведь если увеличивать смещение w₀ линии тренда, и одновременно, на такую же величину уменьшать все недельные коэффициенты, результирующая функция не меняется.

Если функция сезонности применяется мультипликативно (формула (7)), то указанная проблема тоже присутствует, но она будет применена не аддитивно, а мультипликативно.

Проблема решается двумя способами.

1. Выбрать некоторый эталонный день недели, в который вес модели должен будет не настраиваемым, а равен нейтральному значению (ноль при аддитивном объединении и единица - при мультипликативном). Остальные другие веса сезонной функции будут задавать отклонения от этого эталонного дня.

2. Потребовать, чтобы среднее значение весов недельной модели равнялось нейтральному значению (ноль при аддитивном объединении и единица - при мультипликативном).

Мне больше нравится второй вариант, поскольку мы не знаем априори, какой именно день недели следует считать эталонным.

Для того, чтобы реализовать это требование при обучении модели я использую функцию регуляризации. Возможно, это несколько нестандартный случай применения регуляризатора (обычно используются функции L1 и L2), но ведь я хочу установить динамическое ограничение на веса модели, а именно это и делают регуляризаторы.

Чуть более подробно объясняю, причем здесь регуляризаторы.

О функциях регуляризации написано немало, но, в большинстве источников рассматривается только два наиболее распространенных ее вида - регуляризаторы L1 и L2. У меня несколько другая ситуация, поэтому позволю себе небольшое лирическое отступление.

Для построения регрессии необходимо две вещи - аппроксимирующая функция, и обучающие данные. Вид аппроксимирующей функции задает априорные теоретические данные об исследуемом процессе. Например, мы знаем, что процесс периодической, поэтому и аппроксимирующая сезонная функция будет периодическая. Обучающие данные - это практические данные. В процессе обучения модели теория соединяется с практикой.

Но иногда у исследователя еще некоторая дополнительная теоретическая информация, которую сложно выразить в виде аппроксимирующей функции, но зато можно выразить в виде некоторого динамического ограничения на веса этой функции. Например в, нашем случае, мы должны потребовать, чтобы в среднем все веса сезонной периодической модели равнялись нулю.

Для этого к обычной функции потерь (вычисленной, в нашем случае по методу наименьших квадратов) добавляется еще одно слагаемое, так называемая регуляризационная функция. Это слагаемое должно равняться нулю при правильных значениях весов модели и быть больше нуля при неправильных значениях весов. И, желательно, чтобы функция была дифференцируемой (хотя, конечно, в случае известного регуляризатора L1 это требование нарушено).

Метод градиентного спуска старается свести функцию потерь к нулю, следовательно он будет стараться свести к нулю и функцию регуляризации, являющуюся слагаемым в функции потерь, следовательно он будет стараться выполнить ограничение на веса модели, заданное этой функцией.

Функция регуляризации будет такая:

(10)

Здесь N - это требуемое нейтральное среднее значение весов, зависящее от способа объединения убывающей и периодической функций. 

Если функции объединяются аддитивно, N = 0.

Если функции объединяются мультипликативно, N = 1.

Как учесть меняющуюся дисперсию

Дисперсия биномиального распределения ретеншна рассчитывается по формуле (5), повторю ее длянаглядгнсти

(5)

где

I_i = I(t_i) - известное нам количество установок приложения в момент времени t_i

M(r_i) - математическое ожидание ретеншена в момент времени t_i

В качестве функции M(r_i) наиболее корректно использовать вычисленное значение аппроксимирущей функции, то есть:

M(r_i) ≈ R(t) = d(t_i) + s(t_i) + p(t_i) (аддитивная модель) (11)

или

M(r_i) ≈ R(t) = (d(t_i) + s(t_i)) * p(t_i) (мультипликативная модель) (12)

На начальных итерациях градиентного спуска такая величина может быть очень далека от реального значения, она может даже не попадать в диапазон [0,1], поэтому можно попробовать использовать величину r(t).

M(r_i) ≈ r_i (на начальных итерациях градиентного спуска) (13)

Итак, дисперсия в каждой точке известна, поэтому можно использовать так называемый взвешенный метод наименьших квадратов. Как указано в [6] функция потерь этого метода имеет вид:

(14)

Где e(t_i) - ошибка, поученная как разность аппроксимирующей и аппроксимируемой функции.

Подставив сюда формулу дисперсии, получим

(15)

Где

{ t_i }_, {r_i } = аппроксимируемый временной ряд значений ретеншена

{ t_i }_, { I_i } - временной ряд соответствующего количества установок

R(t) - аппроксимирующая функция, вычисляемая по формуле (8) или (9)

M(r_i) - математическое ожидание ретеншена в точке i.

Заметим, что значение, получаемое с помощью этой формулы, будет зависеть от количества установок каждый день и от количества дней наблюдений. Это плохо, поскольку мы будем суммировать эту величину с функцией регуляризации весов сезонной функции, не зависящей от количества установок и рассматриваемого промежутка времени. Поэтому разделим полученное значение на общее количество установок за рассматриваемый период и получим итоговую нормированную формулу:

(16)

Реализация. Пользовательский интерфейс

Рассмотренная выше методика реализована мною на Python + PyTorch в виде колаб ноутбука retention-rate-approximator, который я выкладываю в открытый доступ на github.

В первую очередь коротко распишу его структуру, чтобы вы могли уже его потрогать руками.

Ноутбук состоит из трех блоков:

1. "В первом блоке "Import libraries, define classes and functions" содержится описание классов и функций, которые используются ниже. Вам следует просто выполнить этот блок, наивно надеясь, что он не сольет в сеть ваши персональные данные.

2. Второй блок - пользовательский интерфейс. Этот блок надо тоже выполнить, чтобы этот интерфейс увидеть. Интерфейс состоит из двух частей. Слева - генератор временных рядов, похожих на ретеншен, справа - настройки аппроксиматора. Все параметры уже инициализированы некоторыми не абсурдными начальными значениями, поэтому при первом запуске вы можете просто нажать кнопку Genetate dataset в левой части экрана. Также можно сохранить сгенерированные данные в csv-файл, нажав кнопку Save dataset и загрузить свои собственные обучающие данные в формате csv, нажав кнопку Upload dataset

3. Третий блок "Create the regression model and train it" - собственно создание и обучение аппроксимирующей модели. Если вы заполнили все параметры параметры аппроксимирующей функции в правом блоке настройки аппроксиматора, просто запустите этот блок.

Реализация. Технические детали

Не буду подробно расписывать код ноутбука, он довольно простой и немного прокомментирован, все видно наглядно.

Конечно в процессе возникло работы, как обычно, множество нюансов, все не упомнишь, перечислю лишь несколько моментов.

1. Использован PyTorch а не Tensorflow/Kers поскольку первый более удобен для создания кастомных регуляризаторов и для создания кастомных моделей со сложными конструкциями внутри функции forward, что мне и требовалось.

2. Хотелось добиться в результате некоего подобия продукта, легко распространяемого и удобного в использовании, поэтому был выбран формат Jupyter-ноутбука, который можно запустить на любом компьютере, в любом браузере в Google Colab. Пришлось помучаться с пользовательским интерфейсом на базе библиотеки ipywidgets, с которой работаю впервые. Результатом не вполне доволен, но, грубо, желаемый результат достигнут - получилось доступно и более-менее удобно.

3. Весь код собран в одном ноутбуке и это конечно очень громоздко и неудобно в разработке. Следовало бы отделить 1-2 py-файла, держать их на гитхабе и загружать их в первом блоке ноутбука. Но я таких больших ноутбуков еще не создавал, и не знаю, как правильнее организовать такой проект. Возможно, кто-нибудь даст дельный совет. 

4. Помимо собственно аппроксиматора в ноутбуке есть генератор процессов, похожих на ретеншен, о чем я уже писал выше. Генератор довольно простой, и использует для генерации ту же самую математическую модель ретеншена, которая используется для аппроксимации.

5. Я попробовал реализовать и поэкспериментировал с несколькими вариантами функции для аппроксимации линии тренда, однако, больше всего работал только дробно-линейными функциями y(x) = w0-w1*x/(w2+x) (класс LinearFractionalFunction) и y(x) = w0-(w0-w1)*x/(1/w2+x)(класс LinearFractionalFunction_new). Именно эти классы реализованы лучше других. В частности, для качественного обучения модели имеет смысл предварительно инициализировать веса примерно подходящими подходящими значениями, используя данные обучающей выборки, этим занимается метод init_weights_from_train_data. Этот метод более-менее неплохо реализован в классах LinearFractionalFunction, LinearFractionalFunction_new и SigmoidFunction, но в остальных классах - совсем грубо.

6. Настройки стратегии обучения модели в текущий момент еще не вынесены в пользовательский интерфейс, их необходимо менять непосредственно в коде третьего блока ноутбука в переменных training_strateg, loss_function, regualizer_lambda,
   number_of_sigmas_for_plotting. Что означают эти переменные написано в комментариях в коде.

7. Я пробовал экспериментировать с различными типами оптимизаторов, в том числе использовал LBFGS, часто рекомендуемый для нелинейной регрессии, однако самый распространенный и популярный оптимизатор Adam, по моим наблюдениям, дает результаты как минимум не хуже большинства других. В настройках обучения, через константы, можно настроить оптимизатор Adam или LBFGS. Любой другой, разумеется тоже можно добавить, поменяв одну строчку в коде.

Эксперименты с функцией линии тренда

Выше был приведен список возможных кандидатов в функции аппроксимации тренда, он не много пугает разнообразием. Сейчас его немного отфильтруем.

Явно следует сразу исключить экспоненту - она убывает слишком быстро, совершенно не соответствует природе процесса. Гиперболический тангенс тоже выглядит не слишком удачным решением - у него есть два практически прямых участка - центральный, наклонный, и горизонтальная асимптота, соединенные перегибом - это слабо похоже на реальный процесс. Насчет других вариантов настолько однозначно сказать уже нельзя, нужно экспериментировать.

У меня нет, к сожалению, доступа к большому объему данных реальных проектов чтобы делать выводы, претендующие на объективность но по моим субъективным наблюдениям неплохо ведет себя обратная функция вида d(t) = w₀ + w₁ / (w₂ +t) и дробно-линейная функция вида d(t) = w₀ + w₁ * t / (w₂ +t).

В процесс экспериментов я остановился на функции вида d(t) = w₀ + w₁ * t / (w₂ +t), но далее, в процессе многих последующих экспериментов немного ее видоизменил, выполнив небольшое преобразование ее коэффициентов и получил итоговую функцию, которую считаю наиболее удачной:

d(t) = w₀ + ( w₀ - w₁ )* t / (1/w₂ +t ) (17)

С точки зрения зависимости от входной переменной это так же формула
d(t) = w₀ + w₁ * t / (w₂ +t), но для градиентного спуска важна зависимость не от входа, а от параметров. По поим наблюдениям, такой вид функции благотворно сказывается на результатах градиентного спуска, и у меня есть для этого объяснение. Аргументы в пользу функции такие:

1. Мало весов (3, а не 4 и не 5), что важно при небольшом объеме сильно зашумленных данных.

2. Не слишком быстро убывает, и, по моим наблюдениям, неплохо аппроксимирует ретеншн в течении довольно продолжительного количества дней после релиза.

3. Хороший физический смысл каждого параметра и они независимы друг от друга:

   1. Коэффициент w₀ задает начальное значение функции в нулевой момент времени.

   2. Коэффициент w₁ задает значений функции на бесконечном значении аргумента.

   3. Коэффициент w₂ задает скорость убывания функции.

4. Все 3 параметра имеют значения примерно одного порядка в диапазоне от 0 до 1. По моим наблюдениям, это свойство позволяет оптимизатору наиболее точно находить минимум функции ошибки.

5. При определенных значениях параметров превращается в линейную функцию или константу.

6. Вычислительная простота.

Вовсе не утверждаю, что такая функция единственно правильная. У меня даже нет достаточного количества данных для анализа, чтобы делать такие выводы. Если исследователь располагает большой историей развития метрики за большой промежуток времени, можно попробовать функции с 4-ю или 5- коэффициентами, можно составить функцию из нескольких звеньев. Если же данных наоборот слишком мало единственным возможным вариантом является линейная функция или, даже, константа.

Также я провел эксперименты с разными видами функции аппроксимации патчей.

При рассматриваемых условиях - ежедневное количество установок приложения ~1000 хороший результат можно получить только при использовании простейшего варианта функции p_i(t) = w0_i

Если использовать линейную функцию с двумя коэффициентами или еще более сложную функцию, получается полный провал получается полная ерунда - график ведет себя совершенно противоестественно. Необходимо гораздо больше инсталлов для такой аппроксимации. 

Грубая оценка доверительного интервала

А сейчас самое интересное - оценка доверительного интервала. Ради чего все и затевалось.

Непростая задача, совсем непростая. Регрессия нелинейная, теория тут сложная. Данных в обучающей выборке мало, поэтому не может быть и речи о том, чтобы ее разделить на обучающую тестовую и валидационную. У меня нет доступа к большому объему данных от множества проектов, чтобы проверить методику на большом объеме данных. Кроме того, существует огромное количеств факторов, влияющих на оценку доверительного интервала, это просто трудоемко. Например, известна или нет заранее форма сезонной функции, сколько дней прошло от момента релиза и т.д. Все эти варианты надо продумывать, прорабатывать. Большая работа, а времени у меня не так и много. И, кроме того, это не научное исследование, а практическая работа, надеюсь, полезная не только мне. Мне важнее опубликовать текущие результаты, чем тянуть время.

Ну вы поняли. Закапываться в жесткую теорию, не имея возможности проверить метод на большом объеме реальных данных совершенно глупо, не буду сейчас этого делать.

Но грубую оценку можно сделать, отталкиваясь от грубой теории и от фактических результатов измерений на моделированных данных.

Теория

Установим лучший идеальный доверительный интервал, который может быть получен при помощи такой регрессии. Доверительный интервал нашего метода не может быть лучше этого идеального.

Предположим, ретеншен не меняется в течении каждого патча, то есть линия тренда состоит из горизонтальных отрезков прямых - каждому патчу соответствует свой горизонтальный участок. Нам известна дисперсия в каждый момент времени (это биномиально-распределенная величина с известными параметрами), следовательно, полагая испытания в каждый день независимыми, мы можем посчитать дисперсию величины для каждого участка кривой (то есть для каждого патча) как взвешенное среднее значений дисперсии в каждый момент времени, а стандартное отклонение - как корень из этой дисперсии

(18)

где

k>=0 - номер патча

σ_k - оценка дисперсии k-го патча

n_k- количество измеренных дней в патче k

σ_kj - дисперсия в j-й день k-го патча (посчитанная как дисперсия биномиального распределения)

I_kj - количество установок приложения в j-й день k-го патча

В эту формулу, конечно, должны входить только те точки, которые использовались для аппроксимации. Точки, отброшенные как аномальные, необходимо исключить.

Таким образом, можно построить доверительный интервал аппроксимирующей кривой в каждой точке наблюдений. На графике это доверительный интервал изображается в виде полосы вокруг полученной аппроксимирующей кривой. Ширину полосы можно установить равной 2 сигмы (что должно соответствовать 95% точности аппроксимации) или 3 сигмы (точность аппроксимации - 99.7%)

Таким образом, лучшая точность, которую мы можем требовать от нашего метода аппроксимации,  - это стопроцентное попадание в такой коридор шириной 3 сигмы. Более высокой точности требовать нельзя.

На рисунке 4 коридор такого доверительного интервала показан полосой серого цвета. Видно, что ширина коридора на различных участках различна. Чем шире участок, тем больше в ней точек наблюдений, поэтому коридор уже.

Практика

Не могу проверить адекватность предложенной модели применительно к большому объему реальных данных. Это плохо, но факт. Но можно проверить модель на аппроксимации искусственных данных, сгенерированных при помощи функции, аналогичной аппроксимирующей с добавлением биномиально-распределенного шума. Преимуществом такого метода является точное знание изначальной формы процесса, следовательно, можно точно оценить результат аппроксимации.

Я  провел моделирование 1000 случайных процессов, и оценил их результаты. Во всех случаях я предполагал, что интервал между патчами - около 30 дней, и с момента релиза прошло довольно много времени, от 32 до 365 дней. Во всех экспериментах среднее количество установок в день равнялось 1000. Во всех экспериментах 30% случайных данных от выборки а также все даты старта патчей и дата релиза было отброшено из выборки как якобы аномальные. Во всех случаях сезонная функция была примерно известна заранее. Довольно много допущений, но они не очень недалеки от практики

В качестве генерирующей и аппроксимирующей функции линии тренда использовал функцию вида d(t) = w₀ + ( w₀ - w₁ )* t / (1/w₂ +t )

В качестве линии каждого патча использовал простейшую функцию вида p_i(t) = w0_i

В качестве критерия оценки качества модели  я рассматривал простой бинарный факт - попадает ли кривая полностью в идеальный доверительный интервал 3-сигмы, полученный при помощи описанного выше метода или нет.

Полученные результаты 

1. Во всех испытаниях аппроксимируемая функция попала в идеальный доверительный интервал 3 сигмы в 100% случаях при аппроксимации всех патчей, кроме релизного.(при условии, что хотя бы одна точка наблюдений использовалась для аппроксимации патча - в редких случаях все точки патча отбрасывались как анамальные).

2. Примерно в 25% экспериментов кривая выходит из коридора доверительного интервала в первом сегменте, в начальных 2-3 точках.

Результаты грубо объясняю так:

1. В начальной стадии кривой на ее форму влияет все 3 основных веса, но, поскольку точек в этой начальной части кривой достаточно мало, такого количества недостаточно для правильного подбора всех 3-х коэффициентов. Кроме того, из рассмотрения всегда была выброшена самая первая точка наблюдений - нулевой день, как, потенциально, анормальная, но именно эта первая точка должна вносить максимальный вклад в нелинейную форму кривой. Таким образом, для моделирования кривой на ее начальном участке, до первого патча, аппроксиматор использовал только ~20 точек, чего, вероятно, мало для правильного подбора 3-х коэффициентов.

2. Что же касается аппроксимации патчей - уже через месяц после выпуска игры наиболее существенное значение на форму тренда имеет ее асимптотическое значение, а также коэффициент патча, и, следовательно, процесс приближен к рассмотренному выше идеальному случаю, и, следовательно, доверительный интервал определяется весьма точно.

Наверняка, имеют место и традиционные проблемы нелинейной регрессии - несколько локальным минимумов функции потерь, а также близость функции градиента к нулю в широком диапазоне значений параметров.

Все это , конечно, не настоящее исследование, это все сделано на глазок, но я не вижу смысла для столь грубой модели проводить исследование более полноценное без экспериментов над большими объемами реальных данных.

Критика

Самое время поговорить о недостатках методики. Они очевидны.

Необходимо весьма много данных, за большое количество дней наблюдений, чтобы аппроксимируемый процесс становился похож на реальный. Если, скажем, мы располагаем только лишь неделей наблюдений после релиза проекта, то во многих случаях о нелинейной регрессии не может идти и речи, лучше использовать обычную линейную. Или даже скалярную оценку. Попытка применить к такому процессу нелинейную регрессию приводит к неожиданным значениям коэффициентов. Ниже я привел 2 примера неверной нелинейной аппроксимации. В первом примере мы видим огромный провал в нулевой точке. Во втором примере нет такого гигантского провала, но наблюдается возрастающий процесс. Оба примера являются следствием того, что для аппроксимации были использованы только 4 из 7-ми точек - нулевая точка и еще 2 случайных были отброшены как аномальные.

Это вполне реалистичные кейсы. После выхода игры в софтланч вполне возможны сбои и действительно из 7-ми дней наблюдений у нас вполне может остаться только 4 хороших дня. Количество инсталлов также вполне может быть около 1000 в день, как в рассмотренных примерах. 

Хотелось бы иметь модель, обеспечивающую более точные измерения даже в таких сложных ситуациях.

Что дальше?

Существует ли возможность повысить точность наблюдений?

Да, несомненно. Во введении я уже упоминал несколько потенциальных способов, и сейчас остановлюсь на этом вопросе чуть подробное.

1. Можно улучшить текущую методику, тщательнее проработав детали реализации, тщательнее подобрать гиперпараметры модели. Например, использовать другую, более точную, аппроксимирующую функцию, поэкспериментировать с разными оптимизаторами и видом функции потерь, более точно инициализировать начальные веса модели, добавить какие-нибудь эвристические ограничения на значения весов и т.д. 

2. Более точно оценить вид исходного распределения и его дисперсию и за счет этого сузить доверительный интервал. Например, известно, что дисперсия суммы нескольких независимых биномиально-распределенных величин меньше дисперсии одного биномиального распределения (статья в википеии), этот факт можно попробовать использовать. Ретеншен действительно легко представить в виде суммы нескольких случайных величин, каждая из которых соответствует одной стране или региону, однако независимость этих величин будет под вопросом, это еще следует доказать. Кроме того, этот метод даст заметный эффект только если матожидания у разных сегментов сильно отличаются между собой, поэтому общий выигрыш будет, вероятно, небольшим. Но зато он достается практически бесплатно, без какого-либо усовершенствования метода, только теоретическим обоснованием.

Кстати, я не встречал никаких работ, посвященных исследованию функции распределения ретеншена. Ни теоретических , ни практических. Если у читателей есть информация о подобной работе, был бы очень рад получить ссылочку.

3. Следующий, потенциальный способ повысить точность - использовать для аппроксимации какие-то дополнительные данные. Например, пусть мы хотим аппроксимировать ретеншен 1 дня, располагая данными за 7 дней. Однако, в таком случае, у нас есть также данные о ретеншене 2-го дня за 6 дней измерений, данные о ретеншене 3 дня за 5 дней… , данные о ретеншене 7-го дня за 1 день. Таким образом, можно говорить о ретеншене как о двумерной функции, в координатах (время от релиза, время от инсталла) или как о некоторой поверхности в 3d пространстве. Форма этой поверхности будет довольно сложной, поскольку на нее будут накладываться и сезонные явления и выпуск патчей и аномалии, кроме того существует явная корреляция между ретеншеном N дня и ретеншеном M дня. Но, зато, такой способ позволяет увеличить количество точек наблюдения квадратично! Линия общего тренда такой поверхности должна быть непрерывной и, следовательно, может быть представлена моделью с небольшим числом коэффициентов, количество которых возрастет не квадратично.

   Можно еще выше поднять размерность функции ретеншена и добавить еще больше обучающих данных, рассматривая не классический ретеншен, а роллинг ретеншен ограниченной глубины. Действительно, роллинг ретеншен - это 3d-функция в координатах (время от релиза, время от инсталла, глубина роллинга). Классический ретеншен тут является частным случаем роллинга, когда глубина роллинга равняется нулю. Такую трехмерную модель еще сложнее построить, необходимо преодолеть еще больше проблем, однако количество исходных данных увеличивается уже в кубе!

4. Следующий потенциальный способ поднять точность измерений - учесть данные, наблюдаемые в момент аномалий. Выше я отбросил все анормальные дни, однако, можно рассмотреть эти данные более детально. Тщательный анализ аномалий может позволить сохранить как минимум часть наблюдений. Например, можно заметить, что технический сбой длится не весь день, и охватывает не всех пользователей, и как минимум часть данных в аномальные дни мы сможем использовать. 

5. Построить сложную модель, использующую данных с других проектов.
   Вероятно, это самый мощный способ создать максимально точную аппроксимацию. Крупные паблишеры и аналитические компании располагают огромными датасетами, накопленными от проектов, ранее выпущенных в релиз. Наличие такого датасета позволяет применять сложные модели с большим количеством весов. С этой точки зрения задача аппроксимации ретеншена представляется уже как задача восстановления из помех некоторого сигнала с известными нам свойствами и подобна, например, задаче восстановления зашумленной картинки. Можно построить модель, решающую подобную задачу и обучить его на этом датасете. Такая модель будет способна спрогнозировать значение метрики проекта, возможно, уже по одной двум точкам наблюдения сразу после релиза.

6. Использовать готовые фрймворки для предсказания временных рядов, например Prophet, Neural Prophet, Nixtla/statsforecast. Я только лишь поверхностно знаком с возможностями Prophet и пока не берусь судить, насколько он хорош для данной задачи. С другими фреймворками я пока совсем не сталкивался. В лоб получить заметный эффект Prophet не получилось, ведь он не знает ничего о природе нашего процесса. Однако, если можно задать ему правильный вид нашей аппроксимирующей функции, мы несомненно получим более качественное предсказание, как минимум потому, что он хорошо умеет учитывать сезонные явления - не только недельные, но и годовые, а также нацелен на нелинейную аппроксимацию и учитывает ее тонкости.  Способен ли Prophet на большее, например, рассматривать процесс как двумерную функцию или трехмерную функцию, я пока не знаю, но планирую разобраться.

Заключение и выводы

Итак, подведу краткие итоги

1. Рассмотрена методика нелинейной регрессии метрики "ретеншен", учитывающий недельную сезонность и выпуск патчей.

2. Произведены эксперименты с различными видами аппроксимирующей функции и сделана субъективная рекомендация о выборе вида функции.

3. Методика реализована в виде ноутбука Jupyter-ноутбука retention-rate-approximator, выложено в открытый доступ на Githab

4. Дальнейшие планируемые направления работ - использовать для решения задачи специализированные фреймворки для анализа временных рядов, построить нейросеть с большим количеством весов, обучить ее на большом объеме данных (возможно, частично, искусственно сгенерированных), построить модель для двумерной и трехмерной аппроксимации процесса.

5. Отдельно следует рассмотреть задачу обнаружения аномалий, а также вопросы использования в аппроксимации метрики данных, полученных в аномальные дни.

6. При разработки более сложных моделей уделить больше внимания построению интервальной оценки.

Источники

1. Binomial sum variance inequality https://en.wikipedia.org/wiki/Binomial_sum_variance_inequality

2. Василий Сабиров: Игра в цифры. Как аналитика позволяет видеоиграм жить лучше

3. В.П. Лисьев ТЕОРИЯ ВЕРОЯТНОСТЕЙ И МАТЕМАТИЧЕСКАЯ СТАТИСТИКА: Учебное пособие/ Московский государственный университет экономики, статистики и информатики. – М., 2006. – 199 с.

4. Дмитрий Балтин. Ретеншен — основная метрика F2P игры, вероятностный подход. https://habr.com/ru/articles/726396/

5. К.В. Воронцов,  Машинное обучение. Нелинейная регрессия. Школа анализа данных, Яндекс. https://www.youtube.com/watch?v=--gJR8pd-jg 

6. Филипп Картаев. Дружелюбная эконометрика https://books.econ.msu.ru/Introduction-to-Econometrics/

7. Биномиальное распределение https://ru.wikipedia.org/wiki/Биномиальное_распределение

8. Метрика Retention. Что означает, как ее рассчитать и как ее улучшить https://gopractice.ru/product/retention/

9. Нормальное распределение https://ru.wikipedia.org/wiki/Нормальное_распределение

Центральная предельная теорема https://ru.wikipedia.org/wiki/Центральная_предельная_теорема