Цена качества модели05.03.2024 05:00

Описание проблемы

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

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

Например, мы хотим выдать автокредиты 100 000 потенциальных клиентов, по 1 млн. рублей каждый. Допустим мы готовы отсечь 10% потока, т.е. выдадим кредиты 9 обратившимся из 10. Тогда при DR = 5%, мы потеряем в результате дефолтов 5% * 100 000×1 000 000 = 5 млрд. рублей, если же удастся снизить уровень DR до 3%, то потери будут уже только 3% * 100 000×1 000 000 = 3 млрд. рублей, т.е. чистая прибыль от такого снижения DR будет равна 2 млрд. рублей. Получается, что если наше потенциальное изменение модели может обеспечить такое снижение уровня дефолтов оно принесет чистую прибыль в 2 млрд. рублей в нашем случае. Это может быть уже достаточным аргументом, чтобы нанять отдельного data scientist-а (или целую команду), который будет работать над этой задачей. Нужно понимать, что зачастую мы никак не можем повлиять исходный уровень дефолтов (качество потока) и на то, какую часть клиентов мы отбрасываем, все эти решения обычно принимаются бизнесом на основе того, какую прибыль они хотят получить, какую долю какого рынка захватить и т.п. Единственное на что мы как data scientist-ы мы можем влиять, это качество модели, которое в целом описывается ее коэффициентом Джини (при условии, что модель устойчива out-of-time и этот показатель вообще имеет смысл).

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

Формальная постановка задачи

Обратим ваше внимание, что за хотя модели классификации формально предсказывают вероятность (число от 0 до 1), зачастую за этой вероятностью стоит скор балл (величина от минус бесконечности до плюс бесконечности), который просто превращается в вероятность с помощью логистического преобразования. Соответственно его можно получить, совершив обратное логистическое преобразование. Во многим ситуациях скор балл распределен нормальным образом, но если это не так нормальности можно добиться совершив монотонное преобразование скора (в нашем случае важен только порядок, в котором скор ранжирует клиентов, а не сам скор). Достаточно стандартными является предположение, что кроме скора полученного моделью, существует скор нам не известный и в совокупности они образуют «идеальный» скор, который позволяет полностью отделить «хороших» клиентов от «плохих» (т.е. у всех «плохих» «идеальный» скор будет меньше, чем у «хороших»).

Более формально, давайте рассмотрим следующую модель.

Пусть и $\varepsilon$ независимые случайно распределенные случайные величины N(0, 1) , это скор генерируемый моделью машинного обучения, $\varepsilon$ это некоторая неучтенная информация, которая тем не менее влияет на результат

$z = \alpha x + \sqrt{1-\alpha^2}\varepsilon$ , где $\alpha$ некоторая константа (в данном случае также $z \sim N(0,1)$ ), а z скрытая переменная, которая определяет результат,

Давайте также возьмем некоторый трешхолд , и когда при z < t мы можем наблюдать некоторое событие (дефолт заемщика, болезнь, разрушение конструкции и т.п.) Вероятность данного события

$p = \Phi(t) = \int_{-\infty}^t\frac{\exp(-t^2/2)}{\sqrt{2\pi}}$ , при этом y = I(z < t) это наблюдаемая случайная

величина, равная 1 если z < t и нулю, если $z \geq t$ .

Используя скор и наблюдаемую величину мы можем рассчитать AUC и коэффициент Джини = 2 * AUC — 1.

Давайте попытаемся понять как $\alpha$ и будут влиять на AUC и коэффициент Джини.

Связь между коэффициентом Джини и корреляцией

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

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

import numpy as np
from scipy.stats import norm
import sklearn
from sklearn import metrics
import matplotlib.pyplot as plt
x = np.random.normal(size = 1000000)
e = np.random.normal(size = 1000000)

alpha_array = (np.arange(20) + 1) / 20
p_list = [0.1, 0.5, 0.7]
gini_list = []
for p in p_list:
    t = norm.ppf(p)
    gini_array = alpha_array * 0
    for i, alpha in enumerate(alpha_array):
        z = alpha * x + np.sqrt(1 - alpha * alpha) * e
        gini_array[i] = 2 * sklearn.metrics.roc_auc_score(z < t, -x) - 1
    gini_list.append(gini_array)

plt.rcParams['figure.figsize'] = [10, 10]
plt.plot(alpha_array, alpha_array, label = 'alpha')
for i, gini_array in enumerate(gini_list):
    plt.plot(alpha_array , gini_array, label='p = ' + str(p_list[i]))
plt.legend(loc="upper left")
plt.xlabel('alpha')
plt.ylabel('Gini')

