close

Вход

Забыли?

вход по аккаунту

?

Vorobev 0D950131DA

код для вставкиСкачать
Федеральное агенТство по образованию
Государственное образовательное учреждение
высшего профессионального образования
Санкт-петербургский государственный университет
аэрокосмического приборостроения
С. Н. Воробьев, Н. В. Гирина, Л. А. Осипов
Гистограммное оценивание
плотности распределения
Учебное пособие
Допущено Учебно-методическим объединением вузов
по университетскому политехническому образованию
в качестве учебного пособия для студентов
высших учебных заведений, обучающихся по направлению
подготовки 230200-Информационные системы.
Санкт-Петербург
2008
1
УДК 004.93
ББК 32.93
В75
Рецензенты:
кафедра информационных управляющих систем Санкт-Петербургского
государственного университета телекоммуникаций;
доктор технических наук, профессор С. А. Яковлев
Утверждено
редакционно-издательским советом университета
в качестве учебного пособия
Воробьев С. Н., Гирина Н. В., Осипов Л. А.
В75 Гистограммное оценивание плотности распределения:
учебное пособие / С. Н. Воробьев, Н. В. Гирина, Л. А. Осипов. – СПб.: ГУАП, 2008. – 92 с.: ил.
ISBN 978-5-8088-0388-6
Рассматриваются ядерные и гистограммные методы оценивания
плотности распределения выборки. Приведены примеры оценивания, подтверждающие работоспособность гистограммных методов
при ограниченных размерах обучающей выборки. Метод нормальной
аппроксимации распространяется на двумерные плотности распределения.
Предназначено для студентов, обучающихся по направлению
230200 «Информационные системы».
УДК 004.93
ББК 32.93
ISBN 978-5-8088-0388-6
2
© ГУАП, 2008
© С. Н. Воробьев,
Н. В. Гирина,
Л. А. Осипов, 2008
Содержание
Введение. ......................................................................... 1. Критерии согласия........................................................ 1.1. Гистограмма...................................................... 1.2. Критерий А. Н. Колмогорова................................ 1.3. Критерий χ2 К. Пирсона....................................... 1.4. Проверка гипотезы однородности.......................... 1.5. Критерий ω2....................................................... 2. Ядерные оценки плотности распределения........................ 2.1. Метод парзеновских окон..................................... 2.2. Метод n ближайших соседей................................. 3. Гистограммное оценивание............................................. 3.1. Нормальная гистограммная аппроксимация .......... 3.2. Нормальная гистограммная интерполяция............ 3.3. Адаптивная нормальная МНК-интерполяция......... 3.4. Моделирование корреляционного обнаружения...... 4. Оценивание двумерных плотностей.................................. 4.1. Двумерная гистограмма....................................... 4.2. Нормальная двумерная гистограммная
аппроксимация ........................................................ Заключение..................................................................... Библиографический список................................................ 4
11
11
13
18
21
25
27
27
33
39
39
50
58
60
72
72
78
90
91
3
Введение
Государственный образовательный стандарт высшего профессионального образования по направлению 230200 «Информационные
системы» предусматривает владение методами статистического
моделирования на ЭВМ с использованием имитационных моделей
информационных процессов. В дисциплине «Моделирование информационных систем» изучаются приемы анализа и интерпретации результатов моделирования на ЭВМ, оценки точности и достоверности результатов моделирования.
Одна из задач оценивания и интерпретации результатов моделирования – оценка закона распределения случайной составляющей наблюдений или результатов их обработки в информационной
системе. В адаптивные системы включается канал слежения за помехами: измерение вероятностных свойств нестационарного шума
позволит, например, оперативно перестраивать систему передачи
данных.
Классическая задача математической статистики – проверка
(x) = F(x) о равенстве (близости) оценки F
(x) фунгипотезы H: F
кции распределения наблюдений X заданной функции распределения F(x). Гипотеза H – предположение о том, что выборка X = {x1,
..., xN} (наблюдения) принадлежит генеральной совокупности с
функцией распределения F(x) – проверяется критериями согласия
[1 –3]. Критерий согласия А. Н. Колмогорова базируется на оценке
(x) в виде эмпирической функции распределения F (x), построF
N
енной по выборке X. Критерий согласия χ2 К. Пирсона базируется
на гистограмме H – последовательности значений
hi =
ni
N∆
(1)
относительного числа попаданий выборочных значений xj в k заданных интервалов числовой оси шириной ∆ каждый (j = 1, ..., N).
Теоретическая основа критерия А. Н. Колмогорова – сходимость
по вероятности эмпирической функции распределения к функции
распределения генеральной совокупности


lim N →∞ p  FN (x) − F(x) < ε



 = 1.

Сходимость относительного числа попаданий ni/N к теоретической
вероятности
4
pi =
∫ f (x)dx
∆i
попадания в i-й интервал гистограммы


lim N →∞ p  hi − pi < ε



 = 1, i = 1, ..., k,

делает гистограмму хорошим статистическим аналогом плотности
распределения f(x) генеральной совокупности.
Критерии согласия подробно обсуждаются в монографии [4], где
эта задача называется задачей проверки согласия.
Гипотезы согласия классифицируются как непараметрические,
то есть не требующие знания параметров проверяемых законов распределения, например моментов нормального распределения F(x).
Это тем более относится к сложным гипотезам, таким как «проверка нормальности» [4] – генеральная совокупность распределена по
нормальному закону с произвольными моментами. Стандартные
непараметрические гипотезы согласия: гипотеза однородности, гипотеза независимости, гипотеза случайности, гипотезы о сдвиге и
масштабе [3 –5]. Критерии Н. В. Смирнова и хи-квадрат – основные
критерии проверки однородности выборок X1 и X2 (или нескольких выборок).
Задача оценивания плотности распределения f(x) в общем случае также непараметрическая – предполагается, что неизвестны
не только параметры распределения, но и о самом распределении
нельзя сформулировать конкретную гипотезу. Тогда оценивание
закона распределения случайной составляющей информационного
сигнала выходит за рамки проверки гипотез и становится самостоятельной задачей. Такая ситуация может возникнуть в системах со
сложными комбинациями нелинейных преобразований сигнала,
например при непараметрическом обучении в распознавании образов [6,7].
Как и любое решение в математической статистике, оценка
плотности – некоторая приемлемая статистика
f (x) = ϕ(X).
Тот или иной алгоритм j можно проверить, оценивая известные
плотности распределения. Проверка алгоритма есть выяснение
близости полученной оценки f (x) к оцениваемой плотности f(x) – в
5
этом состоит смысл применения критериев согласия (разд.1) в данном пособии.
В прикладной статистике задача оценки плотности распределения f(x) эргодического процесса x(t) решается на основании вероятности p∆ попадания мгновенного значения процесса в интервал
∆ [8]


p∆ = p  x ∈ ∆



 = f (x)dx.
 ∆
∫
Вычисляется время T∆ пребывания процесса в интервале ∆, и оценка вероятности 
p ∆ = T∆ /T дает оценку плотности f (x) в точке x =
= ∆/2
T
∆
f  x =  = ∆ , 
2  T∆
(2)
T – длительность реализации. Рис. 1 иллюстрирует эту процедуру
для случая ∆ = (0,5;1,0), T = 50, T = 9, так что оценка равна 
p =
∆1
1
∆1
= f (x = 0,75) = 0,18. Здесь время измеряется в отсчетах дискретизированного процесса X. Номера j отсчетов xj, попавших в интервал ∆1, и их значения – в табл. 1.
x
2
1.5
1
delta 1
0.5
0
T
-0.5
-1
-1.5
-2
-2.5
0
5
10
15
20
25
30
35
40
Рис. 1. К оценке плотности распределения
6
45
50
Таблица 1
j
xj
1
3
17
29
32
34
36
38
49
0.635 0.551 0.944 0.587 0.668 0.889 0.525 0.913 0.680
В интервал ∆1 = 0,5 попадает n∆ = 9 значений процесса X из пятидесяти.
Если весь диапазон значений процесса разбит на интервалы ∆
(рис. 1), относительное число попаданий ni/N = Ti/T процесса в
каждый интервал дает гистограмму (1), а оценка (2), записанная
в виде
k
1
f (x) =
ni , N∆ i =1
∑
(3)
есть гистограммная оценка плотности распределения.
Другой подход к оцениванию плотности распределения – построение ядерной оценки [3]
f (x) =
1
Na N
N
 x − xj 
, aN 
∑ ϕ
j =1
(4)
которая при соответствующем выборе ядра j(x) и коэффициента aN
сходится по вероятности к непрерывной плотности распределения
генеральной совокупности f(x):


P  lim N →∞ sup x <∞ f (x) − f (x) = 0



 = 1.

Ядерные оценки применяются в распознавании образов – оценки методом Парзена и методом n ближайших соседей [6,7].
Оценка Парзена имеет вид (4), ядро j(x;xj) называется окном и
является сглаживающей функцией со свойствами плотности распределения. Если ядро – нормальная плотность, оценка Парзена
2
N
 (x − x j ) 
1
1
f (x) =
exp −

N j =1 2πλ N
2λ 2N 

∑
(5)
есть усредненная сумма нормальных плотностей с модами в случайных точках xj. Коэффициент aN задается зависимым от размера
выборки: λ N = L / N . Число L подбирается: при больших значениях λ суммируются медленно меняющиеся функции и получается сильно сглаженная «несфокусированная» оценка; при малых λ
7
складываются узкие слабо пересекающиеся выбросы с центрами в
точках xj – оценка плотности «зашумляется».
В оценке методом n ближайших соседей используются сглаживающие функции с переменной шириной. Ширина задается не функцией размера выборки, а функцией наблюдений – такой, чтобы
захватывались k соседних выборочных значений xj. Метод n ближайших соседей позволяет отслеживать характер плотности выборки X. Вблизи моды плотности f (x) значения xj выпадают более
компактно, чем на «хвостах», поэтому вблизи моды разрешение
оценки выше. Если плотность многомодальна, то ширина сглаживающих функций уменьшается вблизи каждой моды [6].
Недостаток ядерных оценок – необходимость их аппроксимации, так как численную оценку в виде суммы N слагаемых трудно
использовать в аналитических расчетах.
Оценки Парзена и n ближайших соседей трактуются как развивающие гистограммный метод [6,7], тогда как в математической
статистике ядерные оценки самостоятельны [3]. Действительно,
суть понятия гистограммы состоит в группировке N выборочных
значений xj в k<<N интервалах и замене выборки k числами – значениями середин интервалов (случайными или неслучайными в
зависимости от расстановки интервалов на числовой оси). Произвольные расширения этого понятия (назначение k = N, попадание
одного и того же значения xj в различные интервалы) вызывают
сомнения в правомерности отнесения ядерных методов к гистограммным. Оценки Парзена и ближайших соседей включены в пособие (разд. 2) для сравнения с гистограммными.
Гистограммными методами оценивания плотностей в данном
пособии считаются методы, базирующиеся на классическом определении гистограммы. Гистограммная оценка плотности – численная: непрерывная функция f(x) заменяется числами hi = ni/N∆,
i = 1, ..., k. Естественное улучшение численной оценки – ее аппроксимация непрерывной функцией
f (x) =
k
n
∑ Ni f0 (x;xi ), i =1
в которой используются нормированные базисные функции
∞
∫ f0 (x;xi )dx = 1,
−∞
8
(6)
привязанные к серединам xj интервалов гистограммы. Формально
оценка (6) отличается от ядерной оценки (4) количеством слагаемых функционального разложения. По существу оценка (6) сглаживает не исходные выборочные значения xj, как (4), а наблюдения, сгруппированные гистограммой, поэтому следует ожидать ее
большей эффективности. Оценка (6) с нормальными базисными
функциями называется нормальной гистограммной аппроксимацией (разд. 3).
Другой способ получения непрерывной оценки f (x) сглаживанием гистограммы – применение метода наименьших квадратов
(МНК) [3,9]. Нормальные базисные функции записываются в виде
симметричной матрицы плана Z, МНК – коэффициенты b рассчитываются по вектору H значений hi = ni/N гистограммы:
b = Z −1H.
В узлах xj значения МНК-оценки плотности
f (x ) = Z T b = Zb = ZZ −1H = H
i
оказываются равными значениям гистограммы, то есть МНК-функция f (x)  – интерполирующая. МНК-оценка плотности распределения называется нормальной МНК-интерполяцией (разд. 3).
Нормальная аппроксимация и МНК-интерполяция гистограмм,
соответствующих гладким плотностям f(x), дают удовлетворительные по критериям согласия результаты при размерах выборки в десятки-сотни наблюдений. Суммирование относительно небольшого числа k базисных функций позволяет записать гистограммные
оценки в виде относительно компактных формул.
Применение гистограммных оценок плотностей распределения
иллюстрируется подробным примером моделирования нелинейного обнаружителя импульсного сигнала (подразд. 3.4). Гистограммные оценки статистики обнаружения позволяют получить достаточно гладкие рабочие характеристики при ограниченном размере
выборки и не требуют чрезмерного уменьшения гистограммных
интервалов. Использование в этом примере оценок Парзена также
приводит к приемлемым результатам, что подтверждает сделанный ранее вывод о том, что основной недостаток ядерных оценок
связан с их аппроксимацией.
Нормальная гистограммная аппроксимация может быть распространена на многомерные выборки, однако компъютерная обработка многомерных массивов может потребовать большого времени.
В пособии рассмотрены двумерные гистограммы и гистограммное
9
оценивание двумерных плотностей распределения (разд. 4). Приведенные примеры показывают,что двумерная гистограммная аппроксимация может оказаться полезной при обработке изображений.
Основные теоретические положения иллюстрируются примерами, выполненными в системе MATLAB.
10
1. Критерии согласия
1.1. Гистограмма
Пусть выборка XT = [x1, ..., xN] извлечена из генеральной совокупности X с известной плотностью распределения f(x). Отрезок
числовой оси x = [xmin, xmax] разделяется на k интервалов шириной
∆ и подсчитывается ni – число попаданий выборочных значений xj в
каждый из интервалов:
x
− x min
∆ = max
.
k
Таков алгоритм MATLAB – функции HIST [10]: k + 1 граница xih
между интервалами и ширина интервалов ∆ случайны (конкретны
для каждой выборки X). Например, выборка размером N = 50, формируемая функцией RANDN, (точки на оси x) группируется в k = 6
интервалах (диапазон от xmin до xmax, ось xx) с числами попадания
в данной реализации n1 = 4, n2 = 6, n3 = 15, n4 = 14, n5 = 5, n6 = 6
(рис. 1.1).
xxx
4
6
15
5
14
6
x
del
xx
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
3
Рис. 1.1. Построение гистограммы
Независимо от значений ∆ относительное число попаданий в интервалы
HT =
1
[n1, ..., nk ], N
(1.1)
11
так что
k
1
ni = 1.
N i =1
∑
Если значение ni/N считать оценкой вероятности попадания в
интервал ∆i, то оценками значений плотности распределения f(x)
генеральной совокупности в точках x1 = (x1h+x2h)/2, …, xk = (xkh+
+ xk+1h)/2 (середины интервалов) становятся отношения hi = ni/N∆,
так что оценка плотности нормирована:
∞
k
1
f (x)dx =
ni ∆ = 1
N
∆ i =1
−∞
∫
∑
(численное интегрирование по формуле прямоугольников). Таким
образом, в зависимости от задачи гистограмма определяется либо
вектором (1.1), либо вектором
H T = [h1, ..., hk ] =
1
[n1, ..., nk ]. N∆
(1.2)
Гистограмма может строиться и с неслучайными интервалами.
Например, если есть основания считать, что генеральная совокупность описывается стандартным нормальным распределением
N(0,σ), естественно границы между интервалами задать равными
xkh = ±iσ/2, i = 1, 2, …, 6, и построить гистограмму с k = 12 интервалами. В некоторых случаях задаются интервалы с переменной
шириной [4].
Основное вероятностное свойство гистограммы следует из теоремы Бернулли о сходимости по вероятности отношения ni/N к истинной вероятности pi i-го события [3,11]:
n

ni p

→ pi , то есть lim n→∞  i − pi < ε  = 1.
N
N

