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

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

Для решения задачи будем использовать Python и сторонний модуль Matplotlib, который поможет легко и наглядно визуализировать полученные результаты и удостовериться в их адекватности. Из Pandas возьмем только DataFrame для удобства.

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

В данном случае предлагаю не забивать голову формулами, а сосредоточиться на понимании сути вопроса и принципа его решения. Если уловить смысл, то несложные формулы вы совершенно легко сможете написать самостоятельно. А вот код на Python мы напишем, будет интереснее математических формул (заодно и вспомним принципы работы с Matplotlib).

Итак, опишем процесс простым и понятным языком. Для решения вопроса нам нужно пройти через следующие этапы:

• определяем количество элементов в рассматриваемой выборке;

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

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

• количество вхождений в каждый интервал делим на произведение ширины интервала на количество чисел в нашей выборке, таким образом получаем значения кривой плотности распределения вероятности по оси ординат;

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

  Проверяем корректность описанного способа следующим образом:

• генерируем четыре выборки случайных чисел для четырех видов распределений соответственно: Релея, гамма, Вейбулла и экспоненциального;

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

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

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

Напишем код на Python. Все должно быть понятно, поскольку есть комментарии (они существенно увеличили количество строк, но с ними все гораздо яснее).

import random
import math
import matplotlib.pyplot as plt
import pandas
import pandas as pd


def rel_pdf(rel_sigma: float, n: int, delta_x: float) -> pandas.DataFrame:
    """
    Вычисляет Релеевскую кривую плотности распределения вероятности по известной формуле
    :param rel_sigma: среднеквадратическое отклонение
    :param n: количество рассчитанных точек
    :param delta_x: шаг расчета по оси Ох
    :return: pandas.DataFrame
    """
    pdf_x = []  # Координаты по оси абсцисс
    pdf_y = []  # Координаты по оси ординат
    for i in range(1, n):
        x = i * delta_x
        pdf_x.append(x)
        pdf_y.append(((2 * x) / rel_sigma) * math.exp((-x ** 2) / rel_sigma))
    return pd.DataFrame({'x': pdf_x, 'y': pdf_y})


def rel_rand(n: int, rel_sigma: float) -> list:
    """
    Генерирует случайные числа с Релеевской плотностью распределения вероятности
    :param n: количество отсчетов
    :param rel_sigma: среднеквадратическое отклонение
    :return: list
    """
    rel_list = []
    for i in range(n):
        rel_list.append((rel_sigma / math.sqrt(2)) * math.sqrt(-2 * math.log(random.uniform(0, 1))))
    return rel_list


def gam_pdf(v: float, b: float, n: int, delta_x: float) -> pandas.DataFrame:
    """
    Вычисляет кривую плотности Гамма распределения вероятности по известной формуле
    :param v: параметр формы
    :param b: масштабный коэффициент
    :param n: количество рассчитанных точех
    :param delta_x: шаг расчета по оси Ох
    :return: pandas.DataFrame
    """
    pdf_x = []  # Координаты по оси абсцисс
    pdf_y = []  # Координаты по оси ординат
    for i in range(1, n):
        x = i * delta_x
        pdf_x.append(x)
        pdf_y.append(((b ** v) / math.gamma(v) * (x ** (v - 1)) * math.exp(-b * x)))
    return pd.DataFrame({'x': pdf_x, 'y': pdf_y})


def gam_rand(n: int, v: float, b: float) -> list:
    """
    Генерирует случайные числа с гамма распределением плотности вероятности
    :param n: количество отсчетов
    :param v: параметр формы
    :param b: масштабный коэффициент
    :return: list
    """
    gam_list = [random.gammavariate(v, 1 / b) for i in range(n)]
    return gam_list


def weib_pdf(a: float, b: float, n: int, delta_x: float) -> pandas.DataFrame:
    """
    Вычисляет кривую плотности распределения вероятности Вейбулла по известной формуле
    :param a: масштабный коэффициент
    :param b: параметр формы
    :param n: количество рассчитанных точех
    :param delta_x: шаг расчета по оси Ох
    :return: pandas.DataFrame
    """
    pdf_x = []  # Координаты по оси абсцисс
    pdf_y = []  # Координаты по оси ординат
    for i in range(1, n):
        x = i * delta_x
        pdf_x.append(x)
        pdf_y.append((b / a) * ((x / a) ** (b - 1)) * math.exp(-(x / a) ** b))
    return pd.DataFrame({'x': pdf_x, 'y': pdf_y})


def weib_rand(n: int, a: float, b: float) -> list:
    """
    Генерирует случайные числа с распределением плотности вероятности Вейбулла
    :param n: количество отсчетов
    :param a: масштабный коэффициент
    :param b: параметр формы
    :return: list
    """
    wei_list = [random.weibullvariate(a, b) for i in range(n)]
    return wei_list


def exp_pdf(l: float, n: int, delta_x: float) -> pandas.DataFrame:
    """
    Вычисляет кривую экспоненциальной плотности распределения вероятности по известной формуле
    :param l: обратный коэффициент масштаба
    :param n: количество рассчитанных точех
    :param delta_x: шаг расчета по оси Ох
    :return: pandas.DataFrame
    """
    pdf_x = []  # Координаты по оси абсцисс
    pdf_y = []  # Координаты по оси ординат
    for i in range(1, n):
        x = i * delta_x
        pdf_x.append(x)
        pdf_y.append(l * math.exp(-l * x))
    return pd.DataFrame({'x': pdf_x, 'y': pdf_y})