Мы можем видеть, что почти нет влияния различных и трешхолдов на Джини, но в тоже самое время Джини почти в точности равен коэффициенту $\alpha$ .

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

Согласно вероятностной интерпретации AUC в нашем случае:
AUC = P(x_1 > x_2 | z_1 > t > z_2)» src=«https://habrastorage.org/getpro/habr/upload_files/a36/100/d7f/a36100d7ff64613fba5a3e6df9dadc08.svg» />, где <img alt= и x_2, z_2 два независимых реализации x, z определенных выше (с теми же параметрами $\alpha, t, p$ ). Давайте предположим, что x_1-x_2 где z_1 > t» src=«https://habrastorage.org/getpro/habr/upload_files/227/1c6/f15/2271c6f157c74ed2869c5da4ae97f6ac.svg» /> и <img alt= нормально распределенная случайная величина. В этом случае мы должны вычислить ее среднее и дисперсию, чтобы вычислить соответствующие вероятности.

Мы можем видеть, что
$(1-p)E(x|\alpha x + \sqrt{1-\alpha^2}\varepsilon > t) = \int dx d\varepsilon I (z > t) x \frac{\exp (-x^2/2-\varepsilon^2/2)}{2\pi} =» src=«https://habrastorage.org/getpro/habr/upload_files/2cd/c89/ca4/2cdc89ca47c2d33ef62fb879cf4c7666.svg» /> <img alt=$ t) x \frac{\exp\Big (-\frac{x^2+z^2–2xz\alpha}{2(1-\alpha^2)}\Big)}{2\pi}=» src=«https://habrastorage.org/getpro/habr/upload_files/cc6/d19/17b/cc6d1917b431624a9e3b6773d989e478.svg» />

$=(1-\alpha^2) \int dx dz I\Big(z> \frac{t}{\sqrt{1-\alpha ^2}}\Big) x \frac{\exp (-x^2/2-z^2/2)}{2\pi} \exp (xz\alpha)\approx» src=«https://habrastorage.org/getpro/habr/upload_files/405/19b/10d/40519b10d37c11695cfc1033ee03d31f.svg» /> <img alt=$ t) \frac{\exp (-x^2/2-z^2/2)}{2\pi} (x+x^2 z \alpha)=\alpha \int dz I (z > t) z \frac{\exp (-z^2/2)}{\sqrt{2\pi}} =» src=«https://habrastorage.org/getpro/habr/upload_files/ef0/533/d6c/ef0533d6c13e7820b6c55def2d34785d.svg» />

$=\alpha\frac{e^{-t^2/2}}{\sqrt{2\pi}}$

так как при достаточно малых $\alpha$ , $1-\alpha^2\approx 1$ , $\exp(x z \alpha)\approx 1 + x z \alpha$

таким же образом

$(1-p)E(x^2|\alpha x + \sqrt{1-\alpha^2}\varepsilon > t) \approx \int dx dz I (z> t) \frac{\exp (-x^2/2-z^2/2)}{2\pi} (x^2+x^3 z \alpha) =» src=«https://habrastorage.org/getpro/habr/upload_files/25f/914/9e4/25f9149e43efb26448fa54ddc11c56e1.svg» /> <img alt=$ t) \frac{\exp (-z^2/2)}{\sqrt{2\pi}} = 1-p» src=«https://habrastorage.org/getpro/habr/upload_files/1c9/d05/e63/1c9d05e636eb0f11d1038b1be38b2217.svg» />

поэтому

$E(x|\alpha x + \sqrt{1-\alpha^2}\varepsilon > t) \approx \frac\alpha{1-p}\frac{e^{-t^2/2}}{\sqrt{2\pi}}» src=«https://habrastorage.org/getpro/habr/upload_files/10b/3be/7e0/10b3be7e039af184e1e88a107e9be524.svg» /> <img alt=$

$D(x|\alpha x + \sqrt{1-\alpha^2}\varepsilon > t) \approx 1» src=«https://habrastorage.org/getpro/habr/upload_files/637/b41/ddd/637b41dddd4308d5c62d60c2eef49117.svg» /> <img alt=$

следовательно

$E(x_1 - x_2|z_1 > t > z_2) \approx \frac{\alpha}{p (1-p)} \frac{e^{-t^2/2}}{\sqrt{2\pi}}» src=«https://habrastorage.org/getpro/habr/upload_files/6a6/1de/647/6a61de647eb2a250e3a0342d510f0f2f.svg» /> <img alt=$ t > z_2) \approx 2» src=«https://habrastorage.org/getpro/habr/upload_files/795/41b/c5d/79541bc5dc52f4a7e8232ca6b7912e15.svg» />