(1.3)
В литературе [3,6,7] отмечаются многие теоретические недостатки
гистограммы, в частности сходимость при ∆ → 0 (при ∆ ≈ 0 и конечном размере N в ряд интервалов может не попасть ни одного выборочного значения xj, и гистограмма может потерять смысл). Тем
не менее отсутствие оснований отвергнуть гипотезу H: f (x) ≈ f(x)
по критерию χ2 Пирсона считается надежным аргументом в пользу
того, что проверяемая выборка извлечена из генеральной совокупности с плотностью распределения f(x).
12
На практике значение hi = ni/N гистограммы (1) при большом
размере выборки N и относительно небольшой ширине интервалов
∆ приближается к значению плотности f(xi). Гистограмму называют статистическим аналогом (оценкой) непрерывной плотности
распределения f(x) генеральной совокупности [3].
1.2. Критерий А. Н. Колмогорова
Выборка X упорядочивается по возрастанию
значений xj – запи
сывается в виде вариационного ряда X. Значениям x числовой оси
соответствует значение эмпирической функции распределения –
относительного числа выборочных значений, не превосходящих x:
FN (x) =
{j : x j ≤ x}
,
(1.4)
N
min(xj) < x < max(xj). Эмпирическая функция распределения – случайная функция, значения которой равны 0, 1/N, 2/N, ..., 1.
Пример 1.1. Вариационный ряд выборки размером N = 9 из генеральной совокупности со стандартным нормальным распределением (функция RANDN; здесь и далее используются функции
MATLAB со ссылками на литературу только в необходимых случаях) приведена в табл. 1.1.
5
6
7
8
9
0,6369
0,7442
1,2549
–1,3354
4
0,4615
xj
3
0,1481
2
–0,1095
1
–0,9184
j
–2,2162
Таблица 1.1
Эмпирическая функция распределения, возрастающая на величину 1/N = 1/9 в точках xj, показана на рис. 1.2. Функция FN (x) обладает всеми свойствами функции распределения: изменяется от
нуля до единицы, не убывает и непрерывна справа.
Эмпирическая функция распределения сходится по вероятности к функции распределения генеральной совокупности [1 –4]:
p
FN (x) 
→ F(x).
Для нее получены также более сильные результаты.
1. Максимальное отклонение эмпирической функции распределения от теоретической
13
FN
1
0.8
0.6
0.4
0.2
0
-2.5
-2
xmin
-1.5
-1
-0.5
0
0.5
1
xmax
1.5
Рис. 1.2. Эмпирическая функция распределения
D N = sup x <∞ FN (x) − F(x) (1.5)
с вероятностью единица сколь угодно мало при большом размере
выборки N на всей оси (теорема В. И. Гливенко, 1933 г.).
2. Если функция F(x) непрерывна, то при любом t > 0


lim N →∞ p 


∞
2 2

ND N ≤ t  = K(t) =
(−1) i e −2i t ,
i =−∞

∑
∞

t 
i −1 −2i 2t 2
lim N →∞ p  D N >
 = 1 − K(t) = 2 (−1) e
N

i =1
∑
(1.6)
(теорема А. Н. Колмогорова, 1933 г.).
3. Разность эмпирической и теоретической функций распределения асимптотически нормальна [12]:
F(x)(1 − F(x))
.
N
На основе теорем (1.5) и (1.6) разработан критерий согласия
Колмогорова: используется статистика (1.5); соответствие между
∆F(x) N →∞ ∈ Ν(0, σ x ), σ 2x =
14
эмпирической и гипотетической функциями распределения оценивается вероятностью

p

∞
−2i 2t02

ND N > t0  = 1 − K(t0 ) = 2 (−1) i −1 e

i =1
∑
того, что, если гипотеза верна, статистика t = ND N превзойдет
критический уровень t0. Функция K(x) =
∞
∑ (−1) ie −2i x
2 2
 – функция
i =−∞
распределения Колмогорова. Другая интерпретация критерия: вероятность
∞
∑ (−1) i−1e −2i t ,
p=2
2 2
t = ND N , (1.7)
i =1
есть вероятность того, что, если гипотеза верна, максимальная
разность ∆F(x) превзойдет наблюдаемое значение |DN|; если вероятность (1.7) велика, нет оснований отвергнуть гипотезу. Существуют
также таблицы квантилей распределения Колмогорова для доверительных вероятностей α = 0,10; 0,05; 0,01 и различных размеров
выборки N [3]. Гипотеза ошибочно отвергается, когда она верна, с
вероятностью α ≈ 1 –p. Это правило не зависит от проверяемой гипотезы.
Пример 1.2. Проверка MATLAB – функции RAND, предназначенной для генерирования чисел x∈R(0,1). Программа строит эмпирическую функцию распределения оператором STEM (рис. 1.3)
и рассчитывает вероятность (1.7):
N = 100 % размер выборки
x = rand(1,N);
X = sort(x) % вариационный ряд
ff = ones(1,N)/N;
FN = cumsum(ff); % эмпирическая функция распределения
F = X % гипотетическая функция распределения
stem(X,FN)
hold on
plot(X,F,’r’)
DN = max(abs(FN-F))
DNN = DN*sqrt(N) % статистика Колмогорова
s=0
for i = 1:50
s = s+(-1)^(i-1)*exp(-2*i^2*DN^2);
15
end
p = 2*s % вероятность (1.7)
FN, F
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
x
Рис. 1.3. Гипотетическая и эмпирическая функции
равномерного распределения
Максимальное отклонение DN < 0,0479 < 0,134 = t0(100;0,05) – не
превосходит квантиль t0, соответствующую N = 100 и α = 0,05.
С вероятностью (1.7) p = 0,9759 статистика t = ND N , если гипотеза верна, превзойдет полученную. Эти результаты не дают оснований отвергнуть гипотезу о том, что функция RAND генерирует
равномерно распределенные числа.
В монографии [4] приведены оценки сверху критических значений DN асимптотического распределения Колмогорова, рекомендуемые для N ≥ 80: статистика DN < 1,3581/ N для α = 0,05 или
DN < 1,6276/ N для α = 0,01 не противоречит гипотезе. В условиях примера 1.2 максимальные разности DN = 0,1358 или DN =
= 0,1628 также не дали бы оснований отвергнуть гипотезу.
Суммирование значений гистограммы дает выборочную функцию распределения FS(x) [2] – возрастающую в узлах гистограммы
ступенчатую функцию, отличающуюся от эмпирической функции
распределения FN(x) (рис. 1.2) тем, что ее приращения становятся неравными. Как функция FN(x), функция FS(x) изменяется от
16
нуля до единицы, не убывает и непрерывна справа. В отличие от
FN(x) функция FS(x) сходится по вероятности к функции распределения F(x) не во всей области определения, а в узлах гистограммы.
Использование выборочной функции распределения в качестве
(x) может привести к ошибочным результатам. Функция
оценки F
FS(x) описывает выборочные значения, сгруппированные гистограммой, то есть дискретную случайную величину, принимающую
k значений. При этом условия сходимости по вероятности к непрерывной функции F(x) не выполняются. Кроме того, часто гистограмма строится с относительно небольшим числом интервалов,
как правило, значительно меньшим рекомендуемых для критерия
Колмогорова нескольких десятков значений эмпирической функции распределения.
Пример 1.3. Проверка гипотезы хи-квадрат распределения с тремя степенями свободы с функцией распределения (рис. 1.4)
x
 t2 


2x
x
1
F(x) = 2Φ  x  −
exp  −  − 1, x ≥ 0, Φ(x) =
exp  − dt  –
π
 2
2π −∞


 2
∫
интеграл вероятности. χ2(3) – выборка X формируется фрагментом
программы
FN,F
а)
FS,F
б)
DN = 0,1181
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
2
4
6
8
x
10 12
DN = 0,1953
1
0
0
2
4
6
8
10 12
x
Рис. 1.4. К примеру 1.3
17
N = 100
x1 = randn(1,N); x2 = randn(1,N); x3 = randn(1,N);
X = x1.^2+x2.^2+x3.^2
На рис. 1.4, а показана эмпирическая функция распределения,
сравнение которой с функцией F(x) дает значение DN = 0,1181 <
0,1358 – с вероятностью p = 0,95 нет оснований отвергнуть гипотезу. На рис. 1.4, б приведены значения выборочной функции распределения, построенной по гистограмме с k = 9. Критерий DN =
= 0,1953 таков, что гипотезу формально следует отвергнуть, хотя
не ясно, применим ли здесь критический уровень, пригодный при
N ≥ 80.
1.3. Критерий χ2 К. Пирсона
Пусть гистограмма (1.1) построена по независимой выборке,
извлеченной из генеральной совокупности с плотностью распределения f(x). Вероятность попадания одного выборочного значения в
интервал ∆i равна
pi =
∫ f (x)dx.
∆i
Вероятность попадания в этот интервал ni независимых значений
имеет биномиальное распределение
ni ni
p ni = C N
pi (1 − pi ) N −ni .
Вероятность получения конкретной гистограммы (совместная вероятность попадания в i-й интервал ровно ni значений, i = 1, ..., k)
имеет полиномиальное (мультиномиальное [3,4]) распределение
[11,12]




p H  =




N!
p  n1, ..., nk  =
p1n1 p2n2 ... p knk . 
 n1 ! n2 !...nk !
(1.8)
Например, вероятность (1.8) гистограммы, показанной на рис. 1.1,
pi ≈ 3·10 –8, так что можно прикинуть, сколько раз надо запустить
программу, генерирующую функцию RANDN, чтобы гистограмма
повторилась.
Важно, что величины ni в полиномиальном распределении имеют асимптотически k-мерное нормальное распределение, причем
18
квадратичная форма в экспоненциальном показателе этого распределения
2
XN
=
k
(ni − Npi ) 2
Npi
i =1
∑
(1.9)
имеет распределение χ2 с k –1 степенью свободы [3,4]. Статистика
(1.9) используется в χ2-критерии согласия (К. Пирсон, 1900 г.).
Ошибка приближения, когда XN2 считается распределенной как
χ2(k –1) при конечном размере N, минимальна, когда вероятности pi одинаковы (pi = 1/k), что достижимо при различной ширине интервалов гистограммы [4]. На практике критерий согласия
χ2 Пирсона применяется с постоянной шириной интервалов при
N > 50, при этом рекомендуется, чтобы величины Npi ≥ 5 [2]. Как и
критерий Колмогорова, критерий Пирсона не зависит от проверяемой гипотезы.
Пример 1.4. Проверка гипотезы H о распределении Релея: если
независимые x1, x2∈N(0,σ), то y =
f (y) =
x12 + x 22 ∈Re(σ2), то есть
 y2 
y
exp
 − 2 . σ2
 2σ 
(1.10)
Пусть выборка Y формируется из двух выборок X1 и X2, задаваемых
функцией RANDN. По гистограмме выборки Y рассчитывается статистика (1.9) и сравнивается с ее критическим уровнем X02(α,k – 1),
вычисляемым функцией CHI2INV(1 –α,k –1) [13] как квантиль
χ2(k –1)-распределения для вероятности 1 –α. Гипотеза H о том, что
сформированная выборка имеет распределение Релея с плотностью
(1.10), отвергается, если статистика XN2 > X02(α,k – 1).
N = 50
x1 = randn(1,N); x2 = randn(1,N);
y = sqrt(x1.^2+x2.^2); Xma = 3.5
k = 8 delta = Xma/k
xx = 0:delta:Xma
x = delta/2:delta:Xma F = raylcdf(xx,1) p = diff(F) n = hist(y,x) % релеевская выборка с параметром sigma = 1
% число интервалов гистограммы
% середины интервалов
% функция распределения Релея
% вероятности попадания в интервалы
% числа попаданий
19
h = n/N/delta subplot(1,2,1),stem(x,h)
hold on
plot(x,raylpdf(x,1)) s=0
for i = 1:k
s = s+(n(i)-N*p(i))^2/N/p(i);
end
x2 = s x20 = chi2inv(0.95,k-1) % гистограмма
% график плотности Релея
% статистика (1.9)
% критический уровень для alfa = 0.05
Два примера показаны на рис. 1.5. В обоих случаях статистика оказывается меньше критического уровня для вероятности α = 0,05.
Это означает, что, если гипотеза H верна, с вероятностью p = 1 – α =
= 0,95 статистика окажется меньше критического уровня. Именно
такими оказались статистики в экспериментах, так что оснований
отвергнуть гипотезу H нет.
а)
h,f
б)
h,f
0.6
N = 50
k=8
0.7
0.5
x2 = 2.3849
x20 = 14.0671
0.6
N = 100
k = 15
x2 = 16.8652
x20 = 23.6848
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0
0.1
0
1
2
3
x
0
0
1
2
3
x
Рис. 1.5. К примеру 1.4
В примерах 1.2 и 1.4 рассматривались простые гипотезы H: границы равномерного распределения a = 0, b = 1 и параметр σ2 = 1
плотности (1.10) были известны. Критерии согласия Колмогорова
и Пирсона применимы и в более сложной ситуации, когда проверяется гипотеза о принадлежности функции распределения выборки
FN(x) некоторому параметрическому семейству функций FN(x;θ),
20
зависящих от некоторых неизвестных параметров θ (сложные гипотезы) [2 –4]. В таком случае неизвестные параметры заменяются
их оценками θ . На практике часто используют оценки максимального правдоподобия θ (МП-оценки), полученные по выборке X.
При замене статистики Колмогорова и Пирсона записываются
D = sup
F (x) − F(x; θ ) ,
N
x <∞
N
2

 
 ni − Npi (θ) 
2 

 .