def exp_rand(n: int, l: float) -> list:
    """
    Генерирует случайные числа с экспоненциальным распределением плотности вероятности
    :param n: количество отсчетов
    :param l: обратный коэффициент масштаба
    :return: list
    """
    exp_list = [random.expovariate(l) for i in range(n)]
    return exp_list


def pdf(k: int, rnd_list: list) -> pandas.DataFrame:
    """
    Получает кривую плотности распределения вероятности
    :param k: количечиво интервалов разбиения гистограммы
    :param rnd_list: случайный процесс
    :return: pandas.DataFrame
    """
    pdf_x = []  # Координаты по оси абсцисс
    pdf_y = []  # Координаты по оси ординат
    n = len(rnd_list)  # количество элементов в рассматриваемой выборке
    h = (max(rnd_list) - min(rnd_list)) / k  # ширина одного интервала
    a = min(rnd_list)  # минимальное значение в рассматриваемой выборке
    for i in range(0, k):  # проход по интервалам
        count = 0
        for j in rnd_list:  # подсчет количества вхождений значений из выборки в данный интервал
            if (a + i * h) < j < (a + (i * h) + h):
                count = count + 1
        pdf_x.append(a + i * h + h / 2)  # координата по оси абсцисс полученной кривой плотности распределения
        # вероятности
        pdf_y.append(count / (n * h))  # координата по оси ординат полученной кривой плотности распределения
        # вероятности
    d = {'x': pdf_x, 'y': pdf_y}
    return pd.DataFrame(d)


rrand = rel_rand(100000, 1)  # генерируем случайные числа с распределением Релея
r_pdf = pdf(100, rrand)  # строим плотность распределения вероятности по полученной выборке
rel = rel_pdf(1, 100, 0.05)  # получаем кривую распределения Релея по известной формуле

grand = gam_rand(100000, 0.5, 0.5)  # генерируем случайные числа с гамма распределением
g_pdf = pdf(400, grand)  # строим плотность распределения вероятности по полученной выборке
gam = gam_pdf(0.5, 0.5, 500, 0.01)  # получаем кривую гамма распределения по известной формуле

wrand = weib_rand(100000, 1, 5)  # генерируем случайные числа с распределением Вейбулла
w_pdf = pdf(100, wrand)  # строим плотность распределения вероятности по полученной выборке
weib = weib_pdf(1, 5, 500, 0.01)  # получаем кривую распределения Вейбулла по известной формуле

exprand = exp_rand(100000, 1.5)  # генерируем случайные числа с экспоненциальным распределением
e_pdf = pdf(100, exprand)  # строим плотность распределения вероятности по полученной выборке
exp = exp_pdf(1.5, 500, 0.01)  # получаем кривую экспоненциального распределения по известной формуле

# Построение графиков
fig, ax = plt.subplots(nrows=2, ncols=2)  # Создаем фигуру и четыре системы координат на фигуре

ax[0, 0].plot(r_pdf['x'], r_pdf['y'], 'r-', label='Оценка PDF')  # Своя оценка плотности распределения
ax[0, 0].plot(rel['x'], rel['y'], 'b-', label='PDF')  # Плотность распределения построенная по формуле
ax[0, 0].set(title='Релеевское распределение')  # Заголовок для системы координат
ax[0, 0].set_xlabel('x')  # Ось абсцисс
ax[0, 0].set_ylabel('PDF(x)')  # Ось ординат
ax[0, 0].set_xlim(xmin=0, xmax=4)  # Минимальное и максимальное значение по оси Х
ax[0, 0].legend()  # Отображаем легенду для данной системы координат
ax[0, 0].grid()  # Отображаем сетку на системе координат

ax[0, 1].plot(g_pdf['x'], g_pdf['y'], 'r-', label='Оценка PDF')  # Своя оценка плотности распределения
ax[0, 1].plot(gam['x'], gam['y'], 'b-', label='PDF')  # Плотность распределения построенная по формуле
ax[0, 1].set(title='Гамма распределение')  # Заголовок для системы координат
ax[0, 1].set_xlabel('x')  # Подписи по оси абсцисс
ax[0, 1].set_ylabel('PDF(x)')  # Подписи по оси ординат
ax[0, 1].set_xlim(xmin=0, xmax=4)  # Минимальное и максимальное значение по оси Х
ax[0, 1].legend()  # Отображаем легенду для данной системы координат
ax[0, 1].grid()  # включение отображение сетки

ax[1, 0].plot(w_pdf['x'], w_pdf['y'], 'r-', label='Оценка PDF')
ax[1, 0].plot(weib['x'], weib['y'], 'b-', label='PDF')
ax[1, 0].set(title='Распределение Вейбулла')
ax[1, 0].set_xlabel('x')  # ось абсцисс
ax[1, 0].set_ylabel('PDF(x)')  # ось ординат
ax[1, 0].set_xlim(xmin=0, xmax=4)
ax[1, 0].legend()
ax[1, 0].grid()  # включение отображение сетки

ax[1, 1].plot(e_pdf['x'], e_pdf['y'], 'r-', label='Оценка PDF')
ax[1, 1].plot(exp['x'], exp['y'], 'b-', label='PDF')
ax[1, 1].set(title='Экспоненциальное распределение')
ax[1, 1].set_xlabel('x')  # ось абсцисс
ax[1, 1].set_ylabel('PDF(x)')  # ось ординат
ax[1, 1].set_xlim(xmin=0, xmax=4)
ax[1, 1].legend()
ax[1, 1].grid()  # включение отображение сетки

plt.show()

В результате выполнения кода имеем следующие кривые (визуализированы с помощью Matplotlib).

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

Ссылка на репозиторий.