и

$x_1 - x_2 \sim N\Bigg(\frac{\alpha}{p(1-p)} \frac{e^{-t^2/2}}{\sqrt{2\pi}}, 2\Bigg)$

тогда

$P(x_1 > x_2 | z_1 > t > z_2) \approx \Phi \Bigg (\frac{\alpha}{p (1-p)\sqrt{2}} \frac{e^{-t^2/2}}{\sqrt{2\pi}} \Bigg)\approx» src=«https://habrastorage.org/getpro/habr/upload_files/027/a90/363/027a903633ee74e8789d59ffafa33930.svg» /> <img alt=$

так как при достаточно малых , $\Phi(x) \approx \frac12 + \frac1{\sqrt{2\pi}}x$

Заметим также, что

$\Bigg(\frac{e^{-t^2/2}}{\Phi(t)\Phi(-t)}\Bigg)' = \frac{e^{-t^2/2}}{\Phi^2(t)\Phi^2(-t)}\Bigg(t\Phi(t)\Phi(-t)+\Phi(t)-\Phi(-t)\Bigg) > 0 \iff t > 0» src=«https://habrastorage.org/getpro/habr/upload_files/256/c4e/cc8/256c4ecc8f4d29a5c368af9ea2a1858a.svg» /> Значит при <img alt=$

Например при $p \in [0.05;0.95], \frac{\exp(-t^2/2)}{p(1-p)\pi\sqrt{2}}\in[0.9;1.225]$

Поэтому

$AUC = P(x_1 > x_2 | z_1 > t > z_2) \approx \frac{1}{2} + \frac{\alpha}{2}» src=«https://habrastorage.org/getpro/habr/upload_files/140/d1d/b5a/140d1db5a7223b88ac68deec3fa6e22c.svg» /> Поэтому мы можем видеть, что коэффициент Джини <img alt=$ . Значит мы можем интерпретировать коэффициент Джини как коэффициент корреляции между скрытой переменной и скором модели машинного обучения. Например, для соответствующей регрессионной модели коэффициент детерминации $R^2=\alpha^2 = Gini^2$ . Это показывает какую часть информации о скрытой переменной мы знаем благодаря нашей модели, в противоположность информации, которую мы не знаем, но она все еще влияет на нашу скрытую переменную.

Формула вычисления стоимости Джини

Вернемся теперь к исходной задаче определения уровня дефолтов по Джини, исходному уровню дефолтов и порогу отсечения. В данном случае, под порогом отсечения будем понимать трешхолд для скора t_x такой, что $p_x = \Phi(t_x)$ — доля самых «плохих» клиентов, которых мы отбросили.

$E(y|x > t_x)= P (z < t | x\geq t_x) = P\Big(\varepsilon < \frac{t-\alpha x}{\sqrt{1-\alpha^2}}| x \geq t_x\Big) =$

$=\frac1{1-p_x}\int_{t_x}^{\infty} \Phi\Big(\frac{t-\alpha x}{\sqrt{1-\alpha^2}}\Big)\frac{e^{-x^2/2}}{\sqrt{2\pi}}dx$

Как нам удалось установить ранее $Gini\approx \alpha$ и поэтому

$E(y|x > t_x)\approx \frac1{1-p_x}\int_{t_x}^{\infty} \Phi\Big (\frac{t-Gini\cdot x}{\sqrt{1-Gini^2}}\Big)\frac{e^{-x^2/2}}{\sqrt{2\pi}}dx» src=«https://habrastorage.org/getpro/habr/upload_files/387/0b4/f04/3870b4f04f09e597b5fbf58f55f4792c.svg» /> Этот интеграл достаточно легко находится численно, например, с помощью следующей функции <pre><code class=$ def DR_new(gini, p, p_x): t = norm.ppf(p) t_x = norm.ppf(p_x) dx = 0.001 x = np.arange(int((5 - t_x) / dx)) * dx + t_x return np.sum((norm.cdf((t - gini * x) / np.sqrt(1 - gini * gini))) * norm.pdf(x) * dx) / (1 - p_x)

Проверим, теперь что формула соответствует прямым вычислениям достаточно точно. В данном случае возьмем p_x = p , т.е. мы отбрасываем такую же пропорцию клиентов, какое количество клиентов у нас есть в потоке. В этом случае при Джини = 100% мы отбросим всех «плохих» клиентов и уровень дефолтов будет равен 0.

import numpy as np
from scipy.stats import norm
import sklearn
from sklearn import metrics
import matplotlib.pyplot as plt
x = np.random.normal(size = 1000000)
e = np.random.normal(size = 1000000)