X N (θ) =
)
Np
(
θ
i =1
i
k
∑
(1.11)
При проверке по критерию χ2 применяют результат, полученный Р.
Фишером (1924 г.), считая, что статистика (1.11) имеет асимптотическое (при N → ∞) χ2(k –1 –r) – распределение (r – число оцениваемых параметров). Однако эта методика проверки сложных гипотез
неточна: МП-оценки, полученные по выборке X, могут исказить
распределения статистик. Необходимо использовать специальные
оценки параметров, что существенно усложняет задачу проверки
гипотез [3,4].
В настоящем пособии проверка гипотез применяется только
для проверки эффективности рассматриваемых методов оценивания. Гипотетическое распределение считается неизвестным, так
что предположения о его зависимости от каких-либо параметров
неуместны. В этом смысле проверяемые далее гипотезы – простые,
поэтому проверка сложных гипотез далее не обсуждается.
1.4. Проверка гипотезы однородности
В этой задаче требуется выяснить принадлежность двух независимых выборок X = {x1, ..., xN} и Y = {x = y1, ..., yM} одной генеральной совокупности или разным. В общем виде задача формулируется следующим образом. Пусть выборка X извлечена из генеральной
совокупности с функцией распределения F1(x), выборка Y – из совокупности с функцией распределения F2(x). Гипотеза H однородности выборок: F1(x) = F2(x).
Классические методы проверки гипотезы однородности базируются на критериях Колмогорова и Пирсона. Разработаны также
другие непараметрические методы решения этой задачи [3,5,12].
Если функции распределения F1(x) и F2(x) непрерывны, для
проверки гипотезы используется статистика вида (1.5)
21
D NM = sup x <∞ FN (x) − FM (x) , (1.12)
в которой FN(x), FM(x) – эмпирические функции распределения, построенные по выборкам X и Y. По теореме Н. В. Смирнова
(1944 г.)
 NM

D NM > t | H  = 1 − K(t) = α,
lim N →∞ p 
N
M
+
M →∞ 

K(x) – функция распределения Колмогорова. Критерий однородности Смирнова: гипотеза однородности выборок X и Y отвергается, если для вычисленного значения t статистики (1.12) выполняется неравенство
t ≥ t0 =
N+M
tα , NM
(1.13)
tα – квантиль распределения Колмогорова, соответствующая вероятности p = 1 – α. При этом вероятность ошибочно отвергнуть правильную гипотезу равна α.
Для проверки данных, имеющих дискретную структуру, используют критерий однородности χ2. Дискретную структуру
приобретают в том числе непрерывные выборки, сгруппированные в гистограммы. Критерий однородности χ2 также не зависит
от проверяемых (неизвестных) функций распределения. Применительно к сравнению гистограмм задачу можно поставить следующим образом. L выборок Xi одинакового размера группируются в гистограммы Hi с одними и теми же интервалами. Если
выборки принадлежат одной и той же совокупности с функцией
распределения (неизвестной) F(x), то гипотеза однородности означает
[pl1, pl2, ..., plk ]= [p1, p2, ..., pk ],
l = 1, ..., L,
P = [p1, p2, ..., p k ]  – некоторый (неизвестный) вектор вероятносk
тей,
∑ pi = 1. Например, в табл. 1.2 приведены количества n1, n2,
i =1
n3 попаданий в интервалы шириной ∆ = 0,5 полученные по трем
выборкам из N(0,1) (функция RANDN) размером N = 256. Математическое ожидание числа попаданий в интервал ∆i = Npi = n0.
22
Таблица 1.2
n1
0
3
7
14
25
43
53
60
24
17
7
3
0
n2
1
2
9
16
29
37
48
57
36
16
4
1
0
n3
3
0
7
16
35
40
49
39
43
14
7
2
1
n0 0,62 2,37 7,13 16,8 31,0 44,7 50,5 44,7 31,0 16,8 7,13 2,37 0,62
В качестве меры отклонений наблюдений от n0 по критерию χ2 следовало бы назначить статистику вида (1.9)
2
XN
=
L
k
(nli − Npi ) 2
,
Npi
l =1 i =1
∑∑
однако в ней неизвестные вероятности pi надо заменить их оценкаpi. МП-оценки вероятностей [3]
ми 
L
1

pi =
nli
NL l =1
∑
приводят статистику к виду
2
=L
XN
L
2
k
L
nli
− L, ni = nli. l =1 i =1 p ini
l =1
∑∑
При L = 2 статистика (1.14) записывается
2
XN
=
∑
(1.14)
k
(n1i − n2i ) 2
.
i =1 n1i + n2i
∑
(1.15)
Гипотеза однородности отвергается, если значение t статистики
(1.15)
t ≥ χ 12−α (k − 1)
превышает значение (1 –α) – квантили хи-квадрат распределения с
k –1 степенью свободы. Вероятность ошибочно отвергнуть правильную гипотезу равна α.
Пример 1.5. Проверка однородности двух выборок размером
N = 100 из N(0,1). В программе два блока: в первом вычисляется
статистика (1.12) проверки гипотезы однородности по критерию
Смирнова, во втором – статистика (1.15) проверки по критерию
χ2. Так как эмпирическая функция распределения возрастает на
величину 1/N в точках xj вариационного ряда, статистика (1.12)
23
вычисляется поиском на оси x максимального значения разности
d = |n1 –n2|, n1, n2 – число точек двух вариационных рядов, предшествующих текущему значению x. Во втором блоке выборки группируются в гистограммы, из которых исключаются совпадающие нулевые значения, которые привели бы при вычислении статистики
(1.15) к делению на ноль.
N = 100
x1 = randn(1,N); x2 = randn(1,N);
X1 = sort(x1); X2 = sort(x2); x = -3:0.05:3;
n = length(x)
for i = 1:n
n1 = length(find(X1<x(i)));
n2 = length(find(X2<x(i)));
d(i) = (n1-n2)/N;
end
DN = max(abs(d)) xx = -3:1/2:3 n1 = hist(x1,xx); h1 = n1/N n2 = hist(x2,xx); h2 = n2/N k = find(h1 = = 0&h2 = = 0) n1(k) = []; n2(k) = []; kk = length(n1) s = sum((n1-n2).^2./(n1+n2)) s0 = chi2inv(0.95,kk-1) % вариационные ряды
% статистика (1.12)
% другая числовая ось
% первая гистограмма
% вторая гистограмма
% парные нулевые значения
% исключение парных нулей
% число степеней свободы
% статистика (1.15)
% критический уровень
Статистика (1.12) равна DN = 0,10. Квантиль распределения Колмогорова, соответствующая вероятности p = 0,95, равна t0,95(99) =
= 0,134, так что критический уровень в (1.13) t0 = 0,134 0,02 =
= 0,019 < t = DN = 0,10 – гипотеза однородности по критерию Смирнова отвергается. Числа попаданий в интервалы с исключенными
совпадающими нулевыми значениями приведены в табл. 1.3.
Таблица 1.3
n1
4
5
13
15
16
21
15
9
2
0
n2
3
8
13
18
21
16
14
1
4
2
Статистика (1.15) равна t = XN2 = 7,667 < χ 20,95 (11) = 19,675 – по
критерию χ2 гипотеза однородности не отвергается. Противопо24
ложность результатов можно объяснить тем, что в гистограммах
наблюдения группируются и сглаживаются, так что вероятность
резких расхождений уменьшается.
1.5. Критерий ω2
Эмпирическая функция распределения (1.4) выборки X
0
n
FN (x) = 
N
1
при
x < x min ,
при x n ≤ x < x n +1,
при
x ≥ x max,
xj – значения вариационного ряда, j = 1, ..., N, есть относительное
число выборочных значений xj ≤ x. Фундаментальное свойство эмпирической функции распределения FN(x) – сходимость к функции
распределения F(x) генеральной совокупности


p  lim N →∞ FN (x) = F(x)



 = 1.

Статистика


W = ∫  FN (x) − F(x)
−∞ 

2
∞
2


 f (x)dx 
(1.16)
была бы естественной квадратичной мерой расхождения между FN(x) и F(x). Однако найти функцию распределения статистики (1.16) сложно. Единственный компактный результат получен
Н. В. Смирновым (1936 г.): асимптотическая характеристическая
функция [11]



2it
.
ϕ(t) = lim N →∞ Μ exp  itNW 2   =
sin
2it



(1.17)
Преобразование Фурье функции (1.17) неизвестно, поэтому предельную функцию распределения F(NW2) можно получить только
численно, что и сделано Т. Андерсоном и Д. Дарлингом (1952 г.)
[4]. Для типовых значений доверительной вероятности α ее квантили t0(1 –α) (значения, которые статистика 1.16 превышает с вероятностями α) приведены в табл. 1.4.
25
Таблица 1.4
α
t0(1 –α)
0,100
0,050
0,010
0,001
0,347
0,461
0,743
1,168
Статистика (1.16) приводится к расчетному виду [2,4]
N 

2j − 1
1
Nω =
+  F (x j ) −
12N j =1 
2N
2
∑
и сравнивается с критическим уровнем
2

 
(1.18)
ω 20 = t0 (1 − α).
В учебнике [2] отмечается, что критерий ω2 целесообразно применять в тех случаях, когда отклонение (1.5) эмпирической функции распределения от теоретической невелико, но функции различаются на достаточно большом промежутке. Последнее замечание
носит неопределенный характер, так как обычно эти функции различаются во всей области определения. На практике критерий ω2
применяется редко.
26
2. Ядерные оценки плотности распределения
Теоретическое достоинство ядерных оценок (4)
f (x) =
1
Na N
N
 x − xj 
−
aN 
∑ ϕ
j =1
(2.1)
сходимость по вероятности к непрерывной плотности распределения генеральной совокупности f(x) [3]:


P  lim N →∞ sup x <∞ f (x) − f (x) = 0



 = 1.

2.1. Метод парзеновских окон
Метод оценивания плотности распределения сглаживанием случайных наблюдений был предложен Е. Парзеном (1962 г.) и получил название парзеновских окон [6] (метода Парзена [7]). Каждое
выборочное значение xj сглаживается функцией j(x;xj) – парзеновским окном. Оценка Парзена непрерывной плотности распределения f(x) записывается в виде (2.1) – является ядерной оценкой.
Функции окон задаются как некоторые плотности распределения,
тогда (2.1) – нормированная оценка. Если ядро – нормальная плотность с с.к.о. σ = λ, оценка Парзена
2
N
 (x − x j )  1
1
f (x) =
exp −
−
N j =1 2πλ
2λ 2 

∑
(2.2)
нормированная усредненная сумма нормальных плотностей с модами в случайных точках xj:
∞
N
1
1 = 1.
f (x)dx =
N j =1
−∞
∫
∑
2
Коэффициент a N
в (2.1) (дисперсия λ2 в (2.2)) задается зависимым от размера выборки: λ2 = L2/N. Значение λ определяет ширину
сглаживающего окна, которая, в свою очередь, задает «гладкость»
функции f (x). На рис. 2.1 показано формирование оценки (2.2):
суммируются сглаживающие функции (ядра, окна Парзена) – нормальные плотности с с.к.о λ = 1/ 8, L = 1, с модами в случайных
точках xj. Выборка X∈N(0,1) (рис. 2.1, кривая 1), оценка плотности – функция, показанная на рис. 2.1 (кривая 2).
27
fi, f
0.4
0.35
1
2
0.3
0.25
0.2
0.15
0.1
0.05
0
-3
-2
xmin
-1
0
1
xmax
2
3
x
Рис. 2.1. Оценивание плотности методом Парзена, λ = 1/ 8
В монографии [6] подчеркивается, что сглаживающие функции применяются для интерполяции, хотя не ясно, что интерполируется. Интерполяция – вычисление промежуточных значений
функции, заданной в узловых точках, причем интерполирующая
функция в узлах равна интерполируемой [14]. Здесь же задаются
только узловые точки xj. Утверждается [7], что парзеновская оценка есть та же гистограммная, только каждый столбик гистограммы «расплывается» и, накладываясь на другие, дает более гладкую
аппроксимацию плотности. Термин «аппроксимация», определяющий поиск наилучшего приближения в некотором классе функций
[14] здесь уместен, однако в методе Парзена гистограмма вообще не
используется. По-видимому, наиболее ясна трактовка оценки Парзена как ядерной.
Пример 2.1. Оценивание методом Парзена стандартной нормальной плотности (xj∈N(0,1)) с нормальными функциями окна. Блок
программы построения оценки для одного значения параметра λ:
N = 16 x = randn(1,N); L = 1/4
lamda = L/sqrt(N) 28
% размер выборки
% выборка
% значение параметра
del = 0.1
xx = -3:del:3; % ось абсцисс для оценки Парзена
nn = length(xx)
ff = zeros(1,nn);
for j = 1:N
ff = ff+normpdf(xx,x(j),lamda);
end
f = ff/N % оценка Парзена
trapz(xx,f) % проверка нормированности оценки
subplot(1,2,1),plot(xx,f,xx,normpdf(xx,0,1))
pause
На рис. 2.2 показаны оценки Парзена стандартной нормальной
плотности при N = 16 для различных значений параметра L. Малая ширина окон (L ≤ 1) за счет их слабого пересечения приводит
к сильно флюктуирующей оценке. Чрезмерное расширение окон
(L = 4) искажает оценку в сторону увеличения дисперсии распределения, так что на интервале ( –3,5 < x < 3,5) плотность теряет
3.5
свойство нормированности (
∫
f (x)dx < 1 ).
−3.5
Оценки (рис. 2.2) подобны оценкам, приведенным в монографии
[6, разд. 4.3.4] в логарифмическом масштабе.
1.4
1.2
1
0.8
0.6
0.4
0.2
0
f
0.4
L = 1/4
L=1
0.3
0.2
0.1
-2
0
2
x
f
0.4
0
0.2
0.1
0.1
0
2
0
x
0
2
x
f
L=4
0.3
0.2
-2
-2
0.4
L=2
0.3
0
f
0.5
-2
0
2
x
Рис. 2.2. Оценка нормальной плотности методом Парзена, N = 16
29
Пример 2.2. Оценивание (2.2) с параметром ширины λ = 2/ 128
стандартной нормальной плотности по выборке размером N = 128.
Программа проверки свойств оценки Парзена f (x) включает проверку гипотезы о нормальности оценки по критерию Колмогорова, а также проверку гипотезы однородности исходной выборки X
 , соответствующей оценке Парзена, по критерию χ2.
и выборки X

Выборка X имитируется числом выборочных значений, которые в
соответствии с плотностью f (x) попали бы в интервалы гистограммы шириной j = de/del = 5 интервалов дискретизации del исходной
числовой оси xx.
NN = 128 % размер выборки
x = randn(1,NN); % выборка
L=2
lamda = L/sqrt(NN) % значение параметра ширины окна
del = 0.1
xx = -3:del:3; % исходная числовая ось
nn = length(xx)
ff = zeros(1,nn);
for j = 1:NN
ff = ff+normpdf(xx,x(j),lamda);
end
f = ff/NN % оценка Парзена
trapz(xx,f) % проверка нормированности оценки
subplot(1,2,1),plot(xx,f,xx,normpdf(xx,0,1))
FP = cumsum(f)*del % оценка функции распределения
F = normcdf(xx,0,1)
subplot(1,2,2),plot(xx,FP,xx,F)
Dn = max(abs(FP-F)) % статистика (1.5)
pause
de = 0.5 % числовая ось в другом масштабе
xxx = -3:de:3;
nnn = length(xxx)
j = de/del % отношение масштабов
for i = 1:nnn-1
FF(i) = FP(i*j); % отсчеты оценки функции распределения через j
end
FF = [FP(1),FF] % “разреженная” функция распределения
N = round(diff(FF)*NN) % имитация ненормированной гистограммы
sum(N)
for i = 1:nnn-1
30
xxxx(i) = (xxx(i)+xxx(i+1))/2; end
n = hist(x,xxxx) NP = [N;n]
plot(xx,normpdf(xx,0,1),’r’)
hold on
bar(xxxx,n/NN/de)
stem(xxxx,N/NN/de)
kk = find(n = = 0&N = = 0) kkk = length(kk)
N(kk) = []; n(kk) = []; t = sum((N-n).^2./(N+n)) t0 = chi2inv(0.95,nnn-2-kkk) % середины интервалов гистограммы
% ненормированная гистограмма выборки
% парные нулевые значения гистограмм
% исключение нулевых значений
% хи-квадрат – статистика однородности
% критический уровень
(x) привеОценки плотности (2.2) f (x) и функции распределения F

дены на рис. 2.3. Оценка F(x) получена интегрированием оценки
f (x) по формуле прямоугольников (функцией CUMSUM). Из сходиp
мости по вероятности оценки Парзена f (x) 
→ f (x) следует сходимость по распределению (сходимость в каждой точке непрерывной
(x) → F(x). Это дает возможность для проверки
F(x) [3]) функции F
оценки применить критерий Колмогорова: полученное значение
f
F
1
0.5
0.9
0.8
0.4
0.7
0.6
0.3
0.5
0.4
0.2
0.3
0.2
0.1
0.1
0
-2
0
2
x
0
-2
0
2
x
Рис. 2.3. Оценка нормальной плотности методом Парзена, N = 128
31
статистики DN = 0,0607 < 1,3581/ N = 0,1201 с вероятностью p =
= 0,95 не дает оснований отвергнуть гипотезу нормальности оценки
[4]. Реальный смысл этого результата: выборка, по которой построена гистограмма, могла быть извлечена из генеральной совокуп(x).
ности с функцией распределения F
Разбиение числовой оси на интервалы ∆, в j = 5 раз более широкие, чем интервал дискретизации del = 0,1, позволяет построить
(x) имитировать гисгистограмму выборки X, а также по оценке F
 . В табл. 2.1 приведены значения чисел попатограмму выборки X
дания ni и Ni в интервалы гистограмм.
Таблица 2.1
n
N
0
0
3
2
6
7
10 17 17
10 15 17
31
29
29
28
12
15
1
2
2
1
0
1
На рис. 2.4 столбиками изображены гистограмма h (исходная выборка, функция BAR), штрихами – имитирующая гистограмма H
(функция STEM) на фоне плотности N(0,1). После удаления нулевого первого столбца гистограмм рассчитывается χ2 – статистика
�:
(1.15) проверки однородности выборок X и X
2
t = XN
= 2,4861 < χ 20,95 (10) = t0 = 18,3070 −
гипотеза не отвергается.
h, H
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
x
-2
-1
0
1
2
Рис. 2.4. Имитация проверки однородности
32
Реальный смысл последнего результата: есть основания полагать, что выборка X1∈N(0,1) и выборка X2 из генеральной совокупности с плотностью распределения f (x) в виде суммы (2.2) ста двадцати восьми нормальных плотностей с модами в точках X были
бы неразличимы по критерию однородности. Непосредственная
проверка этого предположения потребовала бы построения специального генератора независимых случайных величин с плотностью
распределения f (x) (рис. 2.3).
Окна j(x) в оценках Парзена могут задаваться не только в виде
нормальных, но и других подходящих плотностей распределения.
Ширина окон подбирается в каждом случае.
Основной недостаток метода – необходимость описания получаемой графической оценки формулой, без которого дальнейшее
использование оценки может стать затруднительным. Кроме того,
оценки Парзена заметно флюктуируют при размерах выборки в сотни-тысячи наблюдений [6].
2.2. Метод n ближайших соседей
Оценка Парзена с фиксированной шириной окна λ = L/ N , зависящей от размера выборки, особенно эффективна для плотностей
типа равномерных, когда концентрация выборочных значений xj
на единице длины числовой оси в среднем изменяется мало. При
неравномерных плотностях распределения естественно задавать
ширину окна функцией плотностей, например обратно пропорциональной значениям плотности. Этот принцип зависимости ширины
окна не от размера, а от характера выборки реализован оцениванием методом n ближайших соседей [6,7].
Пусть некоторый интервал ∆ расположен в области моды плотности и его размер задан таким, чтобы вместить n выборочных значений xj. Если центр интервала в точке x, то n значений будут ближайшими соседями x. Конечно, ширина интервала в среднем есть
функция заданного числа n. Вблизи моды ширина интервала относительно невелика, что приведет к хорошему разрешению. По мере
удаления от моды плотности ширина интервала растет. Если плотность многомодальна, то вблизи следующей моды ширина уменьшается и т. д. В монографии [7] метод n ближайших соседей считается
предназначенным для оценивания многомодальных плотностей.
Вероятность попадания в i-й интервал оценивается отношением
pi =
ni
.
N∆ i
(2.3)
33
Необходимые и достаточные условия сходимости вероятности (2.3)
к непрерывной оцениваемой плотности f(x) сформулированы в [6]:
lim N →∞ ni = ∞, lim N →∞
ni
= 0.
N
Смысл первого условия в том, что при N = ∞ количество ni выборочных значений, попадающих в i-й интервал, не должно быть конечным, иначе pi = 0. В то же время, при N → ∞ число ni должно увеличиваться достаточно медленно, чтобы ширина интервала уменьшалась до нуля – в этом смысл второго условия.
Если задать ni = N , ширина интервала
∆i =
1
,
N pi
как и в методе Парзена, имеет вид ∆i = L/ N , но значение L = 1/pi
теперь обратно пропорционально вероятности pi.
Если бы в соседние интервалы ∆i и ∆i+1 включались группы по n
различных выборочных значений, то есть каждое значение xj принадлежало бы одному и только одному интервалу, разбиение выборки X было бы гистограммой с равновероятными интервалами.
Как указывалось в подразд. 1.3, такая гистограмма обеспечила бы
минимальную ошибку описания статистики (1.8) χ2(k –1)-распределением [4]. Однако метод n ближайших соседей такую группировку не предусматривает, одно и то же xj включается в несколько
интервалов, поэтому отнесение этого метода в монографиях [6,7]
к гистограммным представляется искусственным.
Группировку по n ближайшим выборочным значениям можно
организовать разными способами. Один из простых алгоритмов:
1) выборка XT = [x1, ..., xN] упорядочивается в порядке возрастания;
2) для группы значений (x1, ..., xn) вычисляются среднее x 1 и
с.к.о. σ1;
3) пункт 2 повторяется для (x2, ..., xn+1),…, (xN –n+1, ..., xN) – вычисляются x 2, ..., x N −n +1, σ2, ..., σN –n+1;
4) вычисляется оценка плотности вида (2.2)
1
f (x) =
Nn
Nn
∑
i =1
fi (x) =
1
Nn
Nn
∑
i =1
 (x − x i ) 2 
