Уточнение процентилей с помощью полиномиальной аппроксимации

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

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

Рассмотрим проблему на примере. Предположим, мы сгенерировали набор данных на 100 наблюдений, состоящий из дискретных значений с распределением от 1 до 10. Затем мы вычисляем процентили для трёх значений: 0.25, 0.50 и 0.75.

# зафиксируем псевдослучайность np.random.seed(45)  # сгенерируем данные с распределением от 1 до 10 data = pd.DataFrame(data=list(np.random.randint(1, 11, size=100)), columns=['values'])  # найдём длину оси X x_length = len(data['values'].value_counts())  # визуализируем полученные данные plt.hist(data, lw=2.5, edgecolor='#e5d5ff', color="#fbf9ff", bins=x_length) plt.xticks(np.arange(1, x_length + 1, step=1))  # визуализируем процентили for percentile in [0.25, 0.50, 0.75]:     plt.axvline(         x=data['values'].quantile(percentile), label=f'{percentile} Процентиль (не точный)', lw=2.5, color='#ffdb4d', linestyle='--')  # выключаем сетку plt.grid(visible=False)      # вызываем легенду plt.legend()  # вызываем график plt.show() 

Как мы видим, процентили приняли значения 4, 7 и 9. Однако отражают ли они реальное распределение вероятностей? Нет, так как исходные значения были дискретными, следовательно, и процентили будут иметь дискретный вид.

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

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

• Модель не объясняет изменчивость данных вообще, то есть предсказанные значения модели не коррелируют с реальными данными. Модель не лучше, чем простое среднее значение.

• Модель идеально объясняет изменчивость данных, т.е. предсказанные значения полностью совпадают с реальными.

• и Модель объясняет часть изменчивости данных. Чем ближе R² к 1, тем лучше модель описывает данные.

# создадим метрику для оценивания расхождений def r2(x, y):     return 1 - (np.mean((x - y) ** 2) / np.var(x))  # найдём наилучший полином по оценкам метрики def get_best_polynome(x, y):          best_result = 0     best_polynome = 0          # начинаем с 1, так как полином степени 0 не имеет смысла     for degree in range(1, 11):                  # создаем полином указанной степени         polynome = np.poly1d(np.polyfit(x, y, degree))                  # вычисляем метрику для текущего полинома         results = r2(y, polynome(x))                  # обновляем результат и степень полинома, если нашли лучший вариант         if results > best_result and results < 1:             best_result = results             best_polynome = polynome          return best_polynome, best_result  # формируем ось X x_data = range(1, x_length + 1)  # формируем ось Y (сортируем индексы, если они идут вразнобой) y_data = list(data['values'].value_counts().sort_index())  # найдём полином, описывающий распределение вероятностей, лучше всего best_polynome, best_score = get_best_polynome(x_data, y_data)  print(f'Лучший полином степени: {best_polynome.order}\nЛучшая оценка: {str(best_score)[:4]}')  # Лучший полином степени: 8 # Лучшая оценка: 0.99 

Насколько видно, полином, который описывает имеющиеся данные лучше всего, оказался со степенью 8. Если переводить на доступный язык, то степень полинома — это количество изгибов у линии, которые позволяют точнее описать данные. Так, если бы полином был первой степени, то получилась бы прямая, проходящая от точки A до точки B без единого изгиба.

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

# вычисляем коэффициенты полинома и создаём объект полинома на основе переданных коэффициентов coefficients = np.poly1d(np.polyfit(x_data, y_data, best_polynome.order))  # генерируем 10000 равномерно распределённых точкек от 1 до 10 polyline = np.linspace(1, x_length, 10000)  # визуализируем полученные данные plt.hist(data, lw=2.5, edgecolor='#ffffff', color="#fbf9ff", bins=x_length) plt.xticks(np.arange(1, x_length + 1, step=1))  # график, отражающий количество значений для каждого измерения plt.scatter(x=x_data, y=y_data, label=f'Изначальные данные', color='#ffdb4d')  # график, отражающий линию аппроксимации plt.plot(polyline, coefficients(polyline), lw=2.5, label=f'Сгенерированный полином', color='#c2a4f4')  # выключаем сетку plt.grid(visible=False)  # вызываем легенду plt.legend()  # вызываем график plt.show() 

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

1. Вначале необходимо вычислить значение полинома в точке x, для этого используем функцию «function()». Внутри цикла for функция суммирует вклад каждого члена полинома. Здесь используется массив coefficients, который содержит коэффициенты.

2. Далее нормализуем функцию «function()» и записываем её в «normalized_function()», чтобы интеграл от этой функции на заданном интервале был равен 1. Сначала интегрируем функцию function() на интервале от lower_limit до upper_limit, чтобы найти значение нормирующего коэффициента. Затем функция function() делится на этот нормирующий коэффициент, чтобы сумма всех вероятностей (интеграл функции) равнялась 1.

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

4. В заключении используем численный метод (алгоритм Брента, brentq) для нахождения такого значения x, при котором CDF(x) равна заданному процентилю. Иными словами, эта функция находит такое x, при котором интеграл функции от lower_limit до x равен значению процентили.

# вычислим значение полинома в точке x def function(x):          result = 0          # суммируем вклад каждого члена полинома     for i in range(len(coefficients) + 1):                  result = result + coefficients[i] * x ** i          return result  # нормализует функцию, чтобы интеграл от этой функции на заданном интервале был равен 1 def normalized_function(x):          integration, error = quad(function, lower_limit, upper_limit)          result = (1 / integration) * function(x)          return result  # вычисляем функцию распределения вероятностей в точке x def CDF(x):          integration, error = quad(normalized_function, lower_limit, x)          return integration  # находим значение, соответствующее искомому процентилю def get_percentile(percentile):          result = brentq(         lambda x: CDF(x) - percentile, lower_limit, upper_limit)          return result  lower_limit = 1 upper_limit = 10  print(f'25-Процентиль: {get_percentile(0.25)}') # 4.148070830469738 print(f'50-Процентиль: {get_percentile(0.50)}') # 6.635723812105114 print(f'75-Процентиль: {get_percentile(0.75)}') # 8.16891812557873 

В итоге, с изначальных дискретных процентилей в 4, 7 и 9, мы получили более точные непрерывные варианты в 4.14, 6.63 и 8.16 благодаря использованию аппроксимации. Построим график, позволяющий визуально оценить различия изначальных процентилей и уточнённых.

# визуализируем полученные данные plt.hist(data, lw=2.5, edgecolor='#ffffff', color="#fbf9ff", bins=x_length) plt.xticks(np.arange(1, x_length + 1, step=1))  # график, отражающий линию аппроксимации plt.plot(polyline, coefficients(polyline), lw=2.5, color='#c2a4f4')  for percentile in [0.25, 0.50, 0.75]:     plt.axvline(         x=get_percentile(percentile), label=f'{percentile} Процентиль (точный)', lw=2.5, color='#ade6b6', linestyle='-')     plt.axvline(         x=data['values'].quantile(percentile), label=f'{percentile} Процентиль (не точный)', color='#ffdb4d', linestyle='--')  # выключаем сетку plt.grid(visible=False)  # вызываем легенду plt.legend()  # вызываем график plt.show() 

Заключение

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

1. Построение полинома: Мы научились определять, какой полином лучше всего подходит для описания данных, используя коэффициент детерминации R².

2. Нормировка полинома: Для того чтобы полином можно было использовать в качестве функции плотности вероятности (PDF), мы нормировали его, что гарантировало, что интеграл на заданном интервале равен 1.

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

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