Об оценке плотности распределений

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

PDF estimation logo

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

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

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

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

Для оценки "качества" получаемых плотностей распределений будем использовать расстояние Кульбака-Лейблера.

Вычислительный эксперимент выполним на базе языка программирования Python и библиотеки научных вычислений SciPy.

Генерация экспериментальных данных

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

В качестве первой (простой) выборки будем рассматривать набор значений, полученный из нормального распределения с параметрами \((m=1, \sigma^2=4)\); примером более сложного распределения может быть объединение выборок из двух смещенных друг относительно друга нормальных  с параметрами \((m=10, \sigma^2=1)\)  и \((m=1, \sigma^2=4)\); 

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

Полученное в результате объединения двух нормальных выборок распределение является бимодальным (т.е.  плотность распределения имеет два локальных максимума); гистограмма 2000 элементной выборки, полученной из такого закона распределения приводится на рисунке слева.

Если \(f_1(x)\) и \(f_2(x)\) являются функциями плотностей распределений, из которых получены первая и вторая выборки, участвующие в объединении, то плотность распределения, которой соответствует полученная выборка, будет выражаться как среднее плотностей \(f_1(x)\) и \(f_2(x)\).

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

Графики теоретических плотностей распределения

В случае параметрического оценивания плотности распределения, необходимо задаться его видом. На практике, результаты измерений включают погрешности, обусловленные суммой большого числа каких-либо малых случайных воздействий. Вследствие этого, на основе центральной предельной теоремы теории вероятностей, вполне правомерно ожидать, что такие погрешности будут представлять выборку из нормального закона распределения. Плотность нормального распределения имеет следующий вид: $$
f(x) = \frac{1}{\sigma\sqrt{2\pi}}\exp^{-\frac{(x-m)^2}{2\sigma^2}},
$$

где \(\sigma^2\) и \(m\) — параметры — дисперсия и математическое ожидание соответственно.

Таким образом, применение параметрического подхода оценки плотности по выборочным данным предполагает отыскание оценок параметров \(\sigma^2\) и \(m\).

 

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
'''
Parametric and Nonparametric probability density estimations.


Author: Dmitry E. Kislov, kislov@easydan.com
Date: 11 Sept. 2015
'''

from __future__ import print_function

#It is assumed that numpy is imported as np 
from pylab import *
from sklearn.neighbors import KernelDensity
import scipy.stats as st


#Set large font size for axes ticks
rcParams.update({'font.size': 21})

#-----------------Creating samples for experiments-----------------------

#Simple normal distribution parameters
sigma1 = 4 
mean1 = 1

sigma2 = 1
mean2 = 10

#Sample from normal distribution
simple_sample = sigma1*np.random.randn(1000) + mean1

#Sample from complex bimodal distribution
complex_sample =np.array((sigma1*np.random.randn(1000) + mean1).tolist() + (sigma2*np.random.randn(1000)+mean2).tolist())


#Interval of interest
a,b = min(complex_sample), max(complex_sample)
domain = np.linspace(a, b, 1000)


#Exact pdf-function
exact_pdf1 = st.norm(mean1, sigma1)
pdf1_vals = exact_pdf1.pdf(domain)

exact_pdf2_aux = st.norm(mean2, sigma2)
pdf2_vals = 0.5*(exact_pdf2_aux.pdf(domain) + pdf1_vals)


#Plot theoretical pdf's
figure()
plot(domain, pdf1_vals, 'r', domain, pdf2_vals, 'b')


#--------------------Parametric estimation-------------------------------

#Parametric estimation for simple (normal) sample
#To exclude shift in estimation of sigma one should use ddof=1
sigma_est1 = np.std(simple_sample, ddof=1)
mean_est1 = np.mean(simple_sample)


#Parametric estimation for complex (combined) sample
sigma_est2 = np.std(complex_sample, ddof=1)
mean_est2 = np.mean(complex_sample)

#Create instances of normal pdfs
est_pdf1 = st.norm(mean_est1, sigma_est1)
est_pdf2 = st.norm(mean_est2, sigma_est2)

#Values of pdf estimations
est_pdf1_vals = est_pdf1.pdf(domain)
est_pdf2_vals = est_pdf2.pdf(domain)


#Comparison with exact pdfs
fig = figure()
fig.suptitle('Comaprison of parametric pdf estimations')
ax = subplot('121')
ax.plot(domain, est_pdf1_vals, 'r', domain, pdf1_vals, 'b')
ax.set_title('Simple dist.: red - est., blue - exact')