1
exp −
, N n = N − n + 1. (2.4)
2σ i2 
2πσ i

Ширина окон переменна по выборке, но число выборочных значений, попадающих в окна, здесь также может изменяться.
34
Пример 2.3. На рис. 2.5 для N = 16 показана реализация X∈N(0,1)
(числа с номерами от 1 до 16). Числа группируются по четыре (тринадцать групп), для каждой группы вычисляются м.о. и с.к.о. mi
и σi и строятся тринадцать нормальных плотностей N(mi, σi) – слагаемые суммы (2.4). В первое окно f1(x) попадают шесть выборочных значений x1, ..., x6, во второе – четыре значения x2, ..., x5,
в третье – семь значений x2, ..., x8, в тринадцатое – четыре значения x13, ..., x16.
f
f2
f13
f1
x
0
1
2 3 4
1
6
7 8
14 15 16
x
3
13
Рис. 2.5. Слагаемые суммы (2.4)
Примеры оценки (2.4) стандартной нормальной плотности приведены на рис. 2.6. Для N = 16 – рис. 2.6, а, для N = 256 (241 группа
по 16 чисел) – рис. 2.6, б. Характер оценок такой же, как у приведенных в работе [6]. Среди N чисел встречаются группы близкорасположенных, с.к.о. в таких группах мало, соответствующие слагаемые дают выбросы оценки. Увеличение размера выборки слабо
компенсирует это явление. Соответствующие оценки функции распределения показаны на рис. 2.6, в, 2.6, г. Несмотря на флюктуирующий характер оценок, статистики (1.4) равны DN = 0,2370,
DN = 0,0696: в обоих случаях гипотеза нормальности не отвергается с вероятностью p = 0,95, так как статистики оказываются меньше критических уровней 1,3581/ N [4].
35
а)
f
1
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
в)
f
б) 0.8
x
-3
-2
-1
0
1
2
0
3
F
1
г)
0.8
-3
0.4
0.2
0.2
-1
0
1
2
3
x
1
2
3
1
2
3
x
0.8
0.6
-2
0
F
0.4
-3
-1
1
0.6
0
-2
0
-3
-2
-1
0
x
Рис. 2.6. Оценивание методом n ближайших соседей
Пример 2.4. Оценивание двумодальной нормальной плотности
f (x) =
2
 (x − 3) 2 
 (x + 1) 
1
1
exp −
exp −
+
. 2  2 2π
2 
2 2π



(2.5)
В программе выделяются группы по k = 8 соседних выборочных
значений (nn = N –k+1 = 57 групп) и вычисляются оценка (2.4)
(x),
плотности f (x), соответствующая функция распределения F
а также статистика (1.5) проверки гипотезы принадлежности вы(x), соответствующей плотборки X распределению с функцией F
ности (2.5), по критерию Колмогорова. Гипотеза однородности по
критерию Смирнова проверяется сравнением эмпирической функ(x). Вследствие
ции распределения Fэ(x) выборки X с функцией F
(x) разности ∆F =
несовпадения узлов задания функций Fэ(x) и F
j

(xx ), в ко= Fэ(xj) – F(x j ) заменяются разностями ∆Fj = Fэ(xj) – F
j
торых xxj – точки числовой оси, наиболее близкие к выборочным
значениям xj.
N = 64 k = sqrt(N) 36
% размер выборки
% масштаб окон
del = 0.1
xx = -4:del:6; % числовая ось
nx = length(x)
z1 = randn(1,N*1/2)-1; z2 = randn(1,N*1/2)+3;
Z = sort([z1 z2]); % упорядоченная выборка
nn = N-k+1 % число слагаемых в (2.4)
for i = 1:nn
for j = 1:k
ZZ(j) = Z(i+j-1);
end
Zmm(i) = mean(ZZ);
sii(i) = std(ZZ);
end
Zm = Zmm % оценки м.о.
si = sii % оценки с.к.о.
ff = zeros(1,length(xx));
for i = 1:nn
ff = ff+normpdf(xx,Zm(i),si(i));
end
f = ff/nn; % оценка плотности распределения (2.4)
F1 = cumsum(f)*del; % оценка функции распределения
F = normcdf(xx,-1,1)/2+normcdf(xx,3,1)/2;
subplot(1,2,1),plot(xx,f,xx,normpdf(xx,-1,1)/2+normpdf(xx,3,1)/2)
subplot(1,2,2),plot(xx,F1,xx,F)
DN1 = max(abs(F1-F)) % статистика (1.5)
for j = 1:N
dd = abs(Z(j)-xx(:));
k = find(dd = = min(dd)); % поиск наиболее близких отсчетов оси xx
% и выборочных значений
DNN(j) = F1(k)-j/N; % имитация разностей в (1.12)
end
DN = max(abs(DNN)) % статистика (1.12)
Оценки плотности и функции распределения приведены на
рис. 2.7. Статистики (1.5) и (1.12), равные DN = 0,0743 и DN =
= 0,0597, не дают оснований отвергнуть гипотезы нормальности и
однородности.
Недостатки метода n ближайших соседей аналогичны недостаткам метода Парзена.
Рассмотренные в этом разделе ядерные оценки плотности распределения – непараметрические. В методах Парзена и n ближай37
f
F
1
0.9
0.3
0.8
0.25
0.7
0.6
0.2
0.5
0.15
0.4
0.3
0.1
0.2
0.05
0.1
0
-4
-2
0
2
4
6
xx
0
xx
-4
-2
0
2
4
6
Рис. 2.7. Оценивание двумодальной нормальной плотности, N = 64
ших соседей не используются никакие предположения о моментах
или других параметрах оцениваемых распределений.
38
3. Гистограммное оценивание
Гистограмма (рис. 1.1) – ступенчатая функция, непараметрический статистический аналог непрерывной плотности распределения
f(x). Она является численной оценкой вероятностей попадания pi
случайной величины x в интервалы ∆i. Гистограммная оценка f (x)
плотности f(x) – некоторая непрерывная функция, описывающая
гистограмму, с ограничениями
f (x) ≥ 0,
∞
∫ f (x)dx = 1.
−∞
3.1. Нормальная гистограммная аппроксимация
Гистограммную оценку можно получить, умножая значения
гистограммы ni/N на базисные нормальные плотности распределения f0(x;xi,σ) с модами в точках xi – серединах интервалов гистограммы:
k
k
ni
ni
 (x − x i ) 2 
1
exp −
(3.1)
f (x) =
f0 (x; x j , σ) =
. 2
2σ 2 

i =1 N
i =1 N 2πσ
∑
∑
Сумма (3.1) имеет обычный для функциональных рядов вид: базисные функции f0(x) умножаются на свои коэффициенты. Оценка
(3.1) аппроксимирует значения гистограммы. Она неотрицательна
и нормирована:
∞
k n
j
f (x)dx =
= 1.
N
=
1
j
−∞
∫
∑
Аппроксимация гистограммы суммой (3.1) нормальных плотностей далее называется нормальной аппроксимацией гистограммы.
Нормальная аппроксимация (3.1) отличается от оценки Парзена (2.2) тем, что моды сглаживающих функций f0 привязываются
к гистограмме, а не к выборочным значениям xj. Число слагаемых
в сумме (3.1) равно k << N, что позволяет записать оценку относительно компактной формулой.
Пример 3.1. Реализация X∈N(0,1) размером N = 128 (функция
 = –0,0495,
RANDN) имела значения xmin = –2,2023, xmax = 2,1832, m
σ = 0,8906. Программа нормальной аппроксимации включает расчет гистограммы, соответствующий алгоритму функции HIST.
N = 128
x = randn(1,N);
39
m = mean(x) % оценка м.о. выборки
si = std(x) % оценка с.к.о. выборки
k = 12 % число интервалов гистограммы
xma = max(x); xmi = min(x);
de = (xma-xmi)/k % ширина интервала
h = hist(x,k)/N/de % гистограмма
xx = xmi:de:xma; % границы k интервалов – k+1 точка
for i = 1:k
xs(i) = (xx(i+1)+xx(i))/2; % середины интервалов
end
xi = xs;
xxx = xmi:de/10:xma; % ось для плотности распределения
subplot(2,1,1),bar(xs,h)
hold on
plot(xxx,normpdf(xxx,m,si))
for i = 1:k
fff = normpdf(xxx,xs(i),de/2); % базисные функции
subplot(2,1,2),plot(xxx,fff,’k’)
hold on
end
stem(xs,h*3)
pause
f = zeros(1,length(xxx));
for i = 1:k
f = f+h(i)*normpdf(xxx,xs(i),de/2)*de; % оценка (3.1)
end
plot(xxx,normpdf(xxx,m,si),xxx,f,’r’)
hold on
stem(xs,h)
Гистограмма hi = ni/N∆, i = 1, ..., k, k = 12, ∆ = (xmax –xmin)/k =
= 0,3655 (табл. 3.1) и базисные функции f0(x;xi,σ), σ = ∆/2, приведены на рис. 3.1 (на рис. 3.1, б значения гистограммы, показанные
функцией STEM, увеличены для лучшего восприятия).
Таблица 3.1
i
1
2
3
4
5
6
xi
 –2,0196
 –1,6541
 –1,2887
 –0,9232
 –0,5578
 –0,1923
hi
0,0641
0,0855
0,2565
0,2993
0,2138
0,4275
i
7
8
9
10
11
12
xi
0,1732
0,5386
0,9041
1,2695
1,6350
2,0005
hi
0,4489
0,4062
0,2779
0,1710
0,0641
0,0214
40
а)
f, h
0.5
0.4
0.3
0.2
0.1
0
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
2.5
0.5
1
1.5
2
x
2.5
f0
б) 2.5
2
1.5
1
0.5
0
-2.5
-2
-1.5
-1
-0.5
0
Рис. 3.1. Гистограмма и базисные функции
f
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
2.5
Рис. 3.2. Нормальная аппроксимация гистограммы, N = 128
41
Оценка плотности (3.1), аппроксимирующая гистограмму на оси с
 , σ , ). Гисшагом ∆/10, показана на рис. 3.2. Там же приведена N( m
тограмма уже сглаживает наблюдения X, группируя их, поэтому
нормальная гистограммная аппроксимация дает менее флюктуирующую оценку плотности, чем парзеновская (рис. 2.3).
Метод нормальной аппроксимации чувствителен к ширине
базисных функций. Варианты неудачного задания с.к.о. σ – на
 )): слишком малая шири, σ
рис. 3.3 (на рис. 3.3, б кривая 1 – N( m
на (σ = ∆/4) приводит к появлению флюктуаций оценки, слишком
большая (σ = ∆) дает оценку с явно завышенным рассеиванием.
h, f
а)
h, f
б)
0.7
0.45
0.6
1
0.4
0.35
0.5
0.3
0.4
0.25
0.3
0.2
0.15
0.2
0.1
0.1
0.05
x
0
-2
-1
0
1
2
0
x
-2
-1
0
1
2
Рис. 3.3. Нормальная аппроксимация гистограммы
Чрезмерное расширение базисных функций приводит также к
нарушению свойства нормированности оценки на конечном интервале ∆x ее задания, так как крайние базисные функции целиком не
вписываются в ∆x.
Выбрать подходящее значение с.к.о. базисных функций можно
из следующих соображений.
В середине i-го интервала значение аппроксимирующей функции (3.1) должно быть приблизительно равно значению гистограммы hi:
42
∆2
4∆ 2


−
−
2

1
2
σ
2σ 2 + ... =
+
(
+
)
+
(
+
)
fi =
n
n
n
e
n
n
e
 i

i −1
i +1
i −2
i +2
N 2πσ 



