Распределения вероятностей в SciPy

генерирование выборок из различных распределений на Python

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

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

Класс для работы с вероятностными раcпределениями в scipy.stats имеет следующие методы:

Методы Назначения
pdf Плотность распределения вероятностей
cdf Функция распределения вероятностей F(x)
sf 1 - F(x) [т.е. функция "хвоста" распределения]
ppf Обратная функция распределения вероятностей
(часто используется при проверке статистических гипотез)
isf обратная функция к 1-F(x)
stats Базовые характеристики: среднее, дисперсия, асимметрия, эксцесс
moment Нецентрованные моменты k-го порядка распределения
rvs Генерация случайных выборок

Структура классов распределений в Scipy

Представленные в модуле stats  распределения являются объектами двух базовых классов распределений — это rv_continuous and rv_discrete — для непрерывных и дискретных случаев соответственно.

Чтобы создать свое собственное непрерывное распределение достаточно создать класс-потомок от rv_continuous и переопределить хотя бы один из двух методов _pdf или _cdf, которые отвечают соответственно функции плотности распределения и функции распределения вероятностей.

Рассмотрим процесс определения собственного распределения на примере. Определим класс Uniform:

#coding: utf-8
#author: Dmitry E. Kislov
#e-mail: kislov@easydan.com

import scipy.stats as st

class Uniform(st.rv_continuous):
    
    # Плотность распределения вероятностей для равномерного
    # распределения на интервале [0,1]
    def _pdf(self, x):
        return 1.0 if 0.0 <= x <= 1.0 else 0.0

    # Нет необходимости вводить параметры a и b, чтобы получить плотности
    # на произвольном интервале [a,b], rv_continuous содержит необходимый
    # функционал, чтобы создавать плотности вида f((x-loc)/scale), где
    # loc, scale -- параметры смещения и масштаба, это обеспечивает возможность
    # создавать в данном случае плотности на произвольных интервалах [a,b],
    # подбирая соответствующие параметры loc и scale 

Через объекты класса Uniform мы можем получить доступ к распределению и методам (pdf, cdf и т.д.), которые определены за счет наследования от rv_continuous.

Важно отличать объекты распределений с фиксированными (frozen) параметрами и неопределенными. Рассмотрим объект uniform_arbitrary — представитель класса Uniform, с неопределенными параметрами распределения, что означает, возможность установки параметров масштаба и сдвига при обращении к функциям pdf, cdf и т.п.

uniform_arbitrary = Uniform()
print('Распределение с произвольными параметрами (интервал распределения:[0, 100]):',
      uniform.pdf(0.5, loc=0.0, scale=100))
print('Распределение с произвольными параметрами (интервал распределения:[0, 50]):',
      uniform.pdf(0.5, loc=0.0, scale=50))

# Равномерное распределение на интервале [0, 100],
# При вызове .pdf(x) будет вычисляться ._pdf((x-loc)/scale),
# что в данном случае означает, что pdf(x) -- функция плотности для равномерного
# распределения на интервале [0, scale]
uniform_frozen = Uniform(loc=0.0, scale=100)

В случае распределения с фиксированными параметрами масштаба и сдвига uniform_frozen, переопределить их при вызове функции плотности (pdf) и\или функции распределения вероятностей (cdf) не получится:

uniform = Uniform()
uniform_frozen = uniform(loc=0.0, scale=10)
uniform_frozen.pdf(0.0, loc=2) # Попытка переопределить параметры
TypeError: __init__() got an unexpected keyword argument 'loc'

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

При создании собственных распределений важно следить, за тем, что возращают переопределяемые функции _pdf и\или _cdf. Легко можно определить несуществующее распределение (например, не удовлетворяющее условию нормировки), задав, например, в качестве функции плотности распределения,  линейную функцию f(x) = x без каких-либо дополнительных условий. Такой функции плотности, очевидно, быть не может (ввиду того, что она хотя бы может принимать отрицательные значения),  однако, построить такое неверное распределение, переопределив _pdf вполне даже удастся, и, возможно, получится даже сгенерировать случайные величины ему "соответствующие", вызвав метод rvs. Но только будет нельзя доверять результатам.

Если мы хотим определить распределение с дополнительным параметром, то можно поступить следующим образом:

class Uniform(st.rv_continuous):
    
    def __init__(self, *args,  **kwargs):
        self._my_a = 1 if 'par_a' not in kwargs else kwargs['par_a']
        kwargs.pop('par_a')
        return super(Uniform, self).__init__(*args, **kwargs)
        
    
    # Плотность распределения вероятностей для равномерного
    # распределения на интервале [0, par_a], по-умолчанию par_a=1
    def _pdf(self, x):
        return 1.0/self._my_a if 0 <= x <= self._my_a else 0.0
    
uniform_arbitrary = Uniform(par_a=10)

print('Случайное число из интервала [0, 10]:', uniform_arbitrary.rvs())  

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

Приведем примеры генерирования выборок заданного объема из некоторых распределений: нормального, Коши, равномерного и Парето.

import scipy.stats as st

# --------- Равномерное распределение ------------
# Равномерное распределение на интервале [-1, 1]
uniform = st.uniform(loc=-1, scale=2)

# Генерация 100 случайных чисел
uniform.rvs(size=100)
# ------------------------------------------------

# --------- Нормальное распределение ------------
# Нормальное распределение с параметрами mean=5, sigma=2
norm = st.norm(loc=5, scale=2)

# Генерация 100 случайных чисел
norm.rvs(size=100)
# ------------------------------------------------


# --------- Распределение Парето------------
# функция плотности pdf(x, b) = b / x**(b+1),
# пример для случая b=2
pareto = st.pareto(2, loc=0, scale=1)

# Генерация 100 случайных чисел
pareto.rvs(size=100)
# ------------------------------------------------


# --------- Распределение Коши------------
# функция плотности pdf(x) = 1 / (pi * (1 + x**2))
cauchy = st.cauchy()

# Генерация 100 случайных чисел
cauchy.rvs(size=100)
# ------------------------------------------------

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

blog comments powered by Disqus