ax = subplot('122')
ax.plot(domain, est_pdf2_vals, 'r', domain, pdf2_vals, 'b')
ax.set_title('Complex dist.: red - est., blue - exact')

print('Parametric estimations comparison:')
print('Kullback-Leibler divergence d(pdf1,est_pdf1)=%s'%st.entropy(pdf1_vals, est_pdf1_vals))
print('Kullback-Leibler divergence d(pdf2,est_pdf2)=%s'%st.entropy(pdf2_vals, est_pdf2_vals))
#------------------------------------------------------------------------



#--------------------Nonparametric estimation----------------------------

#Scott rule implementation for bandwidth estimation 
#See details in D.W. Scott, Multivariate Density Estimation: 
#Theory, Practice, and Visualization, John Wiley&Sons, New York, Chicester, 1992.
dimensions = 1
bandwidth1 = float(len(simple_sample))**(-1./(dimensions+4))
bandwidth2 = float(len(complex_sample))**(-1./(dimensions+4))
#Previously Scott's rule have been tested, but constant bandwidth lead to more accurate results.

kde_pdf1_est = KernelDensity(kernel='gaussian', bandwidth=0.8).fit(simple_sample[:, np.newaxis])
kde_pdf2_est = KernelDensity(kernel='gaussian', bandwidth=0.8).fit(complex_sample[:, np.newaxis])

#One should use np.exp, see. documentation of kde-procedure in scikit-learn!
kde_pdf1_vals = np.exp(kde_pdf1_est.score_samples(domain[:,np.newaxis]))
kde_pdf2_vals = np.exp(kde_pdf2_est.score_samples(domain[:,np.newaxis]))


fig_kde = figure()
fig_kde.suptitle('Comaprison of nonparametric pdf estimations')
ax = subplot('121')
ax.plot(domain, kde_pdf1_vals, 'r', domain, pdf1_vals, 'b')
ax.set_title('Simple dist.: red - est., blue - exact')

ax = subplot('122')
ax.plot(domain, kde_pdf2_vals, 'r', domain, pdf2_vals, 'b')
ax.set_title('Complex dist.: red - est., blue - exact')

print('Parametric estimations comparison:')
print('Kullback-Leibler divergence d(pdf1,est_pdf1)=%s'%st.entropy(pdf1_vals, kde_pdf1_vals))
print('Kullback-Leibler divergence d(pdf2,est_pdf2)=%s'%st.entropy(pdf2_vals, kde_pdf2_vals))

#------------------------------------------------------------------------

#Plot all figures
show()
Результаты сравнения оценок (параметрический случай) Результаты сравнения оценок (непараметрический случай)

Рассчеты значений дистанции Кульбака-Лейблера для оценок и точных распределений вычисляются по ходу выполнения программы. Они приводятся ниже, но могут несколько отличаться по значениям, от тех, которые получатся если запустить код снова. Это связано с тем, что генерация осуществляется при помощи функции numpy.random.randn, которая по умолчанию создает при каждом новом запуске разные наборы псевдослучайных чисел (фиксировать значение seed в генераторе псевдослучайных чисел было принято нецелесообразным).
 

Parametric estimations comparison:
Kullback-Leibler divergence d(pdf1,est_pdf1)=0.00118575998792
Kullback-Leibler divergence d(pdf2,est_pdf2)=0.31338796242
Nonparametric estimations comparison:
Kullback-Leibler divergence d(pdf1,est_pdf1)=0.00331357721837
Kullback-Leibler divergence d(pdf2,est_pdf2)=0.0226708286761

Таким образом, можно видеть, что параметрическая оценка плотности в случае соответствия выборочных данных гипотетическому распределению, является весьма точной (расстояние Кульбака-Лейблера (d) в этом случае является минимальным (d = 0.001185)). Результаты непараметрического оценивания в этом случае оказываются несколько хуже (d = 0.00331).

Если данные подходы используются для оценки плотности более сложного (существенного отличного от нормального) распределения, то параметрический подход дает весьма плохие результаты (d = 0.313387). Использование непараметрического подхода — ядерной оценки плотности с гауссовым ядром и фиксированным размером скользящего окна, позволяет получить гораздо лучшие результаты (d = 0.022670).

 

blog comments powered by Disqus