=
ni
n
1
+
S(σ) ≈ i , N∆
N 2πσ N 2πσ
S(σ) = (ni −1 + ni +1)e
−
∆2
2σ 2
+ (ni −2 + ni +2 )e
−
(3.2)
4∆ 2
2σ 2
+ ....
Равенство (3.2) – трансцендентное уравнение относительно σ
L=
1  S(σ) 
2π
,
1 +
=
σ
ni 
∆
(3.3)
показанное для одного из центральных интервалов гистограммы, построенной по выборке из N(0,1) с N = 5000, ∆ = 0,6111, на
рис. 3.4.
8
L
6
4
2
0
0.3
0.5
0.7
0.9
si/del
Рис. 3.4. Трансцендентное уравнение (3.3)
Решение уравнения (3.3) σ ≈ ∆/2 определяет такую ширину базисных функций, что каждая из них «захватывает» три соседних интервала гистограммы (рис. 3.1). Оно мало зависит от числа слагаемых в сумме S(σ) и характерно для центральных интервалов. На
43
краях гистограммы, а также при малых размерах выборки, когда
гистограмма резко меняется от выборки к выборке, оптимальное
значение с.к.о. может увеличиться.
Пример 3.2. Предельные свойства нормальной гистограммной
аппроксимации. Пусть f(x) – плотность нецентрального χ2-распределения [3,12] с четырьмя степенями свободы и параметром нецентральности λ2 = 1 (рис. 3.5, а, функция NCX2PDF [13]), числовая
ось для гистограммы задана с интервалом ∆ = 0,5, середины интервалов xj = 0,25; 0,75;...; 24,75 (k = 50). Если бы выборка была велика (N → ∞), в i-й интервал гистограммы попала бы часть выборки
p
(ni / N) 
→∆Fi = F(xi+1) –F(xi) (рис. 3.5, а, штриховая диаграмма),
F(x) – функция распределения (NCX2CDF). Замена в аппроксимирующей сумме (3.1) коэффициентов ni/N на предельные значения
∆Fi определяет предельное качество оценки при данной ширине
интервалов ∆. Предельная аппроксимирующая функция при σ =
= ∆/2 визуально мало отличается от плотности f(x) (рис. 3.5, б).
Этот результат получен при относительно большой ширине интервалов гистограммы ∆. Уменьшение ширины интервалов еще улучшает оценку (рис. 3.6).
а)
f, dF
б)
f
0.14
0.14
0.12
0.12
0.1
0.1
0.08
0.08
0.06
0.06
0.04
0.04
0.02
0.02
0
0
5
10
15
x
20
0
0
5
10
15
x
20
Рис. 3.5. Плотность нецентрального χ2-распределения,
предельная нормальная гистограммная аппроксимация, ∆ = 0,5
44
f, dF
f
0.14
0.14
0.12
0.12
0.1
0.1
0.08
0.08
0.06
0.06
0.04
0.04
0.02
0.02
0
0
5
10
15
x
20
0
0
5
10
15
x
20
Рис. 3.6. Плотность нецентрального χ2-распределения,
предельная нормальная гистограммная аппроксимация, ∆ = 0,25
Пример 3.2 иллюстрирует приближение нормальной гистограммной аппроксимации f (x) к плотности распределения генеральной
совокупности f(x) при N → ∞, ∆ → 0. При ∆ → 0 с.к.о. базисных функций в (3.1) σ → 0, а нормальные базисные функции приближаются
к δ-функциям Дирака [15]. В пределе конечное множество точек –
середин интервалов xi становится континуумом точек γ числовой
оси (k = ∞), базисная функция f0(x;xj,σ) = δ(x –γ) становится δ-функцией в точке x = γ. Если при этом отношение hi = ni/N∆ → f(γ), то
сумма (3.1) записывается сверткой
∞
∫ f (γ)δ(x − γ)dγ = f (x), (3.4)
−∞
равной плотности f(x) по фильтрующему свойству δ-функции [15].
В математической статистике сходимость гистограммы к плотности распределения (по вероятности или другого вида) не доказана,
поэтому к свертке (3.4) следует относиться в смысле статистической аналогии гистограммы и плотности (подразд. 1.1). Реальный
смысл аппроксимации (3.1) состоит в том, что увеличение размера
45
выборки и уменьшение интервалов гистограммы приближает f (x)
именно к f(x), а не к какой-то другой функции. Статистическая однородность выборок с плотностями f (x) и f(x) – достаточное условие качества оценки f (x).
Пример 3.3. Гистограмма выборки X (N = 100) из нецентрального χ2(4,1)-распределения c k = 30 интервалами показана на рис. 3.7,
а. Аппроксимация (3.1) такой флюктуирующей гистограммы требует большей ширины базисных функций, чем полученная из уравнения (3.3) для выборок большого размера (с.к.о. базисных функций σ ≈ 0,5∆). На рис. 3.7, б кривая 1 – плотность χ2(4,1), кривые 2,
3 – нормальная гистограммная аппроксимация с с.к.о. σ ≈ 0,75∆ и
σ ≈ 1,5∆: незначительный уровень флюктуаций оценки достигается
при увеличении с.к.о. в три раза.
Гистограмма h другой выборки (xmax = 18,4156; ∆ = 0,6139) из
той же генеральной совокупности и оценка плотности f (x) при σ =
= 1,5∆ показаны на рис. 3.8, а и 3.8, б (кривая 1). Размер выборки
N = 100 при k = 30 интервалах недостаточен: по критерию χ2 гипотеза принадлежности выборки нецентральному хи-квадрат-распределению с вероятностью p = 0,95 отвергается, так как статистика критерия χ2 = 52,2168 > χ02 (0,95;29) = 42,5570.
h, f
б)
а)
f
2
0.18
0.2
0.16
0.14
0.15
0.12
1
0.1
0.1
0.08
3
0.06
0.05
2
0.04
0.02
0
0
5
10
15
x
0
0
5
10
15
Рис. 3.7. Нормальная гистограммная аппроксимация, χ2(4,1)
46
x
h. f
а)
f
б)
0.16
0.14
0.12
0.1
1
2
0.1
0.08
0.06
0.05
0.04
0.02
0
0
5
10
15
x
0
x
0
5
10
15
Рис. 3.8. Нормальная гистограммная аппроксимация, χ2(4,1)
(x) (табл. 3.2), полученной интегПо функции распределения F
рированием оценки плотности f (x) по формуле прямоугольников
(функция CUMSUM) вычисляются вероятности попадания 
p i в те
же интервалы ∆i, на которых построена гистограмма h.
Таблица 3.2
№

F
№

F
№

F
1
3
0,0032 0,0990
11
12
0,6314 0,6815
19
20
0,8734 0,8841
5
6
7
8
0,2327
0,3068
0,3832
0,4566
13
14
15
16
0,7277
0,7676
0,8002
0,8263
21
22
25
27
0,8946
0,9056
0,9364
0,9432
9
10
0,5217 0,5788
17
18
0,8463 0,8615
29
30
0,9512 0,9546
 i случайной величины с плотностью f (x) в средЧисло попаданий n

нем равно p i N. В табл. 3.3 приведены значения ni (число попада i, округленное до ближайшего цений выборочных значений) и n
лого функцией ROUND.
47
Таблица 3.3
№
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
n

n
7
5
7
8
5
10
9
6
4
6
4
7
3
3
3
4
5
6
7
7
8
7
7
6
5
5
5
4
3
3
№
16
17
18 19
20
21
22
23
24
25
26
27
28
29
30
n

n
2
1
1
1
1
1
1
2
1
0
0
0
0
0
2
2
2
1
1
1
1
1
1
1
0
0
0
0
0
1
По данным табл. 3.3, из которой исключаются интервалы с номерами 25 –29, вычисляется статистика (1.15) проверки гипотезы однородности
2
t = XN
−5 =
k −5
∑
i =1
 i)2
(ni − n
.
i
n +n
(3.5)
i
В программе
N = 100;k = 30;
x1 = randn(1,N)+1; x2 = randn(1,N); x3 = randn(1,N); x4 = randn(1,N);
x = x1.^2+x2.^2+x3.^2+x4.^2; % выборка с нецентральным хи-квадрат
% распределением
xma = max(x);
de = xma/k; % размер интервала
xx = de/2:de:xma; % середины интервалов
n = hist(x,k); % число попаданий в интервалы гистограммы
h = n/N/de % гистограмма
xxx = 0:de/10:xma; % ось для аппроксимации (3.1)
subplot(1,2,1),bar(xx,h)
hold on
plot(xxx,ncx2pdf(xxx,4,1))
X = 0:de:xma; % ось для функции распределения
F = ncx2cdf(X,4,1);
p = diff(F) % гипотетические вероятности попадания
s=0
for i = 1:k
s = s+(p(i)*N-n(i))^2/p(i)/N;
end
hi2 = s % хи-квадрат-статистика
t0 = chi2inv(0.95,k-1) % критический уровень
48
% pause
fff = zeros(1,length(xxx));
for i = 1:k
fff = fff+n(i)/N*normpdf(xxx,xx(i),1.5*de); % аппроксимация (3.1)
end
P = trapz(xxx,fff) % проверка нормированности аппроксимации
% ffn = fff/P; % нормировка оценки (3.1)
% fff = ffn;
F = cumsum(fff)*de/10 % аппроксимация функции распределения
for i = 1:k+1
FFF(i) = F((i-1)*10+1);
end
FF = FFF; % значения функции распределения на оси
гистограммы
pp = diff(FF); % оценки вероятностей попадания
subplot(1,2,2),plot(xxx,ncx2pdf(xxx,4,1), xxx,fff,xxx,ffn)
nn = round(pp*N); % округление
% nn = pp*N
NN = [n;nn]
kk = find(n = = 0&nn = = 0) % парные нулевые значения гистограмм
kkk = length(kk)
n(kk) = []; nn(kk) = []; % исключение нулевых значений
t = sum((n-nn).^2./(n+nn)) % хи-квадрат-статистика однородности
t0 = chi2inv(0.95,k-1-kkk) % критический уровень
проверяются гипотезы принадлежности выборки распределению
χ2(4,1) и однородности. Статистика (3.5) t = 3,9225 < χ2 (0,95;124) =
= 36,4150 не дает оснований отвергнуть гипотезу однородности.
Аппроксимирующая сумма f (x) оказывается близкой к плотности
f(x), несмотря на неадекватный характер исходной гистограммы.
Это означает, что нормальная гистограммная аппроксимация может оказаться эффективной при небольших размерах выборки.
В этом примере оценка f (x) ненормирована (табл. 3.2). Как отмечалось выше, это следствие расширения базисных функций, которое предпринято для уменьшения флюктуаций оценки. Нормированная оценка вычислена делением f (x) на максимальную вероятность pmax = 0,9546 (табл. 3.2) и показана на рис. 3.8, б (кривая 2).
Статистика (3.5) при этом изменилась незначительно: t = 5,7518.
 (табл. 3.3):
Кроме того, проверено влияние округления значений n
использование дробных значений дает близкое значение статистики t = 5,7518.
49
h,f
f
0.4
0.4
0.4
0.3
0.4
0.2
0.4
0.1
0
x
–4
–2
0
2
0.4
–4
–2
0
2
x
Рис. 3.9. Нормальная гистограммная аппроксимация, t(5)
Пример 3.4. Гистограмма (k = 20) выборки из распределения
Стъюдента [3,12] (t-распределения) с пятью степенями свободы и
оценка плотности TPDF(5) [13] показаны на рис. 3.9.
N = 128;k = 20;
x1 = randn(1,N);x2 = randn(1,N);x3 = randn(1,N);x4 = randn(1,N);x5 = randn(1,N);
x6 = randn(1,N);
x = x6./sqrt((x1.^2+x2.^2+x3.^2+x4.^2+x5.^2)/5); % выборка с t-распределением
Как и в предыдущем примере, размер выборки N = 128 мал для
получения «гладкой» гистограммы с k = 20 интервалами: χ2-критерий предписывает отвергнуть гипотезу принадлежности выборки t(5)-распределению (t = 36,5527 > t0 = χ2 (0,95;19) = 30,1435).
В то же время гипотеза однородности выборок с плотностями f(x) =
= t(5) и f (x) (рис. 3.9) не отвергается (t = 8,0701 < t0).
3.2. Нормальная гистограммная интерполяция
Аппроксимация наблюдений Y методом наименьших квадратов (МНК-аппроксимация) предписывает минимизировать сумму
квадратов разностей наблюдений yi и значений аппроксимирующей функции f (x) в точках xi [3,9,12]
50
s=
∑ 
i
2
yi − f (x i )  .

В стандартной МНК-процедуре вычисляется аппроксимирующий многочлен заданного порядка p. Матрица плана Z формируется
построчно: первая строка – значения многочлена нулевого порядка
в точках xj (z1j = xj0 = 1 – строка единиц); вторая строка – значения
многочлена первого порядка в точках xj (z2j = xj1 = xj); третья строка – значения многочлена второго порядка (z3j = xj2), последняя
строка – z(p+1)j = xjp. Например, для p = 4 матрица плана, заданная
на оси x = –3: 1: 3, записывается
Z=
1
–3
9
–27
81
1
–2
4
–8
16
1
–1
1
–1
1
1
0
0
0
0
1
1
1
1
1
1
2
4
8
16
1
3
9
27
81
Коэффициенты bk аппроксимирующего многочлена
p
f (x) =
∑ bkzk (x) (3.6)
k =0
рассчитываются по формуле [9]
b = (ZZ T ) −1 Y. В точках xi значения многочлена
f (x i ) = Z T b. (3.7)
(3.8)
МНК-функции MATLAB (POLYFIT, POLYVAL) позволяют получить коэффициенты (3.7) на оси вычисления гистограммы (xi – середины интервалов, наблюдения yi = hi) и воспроизвести многочлен
в другом масштабе.
Пример 3.5. Программа МНК-аппроксимации гистограммы,
построенной по выборке стандартной нормальной величины (N =
= 1000, p = 5)
N = 1000
x = randn(1,N);
k = 12
xma = max(x); xmi = min(x);
51
de = (xma-xmi)/k
xx = xmi+de/2:de:xma-de/2;
xxx = xmi:de/10:xma;
n = hist(x,k)
h = n/N/de
p = polyfit(xx,h,5) % коэффициенты МНК-аппроксимации
f = polyval(p,xxx) % МНК-аппроксимация на оси xxx
bar(xx,h)
hold on
plot(xxx,normpdf(xxx,0,1),xxx,f)
stem(xx,h)
Аппроксимирующий многочлен (3.6)
f (x ) = 0,3797 + 0,0131x − 0,1267x 2 − 0,0092x 3 + 0,0104x 4 + 0,0012x 5
показан на рис. 3.10, кривая 2; плотность N(0,1) – рис. 3.10, кривая 1.
На практике используют аппроксимацию многочленом невысокого порядка, иначе аппроксимирующая функция становится громоздкой. При этом f (x) на краях может принимать отрицательные
и другие недопустимые для плотности распределения значения,
как на рис. 3.10.
0.4
0.35
0.3
0.25
0.2
0.15
0.1
2
1
0.05
0
-0.05
-4
x
-3
-2
-1
0
1
2
3
Рис. 3.10. МНК-аппроксимация гистограммы многочленом
пятого порядка
52
Строки матрицы плана могут задаваться не только значениями
многочленов, но и другими линейно-независимыми функциями.
Как в подразд. 3.1, базисными функциями (строками матрицы
плана) можно выбрать нормальные плотности с модами в точках xj
(рис. 3.1). При этом получается квадратная матрица плана. Оценка
плотности (3.8) записывается суммой, аппроксимирующей гистограмму с k интервалами:
f (x) =
k
∑
j =1
b j z j (x), z j (x) =
 (x − x j ) 2 
exp −
, 2σ 2 
2πσ 2

1
(3.9)
bj – МНК-коэффициенты (3.7). Функция (3.9) нормирована, если
сумма коэффициентов (3.7) равна единице:
∞
∫
−∞
f (x)dx =
k
∑ b j = 1,
j =1
Следует отметить, что обычно матрица плана задается с числом
строк, меньшим числа наблюдений k. Квадратная матрица плана –
предельный случай.
Пример 3.6. На рис. 3.11 (кривая 1) показана оценка (3.9) в узлах x = xi для гистограммы, полученной по выборке из генеральной
совокупности со стандартным нормальным распределением (N =
= 1000, k = 12, ∆ = 0,5, σ = ∆/2). На рис. 3.11 (кривая 2) – нормальная плотность. Там же изображена базисная функция z4. В узлах
аппроксимации значения оценки равны значениям гистограммы:
f (x i ) = hi , i = 1, ..., k, (3.10)
то есть МНК-сумма (3.9) оказывается интерполирующей, а не аппроксимирующей.
Можно показать, что равенство (3.10) – следствие общего свойства метода наименьших квадратов: МНК с квадратной симметричной матрицей плана дает интерполирующую функцию f (xi ) = Zb
[9]. Действительно, квадратная матрица Z с базисными функциями (3.9) симметрична: ZT = Z; коэффициенты (3.7)
b = (ZZ T ) −1ZH = (ZZ) −1ZH = Z −1Z −1ZH = Z −1H, в соответствии с (3.8)
(3.11)
f (x i ) = Z T b = Zb = ZZ −1H = H,
53
h, f
2
0.4
1
0.3
0.2
0.1
0
-3
-2
-1
0
1
2
3
1
2
3
x
z4
2
1.5
1
0.5
0
-3
-2
-1
0
x
Рис. 3.11. МНК-интерполяция гистограммы, σ = ∆/2
равенство (3.10) выполняется, то есть МНК-функция (3.9) – интерполирующая. МНК-коэффициенты вычисляются по формуле
(3.11), более простой, чем (3.7). Далее МНК-метод называется нормальной МНК-интерполяцией гистограммы.
Пример 3.7. Интерполяция гистограммы выборки из N(0,1), N =
= 128, k = 12, ∆ = 0.5, σ = ∆/2. Воспроизведение значений оценки
(3.9) в точках между узлами – серединами интервалов гистограммы реализуется программой
N = 128
x = randn(1,N); m = mean(x); si = std(x);
k = 12
xma = max(x); xmi = min(x);
de = (xma-xmi)/k
xx = xmi+de/2:de:xma-de/2; % середины интервалов гистограммы
xxx = xmi:de/10:xma; % числовая ось для интерполяции
n = hist(x,k)
h = n/N/de
for i = 1:k
for j = 1:k
54
z(i,j) = normpdf(xx(j),xx(i),de/2);
end
end
Z = z % матрица плана
b = inv(Z)*h’ % МНК-коэффициенты (3.11)
ff = 0
for j = 1:12
ff = ff+b(j)*normpdf(xxx,xx(j),de/2);
end
f = ff % интерполяция (3.9)
p = trapz(xxx,f)
p0 = normcdf(max(xxx),m,si)-normcdf(min(xxx),m,si)
bar(xx,h,’--’)
hold on
plot(xxx,normpdf(xxx,m,si),xxx,f,’r’)
 = –0,0495, σ
 = 0,8906. МНК-коОценки м.о. и с.к.о. выборки m
эффициенты b = 0,0263; 0,0222; 0,0988; 0,1156; 0,0598; 0,1656;
0,1628; 0,1506; 0,0984; 0,0623; 0,0200; 0,0071. Нормальная МНКинтерполяция f (x) показана на рис. 3.12. В серединах интервалов
f (xi ) = hi.
h, f
0.4
0.3
0.2
0.1
0
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x
Рис. 3.12. Нормальная МНК-интерполяция, σ = ∆/2
55
 ) попадает в интервал ∆x = (x ,
, σ
Случайная величина x∈N( m
min
xmax) с вероятностью p = 0,9861, величина x ∈ f (x)  – с вероятностью 
p = 0,9840. Гипотеза однородности проверяется как в примере
3.3: в табл. 3.4 приведены значения ni (число попаданий выбороч i, округленное до ближайшего целого функцией
ных значений) и n
ROUND. Статистика (1.15) проверки гипотезы однородности t =
= 0,2484 < χ2(0,95;11) = 19,6751 не дает оснований отвергнуть гипотезу.
Таблица 3.4
№
1
2
3
4
5
6
7
8
9
10
11
12
n
3
4
12
14
10
20
21
19
13
8
3
1

n
3
5
12
13
11
19
21
18
13
8
3
1
При неотрицательных коэффициентах (3.11) интерполирующая
функция (3.9) также неотрицательна. Если в гистограмме имеются
соседние резко различающиеся значения, некоторые коэффициенты
bj нормальной МНК-интерполяции могут стать отрицательными.
Пример 3.8. На рис. 3.13 показаны значения некоторой гистограммы HT = [0,8439 0,0084 0,8439 0,0506 0,0422 0,0844 0,0844
0,0422], расположенной в точках t = –1,75:0.5:1,75 (∆ = 0,5, k = 8).
h,f
0.9
0.8
0.7
0.6
0.5
0.4
0.3
2
0.2
1
0.1
0
-0.1
2
-3
-2
-1
0
1
2
Рис.3.13. Нормальная МНК-интерполяция
56
3
Значения h1 и h2, h2 и h3 различаются на два порядка. Заданы две
матрицы плана с нормальными базисными функциями, имеющими различные с.к.о.: σ1 = 9∆/20 = 0,225; σ2 = ∆/2 = 0,250. МНКкоэффициенты (3.11) b = 0,4824; –0,0770; 0,4836; –0,0142; 0,0214;
0,0422; 0,0423; 0,0202 с суммой s1 = 1,00 при σ1, b = 0,5481; –0,1439;
0,5545; –0,0469; 0,0267; 0,0433; 0,0442; 0,0204 с суммой s2 = 1,05
при σ2. В первом случае (рис.3.13, кривая 1) оценка (3.9) нормирована, во втором (рис.3.13, кривая 2) – не только ненормирована, но
в ней появляются недопустимые отрицательные значения. В обоих случаях среди коэффициентов имеются отрицательные, так что
интерполяция возможна и при наличии отрицательных коэффициентов. Успех интерполяции обуславливается правильным выбором
ширины базисных функций.
Расчет коэффициентов (3.11) предусматривает умножение матрицы, обратной матрице плана, на вектор H. В этом примере
0,6386
-0,0880
0,0119
-0,0016
 –1
Z =
0,0002
-0,0000
0,0000
-0,0000
-0,0880
0,6507
-0,0897
0,0121
-0,0016
0,0002
-0,0000
0,0000
0,0119
-0,0897
0,6509
-0,0897
0,0121
-0,0016
0,0002
-0,0000
-0,0016
0,0121
-0,0897
0,6509
-0,0897
0,0121
-0,0016
0,0002
0,0002
-0,0016
0,0121
-0,0897
0,6509
-0,0897
0,0121
-0,0016
-0,0000
0,0002
-0,0016
0,0121
-0,0897
0,6509
-0,0897
0,0119
0,0000
-0,0000
0,0002
-0,0016
0,0121
-0,0897
0,6507
-0,0880
-0,0000
0,0000
-0,0000
0,0002 ,
-0,0016
0,0119
-0,0880
0,6386
отрицателен коэффициент
b2 = z −1(2,1)h1 + z −1(2,2)h2 + z −1(2,3)h3 + ... =
= –0,0880×0,8439+0,6507×0,0084 –0,0897×0,8439+… = –0,0770
в основном за счет первого и третьего слагаемых. Если бы значение
гистограммы h2 было больше, коэффициент b2 мог стать положительным. Таким образом, условие неотрицательности коэффициентов b формируется структурой строк матрицы Z –1: отношение
η=
z −1(i, j)
z (i, j ± 1)
−1
соседних разнополярных элементов матрицы определяет допустимое отношение соседних значений гистограммы. Так, отношение соседних значений гистограммы h1 / h2 ≈ 100 >> η = |z –1(2,2)/
z –1(2,1)| = 7,39, поэтому коэффициент b2 отрицателен. На практике
вероятность появления резких перепадов соседних значений гис57
тограммы, соответствующей непрерывной (как правило, гладкой)
плотности распределения, невелика.
3.3. Адаптивная нормальная МНК-интерполяция
В подразд. 3.2 базис задавался множеством одинаковых нормальных плотностей с модами в серединах интервалов гистограммы. Следуя принципу адаптации сглаживающих функций к наблюдениям [6], можно с.к.о. нормальных базисных функций связать со значениями гистограммы hi:
σi = α
hmax
.
hi
(3.12)
Тогда с ростом hi ширина базисных функций уменьшается, и плот′
ность оценивается с большим
разрешением, с уменьшением hi (на
′
краях наблюдения следуют с большими
интервалами) разрешение
уменьшается. Оценка плотности аналогично (3.9) записывается
f (x) =
N
∑
j =1
b j z j (x), z j (x) =
 (x − x j ) 2 
exp −
. 2σ 2j 
2πσ 2j

1
(3.13)
Коэффициенты bj рассчитываются методом наименьших квадратов.
Пример 3.9. Генеральная совокупность имеет двумодальную
нормальную плотность распределения
f (x) =
 (x − m1) 2 
 (x − m2 ) 2 
1
+
exp −
exp −

, (3.14)
2
2σ 1  2 2πσ 22
2σ 22 

2 2πσ 12

1
m1 = –5, m2 = 1, σ12 = 1, σ22 = 4. Для одной из реализаций X (N =
= 128) число попаданий ni выборочных значений xi∈f(x) в интервалы (k = 14, ∆ = 1, середины интервалов xi = (i –8,5)∆) приведено в
табл. 3.6.
Таблица 3.6
№
1
2
3
4
5
6
7
9
10
11
12
13
14
ni
1
9
26
18
13
1
6
8
10
11
6
4
1
σi 2,55 0,85 0,50 0,60 0,71 2,55 1,05 0,90 0,81 0,77 1,05 1,27 2,55
С.к.о. (3.12) задано с коэффициентом α = 1/2 (табл. 3.6, рис. 3.14):
σ i = hmax /4hi .
58
h
0.2
0.1
0
-8
-6
-4
-2
0
2
4
6
2
4
6
x
z
0.8
0.6
0.4
0.2
0
-8
-6
-4
-2
0
x
Рис. 3.14. Гистограмма и соответствующие базисные функции
матрицы плана
h, f
а)
2
0.2
0.15
1
0.1
0.05
0
б)
-6
-4
-2
0
2
4
2
4
x
F
1
0.8
0.6
0.4
0.2
0
-6
-4
-2
0
x
Рис. 3.15. Адаптивная нормальная МНК-интерполяция
59
Адаптивная оценка (3.13) показана на рис. 3.15, а (кривая 1), исходная плотность (3.14) – на рис. 3.15, а (кривая 2). Эмпирическая
функция распределения FN(x) выборки X (рис. 3.15, б, штриховая
, полученная интегрированием оценки
диаграмма) и функция F
f (x) по формуле прямоугольников (функция CUMSUM, рис. 3.15,
б, непрерывная кривая) дают статистику (1.5)
(x) = 0,0274 < 1,3581/ N = 0,1200,
D = sup F (x) − F
N
N
так что по критерию Колмогорова гипотеза о том, что выборка X
принадлежит генеральной совокупности с плотностью распределения f (x), не отвергается.
Как видно из рис. 3.15, б, значения эмпирической функции рас(x) отличаются незначительно. Это
пределения FN(x) и функции F
является предпосылкой применения критерия ω2 [2]. Статистика
(1.18), вычисляемая с использованием значений F(xj) функции распределения генеральной совокупности
x −1 
1
1
F(x) = Φ(x + 5) + Φ 
,
2
2  2 
при N = 128 не превышает значения Nω2 = 0,25. В соответствии с
табл. 1.4 проверка по критерию ω2 также не дает оснований отвергнуть гипотезу.
3.4. Моделирование корреляционного обнаружения
Более подробно рассмотреть применение оценок плотностей распределения можно на примере исследования корреляционного обнаружения [16] – построения отношения правдоподобия и рабочей
характеристики обнаружителя [3,16,17]. Стандартная радиотехническая задача – обнаружение импульсного сигнала s(t), маскируемого аддитивным стационарным гауссовым шумом n(t). В терминах математической статистики проверяются гипотезы H0 и H1
о входном сигнале
при гипотезе H 0,
 n(t)