alpha_array = (np.arange(20) + 0.999) / 20
p_list = [0.1, 0.5, 0.7]
gini_list = []
DR_real_list = []
DR_calculated_list = []
for p in p_list:
    t = norm.ppf(p)
    gini_array = alpha_array * 0
    DR_real_array = alpha_array * 0
    DR_calculated_array = alpha_array * 0
    for i, alpha in enumerate(alpha_array):
        z = alpha * x + np.sqrt(1 - alpha * alpha) * e
        gini_array[i] = 2 * sklearn.metrics.roc_auc_score(z < t, -x) - 1
        DR_real_array[i] = np.mean(z[x > t] < t)
        DR_calculated_array[i] = DR_new(gini_array[i], p, p)
    gini_list.append(gini_array)
    DR_real_list.append(DR_real_array)
    DR_calculated_list.append(DR_calculated_array)

plt.rcParams['figure.figsize'] = [10, 10]
for i, gini_array in enumerate(gini_list):
    plt.plot(gini_array, DR_real_list[i], label='DR real p = ' + str(p_list[i]))
    plt.plot(gini_array, DR_calculated_list[i], label='DR calculated p = ' + str(p_list[i]))
plt.legend(loc="upper left")
plt.xlabel('Gini')
plt.ylabel('DR')

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

Приведем практический пример, предположим нам нужно выдать 100 000 автокредитов, по 1 млн. рублей каждый, т.е. всего нужно выдать 100 млрд. рублей. Допустим изначальный DR по данному потоку клиентов 10%, при этом мы готовы отсеивать 10% самых худших клиентов (получается, что в идеальном мире у нас вообще не было бы дефолтов). Пусть у нас есть модель, согласно которой мы ранжируем клиентов, и она имеет Джини = 30%. В этом случае прямые потери от дефолтов будут 8,71 млрд. рублей. Допустим у нас есть план, согласно которому мы можем увеличить Джини на 1%, 5%, 10%, 50%, каким финансовым результатам это приведет? Согласно приведенным выше расчетам, это принесет экономию в результате уменьшения потерь от дефолтов в размере 52 млн. рублей, 270 млн. рублей, 559 млн. рублей и 3,846 млрд. рублей соответственно.

print(DR_new(0.30, 0.1, 0.1))
print(DR_new(0.30, 0.1, 0.1)-DR_new(0.30 + 0.01, 0.1, 0.1))
print(DR_new(0.30, 0.1, 0.1)-DR_new(0.30 + 0.05, 0.1, 0.1))
print(DR_new(0.30, 0.1, 0.1)-DR_new(0.30 + 0.10, 0.1, 0.1))
print(DR_new(0.30, 0.1, 0.1)-DR_new(0.30 + 0.50, 0.1, 0.1))

0.0871097134754345
0.0005272288318745877
0.0027060299128032483
0.005594040703079159
0.03845790575519499

Доказательство на реальном датасете

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

from ucimlrepo import fetch_ucirepo 
  
default_of_credit_card_clients = fetch_ucirepo(id=350) 
  
X = default_of_credit_card_clients.data.features 
y = default_of_credit_card_clients.data.targets 

from catboost import CatBoostClassifier
model = CatBoostClassifier(verbose = 0)
DR_initial = np.mean(y.values)
number_of_features = [4, 5, 23]
p_x_array = (np.arange(49) + 1)/100
DR_real_list = []
DR_calculated_list = []
gini_list = []
for n in number_of_features:
    model.fit(X[X.columns[:n]], y)
    pred = model.predict_proba(X[X.columns[:n]])[:, 1]
    gini = 2 * sklearn.metrics.roc_auc_score(y, pred) - 1
    gini_list.append(gini)
    DR_real_array = p_x_array * 0
    DR_calculated_array = p_x_array * 0
    for i, p_x in enumerate(p_x_array):
        t_pred = np.sort(pred)[-int(pred.shape[0] * p_x)]
        DR_real_array[i] = np.mean(y.values[pred < t_pred])
        DR_calculated_array[i] = DR_new(gini, DR_initial, p_x)
    DR_real_list.append(DR_real_array)
    DR_calculated_list.append(DR_calculated_array)  

plt.rcParams['figure.figsize'] = [10, 10]
for i, gini_array in enumerate(gini_list):
    plt.plot(p_x_array, DR_real_list[i], label='DR real gini = ' + str(gini_list[i]))
    plt.plot(p_x_array, DR_calculated_list[i], label='DR calculated gini = ' + str(gini_list[i]))
plt.legend(loc="upper left")
plt.xlabel('$p_x$')
plt.ylabel('DR')

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

Выводы

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