x(t) = 
n(t) + s(t) при гипотезе H1,
n(t) – белый шум; s(t) – прямоугольный сигнал с амплитудой U
длительностью T с известным временем прихода (простые гипотезы). Входной сигнал подвергается преобразованию линейным
фильтром с весовой функцией
h(t) = exp(−αt)sin βt.
60
При этом белый шум окрашивается – приобретает функцию корреляции


α
R(t) = exp(−ατ)  cos βτ + sin βτ , β


(3.15)
сигнал принимает форму



s (t) = U  1 − R(t)  при 0 ≤ t < T,
s(t) =  0


 s0 (t) − s0 (t − T)
при
t ≥ T.
(3.16)
Далее сигнал дискретизируется во времени, и его отсчеты xj подвергаются нелинейному преобразованию (по правилу оценивания
корреляционных моментов)
z1 =
1  2
x1 + x 22 + ... + x n2
n

z3 =

, z2 = 1 (x1x 2 + x 2x 3 + ... + x n −1x n ), (3.17)

n −1

1
(x1x 3 + x 2x 4 + ... + x n −2x n ), ..., zn = x1x n.
n−2
Статистика корреляционного обнаружения
α = Z TG (3.18)
формируется умножением вектора Z на вектор G, минимизирующий дисперсию статистики при гипотезе H0 при условиях M[α|H0] =
= 0, M[α|H1] = 1. Теоретический расчет плотностей распределения
нелинейной статистики (3.18) f0(α) = f(α|H0) и f1(α) = f(α|H1), необходимых для получения отношения правдоподобия
Λ(α) =
и рабочей характеристики
f1(α)
f0 (α)
(3.19)
D = ϕ( F), (3.20)
D, F – вероятности обнаружения и ложной тревоги, весьма сложен
[16]. Практическое решение – моделирование обнаружителя и оценивание плотностей.
Моделирование включает:
– генерирование массива значений гауссова шума [16,17] с корреляционной матрицей, соответствующей функции корреляции
(3.15), и сигнала (3.16);
61
– расчет вектора G [16] с использованием корреляционной матрицы BB оценок (3.17);
– формирование статистик (3.18);
– построение гистограмм статистики при гипотезах H0 и H1.
Гистограммное оценивание плотностей f0(α) и f1(α) позволит пос (α) = f (α)/ f (α) и
троить оценки отношения правдоподобия Λ
1
0
рабочей характеристики
 = ϕ(F
).
D
Дисперсия белого шума задана равной единице, так что при амплитуде U и протяженности сигнала n = 5 отсчетов отношение сигнал-шум на входе d0 = U/σ = U дает эффективное отношение сигнал-шум d = U n. При α = 1/2, β = 2π и интервале дискретизации
σ = 1 задача решается программой
U = 3/4 N = 128
n = 5 rn = [1;-0.6065;0.3679;-0.2231;0.1353] s = U*[1.6065;0.6321;1.2231;0.8647;1.0821]
B = [1.0000 -0.6065 0.3679 -0.2231 0.1353
-0.6065 1.0000 -0.6065 0.3679 -0.2231
0.3679 -0.6065 1.0000 -0.6065 0.3679
-0.2231 0.3679 -0.6065 1.0000 -0.6065
0.1353 -0.2231 0.3679 -0.6065 1.0000]
[u,v] = eig(B) x = randn(n,N);
bx = cov(x’)
[ux,vx] = eig(bx) A = u*v^(1/2)*u’*ux*vx^(-1/2)*ux’ y0 = A*x; for i = 1:n
for j = 1:N
S(i,j) = s(i);
end
end
y1 = y0+S; 62
% амплитуда сигнала
% число отсчетов
% линейное преобразование
% отсчеты функции корреляции (3.15)
% отсчеты сигнала (3.16)
% генератор шума
% корреляционная матрица окрашен
ного шума
% собственные векторы и значения
% собственные векторы и значения
% оператор окрашивания
% окрашенный шум (гипотеза H0)
% шум + сигнал (гипотеза H1)
% расчет вектора G
% корреляционная матрица чисел z (3.17)
BB = [0.7193 -0.6580 0.5269 -0.3898 0.2707
-0.6580 0.6980 -0.5909 0.4451 -0.3052
0.5269 -0.5909 0.6622 -0.5608 0.4177
-0.3898 0.4451 -0.5608 0.7337 -0.6367
0.2707 -0.3052 0.4177 -0.6367 1.0183]
s1 = s(1);s2 = s(2);s3 = s(3);s4 = s(4);s5 = s(5)
% преобразования (3.17) отсчетов сигнала
rs1 = 1/5*(s1^2+s2^2+s3^2+s4^2+s5^2);rs2 = 1/4*(s1*s2+s2*s3+s3*s4+s4*s5);
rs3 = 1/3*(s1*s3+s2*s4+s3*s5);rs4 = 1/2*(s1*s4+s2*s5);rs5 = s1*s5;
rs = [rs1;rs2;rs3;rs4;rs5]
% коэффициенты для расчета вектора G
bb = inv(BB);
a00 = rn’*bb*rn; a01 = rn’*bb*rs; a11 = rs’*bb*rs;
c = a00*a11-a01*a01;
rnn = bb*rn;
rss = bb*rs;
G = (a00*rss-a01*rnn)/c % вектор G
pause
% формирование статистик (3.18)
Y = y0;
z1 = 1/5*sum(Y.^2)
for j = 1:N
z2(j) = 1/4*(Y(1,j)*Y(2,j)+Y(2,j)*Y(3,j)+Y(3,j)*Y(4,j)+Y(4,j)*Y(5,j));
z3(j) = 1/3*(Y(1,j)*Y(3,j)+Y(2,j)*Y(4,j)+Y(3,j)*Y(5,j));
z4(j) = 1/2*(Y(1,j)*Y(4,j)+Y(2,j)*Y(5,j));
z5(j) = Y(1,j)*Y(5,j);
end
z = [z1;z2;z3;z4;z5];
a0 = G’*z; % статистика при гипотезе H0
Y = y1;
z1 = 1/5*sum(Y.^2)
for j = 1:N
z2(j) = 1/4*(Y(1,j)*Y(2,j)+Y(2,j)*Y(3,j)+Y(3,j)*Y(4,j)+Y(4,j)*Y(5,j));
z3(j) = 1/3*(Y(1,j)*Y(3,j)+Y(2,j)*Y(4,j)+Y(3,j)*Y(5,j));
z4(j) = 1/2*(Y(1,j)*Y(4,j)+Y(2,j)*Y(5,j));
z5(j) = Y(1,j)*Y(5,j);
end
z = [z1;z2;z3;z4;z5];
a1 = G’*z; % статистика при гипотезе H1
% гистограммы статистик
63
de = 0.4
xx = -1.2:de:4.4; k = length(xx)
n0 = hist(a0,xx);
h0 = n0/N/de n1 = hist(a1,xx);
h1 = n1/N/de subplot(2,1,1),bar(xx,h0)
hold on
stem(xx,h1)
s0 = sum(h0*de)
s1 = sum(h1*de)
pause
% середины интервалов
% гистограмма при гипотезе H0
% гистограмма при гипотезе H1
При амплитуде U = 3/4 вектор GT = [0,1777;0,5056;0,4823;0,2766;
0,0973]. Гистограммы статистик корреляционного обнаружения
показаны на рис. 3.16, а (штрихами – при гипотезе H1).
% интерполяция и аппроксимация гистограмм
si = de/2
for i = 1:k
for j = 1:k
zz(i,j) = normpdf(xx(j),xx(i),si);
end
end
Z = zz; % матрица плана
xxx = -2.6:de/10:6; % числовая ось в другом масштабе
nn = length(xxx);
b0 = inv(Z)*h0’ % МНК-коэффициенты, гипотеза H0
f = zeros(1,nn);
for j = 1:k
f = f+b0(j)*normpdf(xxx,xx(j),si);
end
f0 = f; % интерполяция, гипотеза H0
sf = trapz(xxx,f0) % проверка нормированности оценки
f0n = f0/sf; % нормировка, гипотеза H0
b1 = inv(Z)*h1’; % МНК – коэффициенты, гипотеза H1
f = zeros(1,nn);
for j = 1:k
f = f+b1(j)*normpdf(xxx,xx(j),si);
end
f1 = f; % интерполяция, гипотеза H1
64
sf = trapz(xxx,f1) f1n = f1/sf; f = zeros(1,nn);
for j = 1:k
f = f+h0(j)*de*normpdf(xxx,xx(j),si);
end
f0a = f; f = zeros(1,nn);
for j = 1:k
f = f+h1(j)*de*normpdf(xxx,xx(j),si);
end
f1a = f; subplot(2,1,2),bar(xx,h1)
hold on
stem(xx,h0)
plot(xxx,f0,’r’,xxx,f1,’k’,xxx,f0a,’k’)
pause
% проверка нормированности оценки
% нормировка, гипотеза H1
% аппроксимация (3.1), гипотеза H0
% аппроксимация (3.1), гипотеза H1
МНК-интерполяция гистограмм показана на рис. 3.16, б, кривые 1
и 3; кривая 2 – аппроксимация (3.1) при гипотезе H0, при гипотезе
а)
h
3.5
3
2.5
2
1.5
1
0.5
0
б)
-1
0
1
2
3
4
3
4
x
h, f
3.5
3
2.5
1
2
1.5
1
3
2
0.5
0
-1
0
1
2
Рис. 3.16. Гистограммы статистики, нормальная интерполяция
и аппроксимация, U = 3/4, N = 128
65
H1 аппроксимирующая и интерполирующая оценки визуально неразличимы.
На рис. 3.17 гистограммные оценки сравниваются с ядерными:
нормальные гистограммные оценки (кривые 1) и оценки Парзена
(2.2) (кривые 2) при амплитуде сигнала U = 1 строятся программой
ami = min(a0);
ama = max(a1);
la = (ama-ami)/N*2;
xp = ami:(ama-ami)/N:ama-(ama-ami)/N;
ff = zeros(1,N);
for j = 1:N
ff = ff+normpdf(xp,a0(j),la);
end
fp0 = ff/N; % оценка Парзена при гипотезе H0
trapz(xp,fp0)
ff = zeros(1,N);
for j = 1:N
ff = ff+normpdf(xp,a1(j),la);
f
7
6
2
5
4
3
2
1
2
1
1
a
0
0
0.5
1
1.5
2
Рис. 3.17. Нормальная гистограммная интерполяция и оценка Парзена,
U = 1, N = 128
66
end
fp1 = ff/N; % оценка Парзена при гипотезе H1
trapz(xp,fp1)
Lp = fp1./fp0; % отношение правдоподобия по Парзену
plot(xp,fp0,xp,fp1,xxx,f0,’r’,xxx,f1,’r’)
pause
Оценки Парзена флюктуируют заметно интенсивнее сглаженных
гистограммных оценок.
На рис. 3.18 (кривая 1) показаны аппроксимирующие гистограммные оценки плотностей статистики для той же выборки. Там
же (кривые 2) изображены нормальные плотности f0(α) = N(m0,σ)
и f1(α) = N(m1,σ) с параметрами m0 = 0, m1 = 1 и дисперсией σ2 =
= 1/d2 = 1/U2n = 0,2 (распределение статистики линейного обнаружения прямоугольного сигнала в белом шуме). Оценки плотностей
корреляционного обнаружения пересекаются в меньшей степени,
что определяет преимущество корреляционного обнаружения над
линейным.
f
3.5
3
2.5
1
2
1.5
1
1
2
2
0.5
alfa
0
-1.5
-1
-0.5
0
m0
0.5
1
m1
1.5
2
2.5
3
Рис. 3.18. Гистограммная аппроксимация плотностей распределения,
U = 1, N = 128
67
Более наглядно сравнение методов обнаружения с использованием отношений правдоподобия (3.19), формируемых программой
n00 = find(n0 = = 0)
n01 = min(find(n0>0))
h0(n00) = [] h1(n00) = []
nop = length(h0)
xop = xx(n01):de:xx(n01)+(nop-1)*de;
Lh = h1./h0; Ln = normpdf(xx,U*sqrt(5),1)./normpdf(xx,0,1);
n00 = find(f0a< = 0.0001)
n01 = min(find(f0a> = 0.0001))
f0a(n00) = [] f1a(n00) = []
nop = length(f0a)
xopp = xxx(n01):de/10:xxx(n01)+(nop-1)*de/10
La = f1a./f0a; subplot(1,2,1),plot(xopp,La,’k’,xop,L,xx,Ln,’r’)
ylim([0 10])
subplot(1,2,2),plot(xopp,La,’k’,xp,Lp)
ylim([0 10*max(la)])
% исключение нулевых значений гистограммы
% отношение правдоподобия по гистограммам
% отношения правдоподобия линей
ного обнаружения
% исключение нулевых значений аппроксимации
% отношение правдоподобия по аппроксимации
Нулевые значения гистограмм и оценок плотностей исключаются,
чтобы избежать деления на ноль. На рис. 3.19 а показаны оценки
отношения правдоподобия корреляционного обнаружения: 1 – отношение La(α) аппроксимирующих плотностей (3.1), 2 – отношение Lh(α) значений гистограмм при гипотезах H1 и H0. Функция
La сглаживает ломаную Lh, являющуюся грубым приближением.
На рис. 3.19, а (кривая 3) – отношение правдоподобия Ln(α) как
отношение нормальных плотностей при линейном обнаружении
сигнала. Крутизна отношения правдоподобия корреляционного
обнаружении выше крутизны отношения правдоподобия линейного обнаружения, что является признаком преимущества первого:
производная отношения правдоподобия (3.19) дает рабочую характеристику (3.20) [16]. На рис. 3.19, б сравниваются оценки отношения правдоподобия La(α) (кривая 1) и Lp(α) – отношение оценок
68
L
L
а)
б)
9
250
8
2
7
200
6
2
3
2
150
5
1
4
100
3
2
50
1
1
0
0
0.5
1
1.5
2
alfa
0
0
0.2 0.4 0.6 0.8
1
alfa
Рис. 3.19. Отношение правдоподобия, U = 1, N = 128
Парзена. Оценки Парзена дают искаженную оценку отношения
правдоподобия по сравнению с гистограммными оценками.
Обнаружитель полностью характеризуется рабочей характеристикой (3.20). Программа расчета рабочих характеристик по оценкам плотностей распределения
Fh = 1-cumsum(h0*de); % к построению РХ по гистограммам
Dh = 1-cumsum(h1*de);
Fa = 1-cumsum(f0a*de/10); % к построению РХ по аппроксимированным гистограммам
Da = 1-cumsum(f1a*de/10);
Fn = 1-cumsum(f0n*de/10); % к построению РХ по интерполированным гистограммам
Dn = 1-cumsum(f1n*de/10);
Fp = 1-cumsum(fp0*dp); % к построению РХ по оценкам Парзена
Dp = 1-cumsum(fp1*dp);
F = 1-normcdf(xxx,0,U/sqrt(5)); % к построению РХ линейного обнаружения
D = 1-normcdf(xxx,U,U/sqrt(5));
На рис. 3.20 показаны рабочие характеристики корреляционного
(построены по оценкам Fa, Da) и линейного обнаружителя (кривые 1 и 2 соответственно) при отношении сигнал-шум на входе
d0 = 1/2 и d0 = 2/3. При таких значениях отношения сигнал-шум
линейное обнаружение эффективнее корреляционного.
69
D
D
1
1
2
0.9
1
0.8
2
0.8
0.7
0.7
0.6
0.6
U = 1/2
0.5
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0.2
0.4
0.6
U = 2/3
0.5
0.4
0
1
0.9
0.8
1
F
0
0
0.2
0.4
0.6
0.8
1
F
Рис. 3.20. Рабочие характеристики
D
D
1
3
5
1
0.9
0.9
0.8
0.8
2
0.7
1
0.7
4
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0
0.1
2
0.2
0.3
0.4
F
0.5
0.2
4
0
0.1
0.2
0.3
0.4
Рис. 3.21. Фрагменты рабочих характеристик, d0 = 3/4
70
F
0.5
При увеличении отношения сигнал-шум до d0 = 3/4 преимущество переходит к корреляционному обнаружению. На рис. 3.21: 1
(ломаная линия) – характеристика, построенная непосредственно
по гистограммам (Fh, Dh); 2 – построена по аппроксимирующим
оценкам (Fa, Da); 3 – построена по нормированным интеполирующим оценкам (Fn, Dn); 4 – рабочая характеристика линейного обнаружения (F, D); 5 – построена по оценкам Парзена. Нормировка
интерполированных оценок плотностей необходима потому, что
их интегрирование дает значения sf = 0,9858 (гипотеза H0) и sf =
= 0,9859 (гипотеза H1). Нормированная интерполяция и аппроксимация гистограмм сглаживают характеристику, построенную
по гистограммам и задают своеобразный доверительный интервал,
внутри которого, по-видимому, и находится истинная рабочая характеристика. Ядерные оценки (Парзена) несмотря на их флюктуации также дают гладкую рабочую характеристику. Последнее
подтверждает, что основной недостаток ядерных оценок плотности
не в их изрезанности, а в необходимости дополнительной аппроксимации.
Преимущество корреляционного обнаружения над линейным
при отношении сигнал-шум на входе d0 ≥ 3/4 – его характерное
свойство [16].
В примере расчета и моделирования корреляционного обнаружения можно выделить два принципиальных момента:
1) сложность теоретического расчета плотностей распределения
нелинейной статистики, которая преодолевается только ее моделированием;
2) использование гистограмм статистики для построения сглаженных характеристик обнаружителя требует малых размеров
гистограммных интервалов, что, в свою очередь, приводит к необходимости увеличения размеров выборок. Использование гистограммных оценок плотностей распределения статистики позволяет получить достаточно гладкие оценки характеристик обнаружителя при ограниченном размере выборки и не требует чрезмерного
уменьшения гистограммных интервалов.
71
4. Оценивание двумерных плотностей
4.1. Двумерная гистограмма
На рис. 4.1 показана двумерная реализация X, состоящая из N =
= 200 точек с независимыми случайными координатами x1∈χ2 (4),
x2∈N(0,1).
N = 200
x1 = randn(1,N); x2 = randn(1,N); x3 = randn(1,N); x4 = randn(1,N);
x = x1.^2+x2.^2+x3.^2+x4.^2;
y = randn(1,N);
plot(x,y,’k*’)
pause
Пары координат xi1, xi2 извлечены из генеральной совокупности с
совместной плотностью распределения (рис. 4.2)
f (x1, x 2 ) =
 x + x 22 
x1
exp  − 1
. 2 
4 2π

(4.1)
Сечения плотности (4.1) – условные одномерные плотности распределения f(x1|x2), f(x2|x1) – в точке ее максимума (x1 = 2, x2 = 0)
показаны на рис. 4.3: а – плотность χ2-распределения с четырьмя
степенями свободы; б – стандартная нормальная плотность.
x2
2
1.5
1
0.5
0
x1
-0.5
-1
-1.5
-2
-2.5
0
2
4
6
8
Рис. 4.1. Двумерная реализация
72
10
12
0.08
0.06
0.04
0.02
0
4
3
2
1
0
15
-1
-2
10
-3
5
-4 0
Рис. 4.2. Двумерная плотность
а)
f(x1|x2)
б)
f(x2|x1)
0.07
0.07
0.06
0.06
0.05
0.05
0.04
0.04
0.03
0.03
0.02
0.02
0.01
0.01
0
0
5
10
x1
15
0
x2
-4
-2
0
2
4
Рис. 4.3. Условные плотности распределения координат
73
Пример 4.1. Двумерная выборка извлечена из генеральной совокупности с плотностью распределения (рис. 4.4)
f (x1, x 2 ) =
 x 2 + x 22 
1
exp  − 1
. 2π
2


(4.2)
0.2
0.15
0.1
0.05
0
4
3
2
1
0
-1
-2
-3
-4
-4
-3
-2
-1
0
1
2
3
4
Рис. 4.4. Двумерная нормальная плотность
Двумерная гистограмма выборки может быть построена следующим образом. Если задать ширину интервала гистограммы ∆ и
сдвиг d, то номер интервала, в который попадает случайное наблюдение xi, равен (округление до ближайшего целого)
n = round(x j / ∆ + d). (4.3)
При ∆ = 0,5, d = 9 число x = –1,8 попадает в пятый интервал, число
x = 0 – в девятый, число x = 1,8 – в тринадцатый интервал. Таким
образом, (4.3) – преобразование выборочных значений xi1 в номер
строки, значений xi2 – в номер столбца гистограммы, а гистограмма записывается в виде матрицы. Относительная сумма числа попаданий в каждую ячейку матрицы задает значения hij.
Нумерация интервалов начинается с n = 1, поэтому при ∆ = 0,5,
d = 9 соотношение (4.3)
74
1 = round(2x + 9)
допускает минимальное значение xmin = –4.25. Вероятность того,
что стандартное нормальное число x < –4,25, равна p ≈ 10 –5, следовательно, в генерируемых выборках размерами сотни-тысячи значения x < xmin могут появляться редко.
Расчет и изображение двумерной гистограммы (рис. 4.5) выполняется MATLAB-программой
N = 1024
d = 1/2 dd = 9 x1 = randn(1,N); x2 = randn(1,N);
X = round(x1/d+dd); Y = round(x2/d+dd); xma = max(X); yma = max(Y);
nn = max(xma,yma) h = zeros(nn);
for i = 1:N
j = X(i);
k = Y(i);
% размер интервала гистограммы
% сдвиг
% номер интервала, в который попадает x1
% номер интервала, в который попадает x2
% число интервалов гистограммы
h
0.2
0.15
0.1
0.05
0
4
2
y
2
0
x
4
0
-2
-4
-2
-4
Рис. 4.5. Двумерная гистограмма
75
h(j,k) = h(j,k)+1; % накапливающееся число попаданий в элемент гистограммы
end
h = h/N/d^2 % двумерная гистограмма
[XX,YY] = meshgrid(((1:nn)-dd)*d,((1:nn)-dd)*d)
meshc(XX,YY,h)
colormap(colorcube)
pause
Как в одномерном случае (см. рис. 1.2), двумерная гистограмма есть ступенчатая оценка совместной плотности, показанной на
рис. 4.4. Ее изрезанность объясняется малым размером выборки
N = 1024 (32 значения по одной координате).
В табл. 4.1 приведены значения
hij =
nij
(4.4)
N∆ 2
нормированной гистограммы (строки и столбцы соответствуют рис.
4.5).
Таблица 4.1
№
4
5
6
7
8
9
10
11
12
3
4
5
6
7
8
9
10
11
12
13
14
15
0
0
0
0,0039
0,0078
0,0039
0,0039
0
0,0078
0
0,0039
0
0
0
0
0,0039
0,0078
0,0078
0,0039
0,0156
0,0273
0,0195
0,0039
0,0078
0
0
0
0,0078
0,0078
0,0195
0,0586
0,0352
0,0469
0,0391
0,0352
0,0234
0,0078
0,0078
0
0
0,0039
0,0078
0,0430
0,0469
0,0703
0,0703
0,0664
0,0508
0,0430
0,0195
0,0039
0
0
0
0,0195
0,0234
0,1016
0,1016
0,1875
0,0898
0,0859
0,0664
0,0234
0,0156
0
0
0,0039
0,0352
0,0625
0,1094
0,1172
0,1250
0,1758
0,0898
0,0469
0,0156
0,0117
0,0039
0,0078
0,0078
0,0117
0,0508
0,0781
0,0703
0,1523
0,0977
0,0742
0,0625
0,0313
0,0039
0
0,0039
0,0117
0,0078
0,0391
0,0742
0,0703
0,0977
0,1211
0,0664
0,0234
0,0195
0
0
0
0,0039
0,0156
0,0195
0,0430
0,0664
0,0586
0,0547
0,0313
0,0156
0,0156
0
0,0039
На рис. 4.6 и 4.7 показаны несколько вертикальных и горизонтальных сечений двумерной гистограммы (функции IMPROFILE и
IMCONTOUR [18]), показывающих ее изрезанность. Изображения
на рис. 4.6 и 4.7 формируются программой
xx = 1:xma;
for j = round(nn/2):nn
yy = ones(1,xma)*j;
improfile(h,xx,yy)
hold on
76
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
Y
X
Рис. 4.6. Сечения гистограммы вертикальными плоскостями
Рис. 4.7. Сечения гистограммы горизонтальными плоскостями
77
2
4
6
8
10
12
14
2
4
6
8
10
12
14
Рис. 4.8. Изображение матрицы-гистограммы
grid on
end
pause
levels = 0:0.01:0.25;
imcontour(h,levels)
pause
Двумерную гистограмму-матрицу можно показать ее полутоновым изображением (функция IMSAGESC [18]). На рис. 4.8 показано «негативное» изображение (imagesc(-h+1)), в котором темные
′
пикселы соответствуют большим
значениям гистограммы.
4.2. Нормальная двумерная гистограммная аппроксимация
Подобно рассмотренной в подразд. 2.2 сумме (2.2) нормальных
плотностей можно аппроксимировать двумерную гистограмму
взвешенной суммой двумерных нормальных плотностей с модами
в узлах гистограммы. Оценка двумерной плотности при этом
78
k k
 (x − x i ) 2 + (y − y j ) 2 
∆2
f (x, y) =
h
exp
−
. ij
2πσ 2 i =1 j =1
2σ 2


∑∑
(4.5)
В сумме (4.5) использованы более простые обозначения: x = x1,
y = x2. Аппроксимирующая плотность (4.5) нормирована: с учетом
(4.4) и k2 = N
∞ ∞
∫∫
f (x, y)dxdy = ∆ 2
k
k
nij
∑∑ N∆ 2
i =1 j =1
−∞ −∞
=
k
k
1
nij = 1.
N i =1 j =1
∑∑
Пример 4.2. Программа
N = 256
x1 = randn(1,N)+1; x2 = randn(1,N)-1; x3 = rand(1,N)*10-5;
r1 = -0.7; r2 = 0.7;
y1 = randn(1,N)*sqrt(1-r1^2)+x1*r1+1; % условная плотность
y2 = randn(1,N)*sqrt(1-r2^2)+x2*r2+1; % условная плотность
y3 = rand(1,N)*10-5;
x = cat(2,x1,x2,x3) % конкатенция массивов
y = cat(2,y1,y2,y3)
plot(x,y,’k*’)
pause
генерирует три группы по N = 256 точек со случайными координатами (рис. 4.9). Точки первой группы имеют координаты x1∈N(1,1),
y1∈N(1,1) с коэффициентом корреляции r1 = –0,7, второй групy
5
4
3
2
1
x
0
-1
-2
-3
-4
-5
-5
-4
-3
-2
-1
0
1
2
3
4
5
Рис. 4.9. Случайные точки
79
пы x2∈N( –1,1), y2∈N(1,1) с коэффициентом корреляции r1 = 0,7.
В третьей группе независимые точки с равномерным распределением координат на интервале ( –5; 5). Точки первых групп расположены относительно компактно, точки третьей группы имитируют
импульсную помеху. Гистограмма (рис. 4.10) и ее контурный график (рис. 4.11) строятся программой
d = 1/2 dd = 11 nx = round(x/d+dd); ny = round(y/d+dd); xma = max(nx)
yma = max(ny)
nn = max(xma,yma) n = zeros(nn);
for i = 1:3*N
j = nx(i);
k = ny(i);
n(j,k) = n(j,k)+1; end
h = n/3/N/d^2 % размер интервала гистограммы
% сдвиг
% номер интервала, в который попадает координата x
% номер интервала, в который попадает координата y
% число интервалов гистограммы
% число попаданий в элемент гистограммы
% гистограмма
h
0.12
0.1
0.08
0.06
0.04
0.02
0
5
5
0
0
x
-5
-5
Рис. 4.10. Гистограмма
80
y
sum(sum(h))
[XX,YY] = meshgrid(((1:nn)-dd)*d,((1:nn)-dd)*d);
mesh(XX,YY,h)
colormap(gray)
pause
levels = 0.002:0.002:0.1;
contour(YY,XX,h,levels,’k’) % контурный график гистограммы
pause
Точки с равномерным распределением координат образуют хаотичный фон. Существование провалов в гистограмме иллюстрируется
контурным графиком: белые пятна внутри темной области соответствуют резким уменьшениям значений гистограммы.
y
5
4
3
2
1
0
x
-1
-2
-3
-4
-5
-5
-4
-3
-2
-1
0
1
2
3
4
5
Рис. 4.11. Контурный график гистограммы
Оценка гистограммы (4.5) и ее изображение выполняются программой
si = d % с.к.о. базисных функций
de = d/4
XXX = ((1:nn)-dd)*d; YYY = ((1:nn)-dd)*d; % центры интервалов гистограммы
xx = -3:de:3; yy = -3:de:3; % оси
nnn = length(xx)
Z = zeros(nnn);
81
for i = 1:nn
for j = 1:nn
for ii = 1:nnn
for jj = 1:nnn
Z(ii,jj) = Z(ii,jj)+h(i,j)*(1/2/pi/si^2*exp(-((xx(ii)-XXX(i)).^2+(yy(jj)-...
YYY(j)).^2)/2/si^2))*d^2; % оценка (4.5)
end
end
end
end
mesh(yy,xx,Z) % рис. 4.12
colormap(gray)
grid on
pause
A = -Z*4+1;
subplot(1,2,1),imagesc(xx,yy,A) % полутоновое изображение
colormap(gray)
levels = 0.02:0.01:0.1;
subplot(1,2,2),contour(xx,yy,Z,levels,’k’) % контурное изображение
На рис.4.12 и 4.13 показаны оценка (4.5), ее полутоновое изображение и контуры горизонтальных сечений (функции IMAGESC и
CONTOUR разворачивают изображения на 90º). Значение с.к.о. базисных функций завышено – равно размеру интервала гистограммы
(σ = d). Поэтому сглаженная гистограмма становится значительно
шире исходной и уменьшается по высоте. Хаотичный фон становится практически незаметным. Просматривается наличие двух групп
компактных частично перекрывающихся точек. Полутоновое и
контурное изображения характеризуют конфигурацию оценки.
На рис. 4.14 показана оценка (4.5) при значении с.к.о., близком
к оптимальному: σ = d/2. Компактные группы точек просматриваются более отчетливо, на контурном изображении (рис. 4.15) проявляются три максимума оценки f.
Применение заниженного значения σ = d/3 уменьшает эффект
сглаживания гистограммы и приводит к фрагментации изображений оценки (рис. 4.16, 4.17). При этом теряется информация о наличии групп компактных точек.
Сумму (4.5) можно использовать для нормальной аппроксимации произвольных поверхностей с приблизительным сохранением
их объема: численное интегрирование по формуле прямоугольников на координатных сетках разного масштаба даст несколько отличающиеся результаты.
82
f
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
3
2
1
0 x
-1
-2
-2
-3 -3
-1
0
3
2
y 1
Рис. 4.12. Нормальная интерполяция гистограммы, σ = d
x
-2
2
-1
1
0
y
0
1
-1
2
-2
3
-2
0
x
3
-3
2
-3
y
-3
-2
-1
0
1
2
3
Рис. 4.13. Полутоновое и контурное изображения оценки
83
f
0.1
0.08
0.06
0.04
0.02
0
3
2
1
0
x -1
-2
-3
-2
-3
-1
3
2
1
0 y
Рис. 4.14. Нормальная интерполяция гистограммы, σ = d/2
x
-3
3
-2
2
-1
1
0
y
0
1
-1
2
-2
3
-2
0
2
x
-3
-3
y
-2
-1
0
1
Рис. 4.15. Полутоновое и контурное изображения
84
2
3
f
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
3
2
1
x
0
-1
-2
-3
-3
0
-1
-2
y
3
2
1
Рис. 4.16. Нормальная интерполяция гистограммы, σ = d/3
x
-3
3
-2
2
-1
1
0
y
0
1
-1
2
-2
3
-2
0
2
x
-3
y
-3
-2
-1
0
1
2
3
Рис. 4.17. Полутоновое и контурное изображения
85
Пример 4.3. Функция (рис. 4.18)
π
π
f (x, y) = sin 2 x sin 2 y, 8
8
(4.6)
заданная на координатной сетке 16×16 с интервалом дискретизации ∆ = 1, маскируется помехой с равномерной плотностью в интервале (0; 1/2). Ее преобразование (4.5) на сетке 64×64 с интервалом
∆ = 1/4 и с.к.о. базисных функций σ = ∆/2 показано на рис. 4.19.
f
1.5
1
0.5
0
15
x
10
10
5
0
y
15
5
0
Рис. 4.18. Сглаживаемая функция
f
1.4
1.2
1
0.8
0.6
0.4
0.2
0
15
x
10
5
0
0
5
10
Рис. 4.19. Сглаженная функция
86
y
15
 =
Исходный объем V = 127,6126, объем сглаженной функции V
= 119,5584. Сглаживание фильтрует помеху и приближает результат к виду (4.6), показанному на рис. 4.20.
f
1.5
1
0.5
0
15
x
10
10
5
y
15
5
0
0
Рис. 4.20. Функция (4.6)
Пример 4.4. Применение преобразования (4.5) для фильтрации
изображения. На рис. 4.21, а изображен черный квадрат (матрица единиц размером 8 × 8) на белом поле (матрица нулей размером
12 × 12). Квадрат искажается помехой (функция RAND), так что
матрица изображения принимает вид (рис. 4.21, б)
000
0
0
0
0
0
0
0
00
000
0
0
0
0
0
0
0
00
0 0 0,050 0,769 0,393 0,514 0,109 0,238 0,544 0,982 0 0
0 0 0,177 0,556 0,385 0,208 0,078 0,262 0,824 0,594 0 0
0 0 0,065 0,083 0,590 0,106 0,942 0,647 0,187 0,990 0 0
0 0 0,861 0,797 0,801 0,396 0,728 0,801 0,985 0,253 0 0
0 0 0,555 0,068 0,534 0,581 0,154 0,475 0,797 0,328 0 0
0 0 0,162 0,980 0,319 0,621 0,168 0,497 0,291 0,571 0 0
0 0 0,695 0,810 0,807 0,318 0,697 0,458 0,849 0,302 0 0
0 0 0,622 0,140 0,146 0,406 0,503 0,100 0,178 0,355 0 0
000
0
0
0
0
0
0
0
00
000
0
0
0
0
0
0
0
00
87
а)
б)
2
2
4
4
6
6
8
8
10
10
12
12
2
4
6
8
10
12
в) 0
2
2
2
4
4
6
6
8
8
10
10
2
4
6
8
6
8
10
12
0
г)
0
4
10
0
2
4
6
8
10
Рис. 4.21.
Преобразование (4.5) на полной координатной сетке с интервалом
∆ = 1/5 (56 ×56 пикселов) и с.к.о. σ = 1 дает расфокусированное
изображение квадрата (рис. 4.21, в). Уменьшение с.к.о. увеличивало бы резкость изображения, но его интенсивность приближалась
бы к показанной на рис. 4.21, б. Так как положение квадрата известно, преобразование (4.5) можно реализовать в пределах границ
квадрата и увеличить с.к.о. для сглаживания интенсивности. На
рис. 4.21, г показан результат такого ограниченного преобразования с ∆ = 1/5 и σ = 5 (36 ×36 пикселов): визуально преобразованное
изображение не отличается от исходного. Конкретно эффект сглаживания показан на рис. 4.22: при средней интенсивности преобразованного изображения B = 1 максимальная интенсивность (один
из центральных пикселов) Bmax = 1,152, центральная часть размером 30 ×30 пикселов имеет интенсивность B ≥ 0,9.
Оценку двумерной гистограммы (4.5) можно распространить на
m-мерный случай. Многомерная гистограмма задается m-мерным
массивом чисел
88
B
1.2
1
0.8
0.6
0.4
0.2
0
0
2
4
6
8
10
12
Рис. 4.22. Интенсивность в центральном горизонтальном сечении
hη =
nη
N∆ m
,
где η – m-мерный номер гиперкуба со стороной ∆. Нормальная многомерная аппроксимация плотности распределения с одинаковым
по всем координатам с.к.о. базисных m-мерных плотностей записывается суммой
f (X ) = ∆ m h Ν(X − µ ; σ),
∑
η
η
η
η
η
в которой μη – центры гиперкубов; Xη – m-мерная координатная
сетка. Аппроксимирующая сумма нормирована:
∞
∞
... f (X η )dX η =
∫ ∫
−∞
−∞
n
∑ Nη = 1.
η
Например, трехмерная аппроксимирующая функция имеет вид
f (x, y, z) =
 (x − x i ) 2 + (y − y j ) 2 + (z − zk ) 2 
hijk exp −
.
2σ 2
i =1 j =1 k =1


kx ky kz
∆3
3
(2πσ 2 ) 2
∑∑∑
В общем случае базисные функции могут задаваться в виде различных плотностей распределения jη(Xη), гиперкубы трансформи89
роваться в m-мерные параллелепипеды с разными сторонами ∆η,
ширина δη базисных функций также может быть различной. Тогда
оценка запишется подобно многомерной оценке Парзена [7]
f (X η ) = ∆ m
∑ hηϕ η(X η − µ η;δ η).
η
Практическая ценность подобных громоздких или обобщающих
формул невелика. При увеличении размерности массивов время на
их обработку возрастает в геометрической прогрессии, так что необходимы специальные меры преодоления этих трудностей. Надо
также иметь в виду невозможность наглядного графического отображения многомерных массивов, а их проекции на координатные
плоскости могут оказаться малоинформативными.
Заключение
Непараметрические гистограммные оценки гладких одномерных плотностей распределения методами нормальной аппроксимации и нормальной МНК-интерполяции удовлетворительны с
позиций критериев согласия и однородности. Гистограммное оценивание применимо при размерах выборки в десятки-сотни наблюдений. В отличие от ядерных оценок Парзена или n ближайших соседей нормальная гистограммная аппроксимация и интерполяция
дают результаты, не нуждающиеся в дополнительной аппроксимации.
Метод нормальной аппроксимации распространяется на двумерные плотности, что может оказаться полезным при обработке
изображений. Его применение к оцениванию многомерных плотностей требует решения проблемы быстрого перебора элементов
многомерного пространства.
90
Библиографический список
1. Кибзун А. И. и др. Теория вероятностей и математическая статистика. М.: Физматлит, 2007. 232 с.
2. Горяинов В. В. и др. Математическая статистика / Под ред.
И. С. Зарубина, А. П. Крищенко М.: МГТУ им. Н. Э. Баумана, 2002.
424 с.
3. Ивченко Г. И., Медведев Ю. И. Математическая статистика.
М.: Высшая школа, 1985. 248 с.
4. Кендалл М, Стьюарт А. Статистические выводы и связи. М.:
Наука, 1973. 899 с.
5. Холлендер М., Вулф Д. Непараметрические методы статистики. М.: Финансы и статистика, 1983. 518 с.
6. Дуда Р., Харт П. Распознавание образов и анализ сцен. М.:
Мир, 1976. 511 с.
7. Фомин Я. А., Тарловский Г. Р. Статистическая теория распознавания образов. М.: Радио и связь, 1986. 264 с.
8. Бендат Дж., Пирсол А. Прикладной анализ случайных данных. М.: Мир, 1989. 540 с.
9. Воробьев С. Н., Осипов Л. А. Регрессионный анализ. СПб.:
ГУАП, 2000. 66 с.
10. Ануфриев И. и др. MATLAB 7. СПб.: БХВ – Петербург, 2005.
1104 с.
11. Чистяков В. П. Курс теории вероятностей. М.: Наука, 1982.
256 с.
12. Королюк В. С. и др. Справочник по теории вероятностей и
математической статистике. М.: Наука, 1985. 640 с.
13. Дьяконов В., Круглов В. Математические пакеты расширения MATLAB. СПб.: Питер, 2001. 480 с.
14. Бронштейн И. Н., Семендяев К. А. Справочник по математике. М.: Наука, 1986. 544 с.
15. Иванов В. А. и др. Математические основы теории автоматического регулирования. М.: Высшая школа, 1971. 808 с.
16. Воробьев С. Н. Эффективное обнаружение детерминированных сигналов. СПб.: ГУАП, 2003. 139 с.
17. Воробьев С. Н., Осипов Л. А. Линейные системы. Расчет и моделирование. СПб.: ГУАП, 2004. 122 с.
18. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов
и изображений. СПб.: Питер, 2002. 608 с.
91
Учебное издание
Воробьев Станислав Николаевич
Гирина Наталия Владимировна
Осипов Леонид Андроникович
ГИСТОГРАММНОЕ ОЦЕНИВАНИЕ
ПЛОТНОСТИ РАСПРЕДЕЛЕНИЯ
Учебное пособие
Редактор А. В. Подчепаева
Верстальщик С. Б. Мацапура
Сдано в набор 09.09.08. Подписано к печати 22.10.08.
Формат 60×84 1/16. Бумага офсетная. Печать офсетная.
Усл.-печ. л. 5,35. Уч.-изд. л. 5,75. Тираж 200 экз. Заказ №
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
Документ
Категория
Без категории
Просмотров
0
Размер файла
4 444 Кб
Теги
vorobev, 0d950131da
1/--страниц
Пожаловаться на содержимое документа