Reconstruction of probability density distribution by histogram method in soil science and ecology

Abstract

Different methods of probability density reconstruction are investigated in that article. We shown that “histogram method with equal intervals of groupment” which is the most frequently using in soil science and ecology for the solving of that task, is really ill-suited. The “histogram method with intervals of equal probability” is recommended among the simplest methods for solving of that task.

Full Text

Введение Экспериментальные данные представляют собой реализации случайных величин или случайных процессов. Из теории вероятностей известно, что исчерпывающее описание случайной величины дается ее законом распределения вероятностей - правилом, позволяющим определять вероятность попадания этой величины в любую заданную область ее значений. Закон распределения случайной величины (распределение) может быть задан с помощью любой из двух взаимно однозначно связанных между собой функций: функции распределения и плотности вероятности [Костылев с соавт., 1991: с. 56-57]. Введём необходимую для дальнейших рассуждений терминологию. Пусть в эксперименте исследуется случайная величина у, причем измеренные значения оказались равными у1, …, уn; измерения независимы между собой. Подмножество из n элементов у1, …, уn из генеральной совокупности называется выборкой объема n. Для получения эмпирического распределения случайной величины у имеющийся набор значений у1, …, уn (выборку объема n) представляют в виде гистограмм: для непрерывной случайной величины у производится разбиение отрезка значений этой величины [C1 Ck+1] на k интервалов [C1 C2], [C2 C3], …, [Ck Ck+1]; здесь C1 равно наименьшему значению в выборке, а Ck+1 - наибольшему, т.е. C1 = min (уi), Ck+1 = max (уi) i=1,…,n i=1,…,n (отрезок [C1 Ck+1] называется «интервалом выборки», а отрезки [C1 C2], …, [Ck Ck+1] - «интервалами группировки» [Костылев с соавт., 1991: с. 78]). Для простоты некоторых дальнейших обозначений мы будем считать, что выборка у1, …, уn представлена в виде вариационного ряда, т.е. упорядочена по возрастанию элементов. Тогда C1 = у1, Ck+1 = уn. Подсчитывается число элементов исходного массива у1, …, уn, лежащих в каждом интервале группировки. Обозначим эти числа fi (i = 1, …, k). При этом предполагается, что каждый элемент (yj), попавший в i-ый интервал, совпадает с «усредненным» значением интервала, то есть равняется этому значению: Či = (Ci + Ci+1)/2. Затем определяются статистические частоты рi путем деления чисел fi на объем выборки: рi = fi/n. Совокупность значений Či и соответствующих им частот рi называют гистограммой [Живописцев и Иванов, 1993: с. 49-50] (иногда вместо рi по оси ординат откладывают просто количества элементов fi). Типичная гистограмма наблюдений (Či, fi) представлена на рис. 1. Глаголев%20и%20Сабреков,%20Рис Рис. 1. Типичная гистограмма наблюдений. При подготовке группированных данных большую значимость приобретает вопрос об обоснованном разбиении интервала выборки на интервалы группировки. В литературе имеется множество рекомендаций по выбору числа интервалов группировки (см., например, [Костылев с соавт., 1991: с. 78-79] и ссылки там, а также [Косарев, 2003; 2008: с. 31-32]). Однако при этом оказывается, что для одной и той же выборки разные методы выбора числа интервалов могут приводить к качественно разным распределениям. В данной работе мы преследовали следующие цели: 1) проанализировать те особенности, которые имеет построение плотности распределения в почвоведении и экологии; 2) предложить метод восстановления (по экспериментальным данным) плотности распределения, лишенный неоднозначности (свойственной методам построения гистограмм), связанной с выбором интервалов группировки. Теория Классические способы разбиения выборки на равные интервалы группировки Здесь мы упомянем лишь некоторые подходы, которые часто рекомендуются для практического использования [Костылев с соавт., 1991: с. 78; Живописцев и Иванов, 1993: с. 50; Косарев, 2008: с. 32]): 1. Формула Старджеса: k = 1 + log2(n). 2. Формула Брукса-Каррузера: k = 5·lg(n). 3. Формула Хайнхольда-Гаеде: k = n½. 4. 1-ая формула Костылева с соавт.: k = 4·lg(n). 5. 2-ая формула Костылева с соавт.: k = 5·lg(n/10). 6. Формула Живописцева-Иванова: k = 10·ln(n). 7. Формула Косарева: k = n1/3. Заметим, что уже на этом этапе возникает неоднозначность. Даже не говоря о том, что разные формулы дают различные результаты, то, что по своему смыслу количество интервалов группировки должно быть целым (а это условие не обеспечивает ни одна из вышеприведенных формул), ставит нас перед выбором процедуры округления до целого. Как правильней провести эту процедуру? Очевидно, что бóльшая однозначность достигается при округлении либо всегда в сторону меньшего целого числа, либо всегда в сторону большего. Однако при округлении в сторону большего целого мы будем получать неудовлетворительные результаты при малых k, поэтому остается рекомендовать округление в сторону меньшего целого. К сожалению, это не сильно улучшает ситуацию. В табл. 1 для разных n приведены k, рассчитанные по различным формулам (с округлением в сторону меньшего целого). Заметим, что сейчас мы не рассматриваем формулу Живописцева-Иванова в силу ее кажущейся очевидной абсурдности: действительно, например, выборку, содержащую лишь 3 элемента, согласно этой формуле следует разделить на… k = 10·ln(4) = 10.986 » 10¸11 интервалов (однако ниже мы увидим, что в некоторых ситуациях именно эта формула дает лучшие результаты). Таблица 1 Количество интервалов группировки (k) для различных объемов выборок (n) k Расчет по формуле Старджеса Брукса-Каррузера Хайнхольда-Гаеде Костылева с соавт. 1-ой 2-ой БО СО БО СО БО СО БО СО БО СО 1 1 1 0 1 1 1 0 1 -5 1 2 2 1 1 1 1 1 1 1 -3 1 4 3 1 3 1 2 1 2 1 -1 1 8 4 2 4 2 2 2 3 2 0 2 16 5 3 6 3 4 3 4 3 1 3 32 6 5 7 5 5 5 6 5 2 3 64 7 5 9 5 8 5 7 5 4 3 128 8 7 10 7 11 7 8 7 5 5 256 9 9 12 11 16 11 9 9 7 7 512 10 9 13 13 22 15 10 9 8 7 1024 11 11 15 15 32 19 12 11 10 9 Примечание: БО - без ограничений; СО - с ограничением (подробности см. далее в тексте). Из табл. 1 видно, что с увеличением объема выборки разрыв между оптимальными значениями k, рекомендуемыми различными формулами, увеличивается (заметим, что пока мы обсуждаем формулы «в чистом виде», без каких-либо ограничений, т.е. мы обсуждаем столбцы «БО» соответствующих формул в табл. 1). В частности, k может быть от 10 до 32 для n ~ 103. Впрочем, даже для выборок малых объемов вариации количества интервалов также существенны. Например, согласно 2-ой формуле Костылева с соавт. выборка объема n = 16 еще слишком мала, чтобы ее можно было представлять в виде гистограммы. Действительно, при n = 16 по 2-ой формуле Костылева с соавт. имеем k = 1, т.е. должен быть лишь один интервал группировки, и смысл построения гистограммы теряется. А по формуле Брукса-Каррузера получается, что гистограмму не только можно строить, но, более того, нужно рассортировать выборку по достаточно большому числу интервалов - по шести! Кстати, обратим внимание на поведение формул при малых объемах выборок. Достаточно очевидно, что довольно бессмысленно строить гистограмму если выборка состоит лишь из двух элементов, однако, например, формула Старджеса (которая, пожалуй, применяется на практике чаще всего), все-таки, рекомендует это делать, разбив интервал выборки на 2 интервала группировки. Не лучше и абсурдная рекомендация разбить выборку из 4 элементов по трем интервалам группировки, даваемая формулой Брукса-Каррузера. Лишь 2-ая формула Костылева с соавт. не рекомендует строить гистограммы при малых объемах выборки, но она в этом «перегибает палку». Представляется достаточно очевидным, что хоть на два-то интервала группировки можно было бы разбить выборку, содержащую 16 элементов. Указанные трудности были давно осознаны, в результате чего на практике стали использовать различные ограничения. Так, согласно [Костылев с соавт., 1991: с. 79], часто полагают, что, во-первых, 0.55·n0.4 ≤ k ≤ 1.25·n0.4 (1) и, во-вторых, k рекомендуют выбирать нечетным. Сразу заметим, что эти два правила иногда противоречат друг другу. Так, например, при n = 8 0.55·n0.4 » 1.3, а 1.25·n0.4 » 2.9, т.е. первому правилу удовлетворяет единственное целое значение k = 2, но оно не удовлетворяет второму правилу. Такое же противоречие мы имеем и при n = 5, 6 или 7. Учитывая, что второе правило носит лишь рекомендательный характер, при возникновении подобного противоречия мы будем следовать только первому правилу. Заметим, что ограничения, предлагаемые различными авторами, очень сильно различаются. Так, Ф.А. Живописцев и В.А. Иванов [1993: c. 50] предлагают примитивное ограничение, не зависящее от объема выборки (ср. с №1): 5 ≤ k ≤ 30 (1а) В столбцах «СО» табл. 1 приведено число интервалов группировки, рассчитанное по каждой формуле с учетом (№1) и требования нечетности. Из табл. видно, что, действительно, теперь достигнута гораздо большая однозначность, хотя и не полная. Например, при n = 32¸64, все формулы кроме последней дают оптимальное значение k = 5, а последняя формула рекомендует k = 3. На это можно было бы не обращать внимание, однако при относительно больших n разнобой возникает и между другими формулами (см. последние три строки в табл. 1). Причем для n ~ 103 оптимальные количества интервалов, рекомендуемые формулами Старджеса и Хайнхольда-Гаеде различаются почти в два раза (мы специально сравниваем здесь только эти две хорошо известные и часто применяемые формулы, если же обратить внимание на более экзотическую 2-ую формулу Костылева с соавт., то мы увидим, что различие может достигать более двух раз). Пока мы занимались только простейшим классом формул разбиения интервала выборки на интервалы группировки: все рассмотренные выше формулы подразумевали, что k - это количество равных интервалов группировки, т.е. если k было найдено по одной из вышеприведенных формул, то каждый интервал [Cj Cj+1] строился таким образом, что Cj+1 = Cj + DС, где DС = (Ck+1 - C1)/k = const. Не говоря уже о том, что существенная неоднозначность разбиения появится тогда, когда мы перейдем к методам неравных интервалов группировки, даже сейчас, все еще оставаясь в классе формул равных интервалов, мы вскроем дополнительный источник неоднозначности разбиения. Способы разбиения, зависящие от типа распределения Понятно, что представление в виде гистограммы (функции плотности вероятности на системе интервалов группировки) с формально-математической точки зрения является кусочно-постоянной аппроксимацией (т.е. аппроксимацией кусочными полиномами 0-го порядка) этой функции на сетке с узлами C1, C2, …, Ck+1 (кстати, заметим, что почти столь же часто используемое, как и гистограмма, представление плотности вероятности в виде «полигона частот» - см., например, [Дмитриев, 1995: c. 45] - представляет собой кусочно-линейную аппроксимацию, т.е. аппроксимацию кусочными полиномами 1-го порядка). Из курса вычислительной математики (см., например, [Петров и Лобанов, 2006: с. 138]) известно, что на равномерной сетке, заданной на отрезке [a, b], абсолютная погрешность интерполяции функции f(С) одной переменной С (под этой погрешностью мы будем понимать абсолютное значение остаточного члена и обозначать ее RN) удовлетворяет оценке: RN ~ DСN + 1·max| f (N + 1)(C)|/(N + 1), (2) СÎ[a, b] где N - степень аппроксимирующего полинома (т.е. N = 0 для гистограммы и 1 для полигона частот), f (N + 1)(C) - производная (N + 1)-го порядка от f. Пусть мы хотим обеспечить некоторую наперед заданную точность RN. Как-либо изменять f (N + 1)(C) мы не можем, поскольку это - производная реально существующего (хотя пока и не известного нам) искомого распределения. Единственное, за счет чего мы можем достигнуть надлежащей точности - это величина интервала DС которую следует выбрать тем меньше, чем больше f (N + 1)(C), согласно оценке DС ~ [(N + 1)·RN /max| f (N + 1)(C)|]1/(N + 1), (2a) СÎ[a, b] Отсюда видно, что разные распределения могут потребовать разного количества интервалов группировки: чтобы наглядно представить такие распределения, у которых функция плотности меняется очень резко (т.е. максимальное абсолютное значение производной велико), понадобится большое количество относительно коротких интервалов DС, а чтобы представить распределения с плавно изменяющейся функцией плотности, можно ограничиться небольшим количеством относительно длинных интервалов (в предельном случае максимально гладкого распределения, которым является равномерное распределение, достаточно одного единственного столбца, занимающего весь интервал выборки от C1 до Ck+1 = C2). Итак, мы видим, что на самом деле количество интервалов группировки должно определяться не только имеющимся в нашем распоряжении объемом выборки, но и свойствами распределения. При этом может оказаться, что мы имеем выборку гораздо меньшего объема, чем необходимо для восстановления данного распределения (тогда несмотря на выполнение рекомендаций табл. 1 мы его восстановить не сможем). Или, напротив, может оказаться, что наша выборка гораздо больше чем необходимо для надежного восстановления некоторого относительно гладкого распределения (тогда, даже если мы сильно нарушаем рекомендации табл. 1, распределение все равно может быть успешно восстановлено). Однако формул, позволяющих рассчитать потребное количество интервалов группировки в зависимости от свойств распределения также известно несколько и они опять порождают неоднозначность. Приведем лишь два примера совершенно различных формул из [Костылев с соавт., 1991: с. 78]: k = (4/æ)·lg(n/10) (3a) и k = (Е{y} + 1.5)·n0.4/6, (3b) где æ = (Е{y})-½ - контрэксцесс распределения. Трудность использования двух последних выражений состоит в том, что число интервалов группировки нужно выбрать до того, как будут определены числовые характеристики распределения, в том числе и эксцесс или контрэксцесс. Обычно полагают, что для эмпирических распределений 1.8 ≤ Е{y} ≤ 6, (4) и получают соответствующие границы числа интервалов [Костылев с соавт., 1991: с. 79]. Именно из формул (№3b) и (№4) были получены ограничения (№1). Очевидно, что если для получения ограничений использовать формулу (№3a) вместо (№3b), то и ограничения были бы совершенно другими, а тогда более или менее стройная система табл. 1 (которая дает однозначность хотя бы для не очень больших n) нарушилась бы! Более того, ограничения эти (даже без всяких изменений - в том виде, в котором они заданы сейчас) выписаны в предположении справедливости условия (№4), а для распределений, у которых, например, Е{y} = 70, рекомендации табл. 1 (по ограничению количества интервалов группировки) будут просто неправильными. Но у способов разбиения на равные интервалы есть и более существенный недостаток. Формула (№2a) требует везде использовать малые интервалы DС, определяемые наибольшим значением производной на всем интервале выборки [C1 Ck+1]. Однако на практике мы имеем выборку ограниченного объема и в соответствие с табл. 1 не можем разбивать интервал выборки на слишком большое количество интервалов группировки. А ведь плотность распределения может на каком-то участке изменяться так быстро, что потребное количество интервалов группировки превысит имеющееся в нашем распоряжении количество измерений n. Разбиения выборки на неравные интервалы группировки Предположим, что мы хотим построить гистограмму (т.е. аппроксимировать функцию плотности распределения) не на всем интервале выборки, а лишь на его левой половине (положим для определенности, что здесь функция меняется достаточно резко). При этом по условию (№2a) будет получено некоторое значение DСL (относительно малое, поскольку резкое изменение функции количественно выражается в том, что max| f '(C)| будет достаточно большой величиной). Предположим теперь, что мы хотим построить гистограмму лишь на правой половине интервале выборки (положим для определенности, что здесь функция меняется очень плавно). При этом по условию (№2a) будет получено некоторое значение DСR (относительно большое, поскольку плавное изменение функции количественно выражается в том, что max| f '(C)| будет достаточно малой величиной). Но ничто не мешало нам рассмотреть не левую и правую части интервала выборки, а разделить его, скажем на три части. Или на четыре… Если довести эту идею до логического завершения, то мы приходим к очевидному требованию: длину каждого интервала группировки следует выбирать в соответствии с поведением функции плотности распределения - на участках быстрого изменения этой функции (т.е. на участках с большим абсолютным значением производной) следует использовать малые интервалы DС, а на участках медленного изменения функции (т.е. на участках с малым абсолютным значением производной) можно использовать большие интервалы DС. Иначе говоря, мы пришли к необходимости использования неравномерной системы интервалов. Внешне это будет выглядеть так, что гистограмма состоит из относительно узких столбиков на участках резкого убывания или возрастания функции распределения, а на участках плавного изменения ширина столбиков увеличится. Строить систему неравномерных интервалов можно различными способами и среди них есть весьма простые. Например, используют не интервалы равной длины, а интервалы с равной вероятностью. Для оценки числа таких интервалов К часто используется формула Уильямса, которая при уровне значимости 0.1 имеет вид [Костылев с соавт., 1991: с. 79]: К = 1.9·n0.4. (5) (естественно, К каким-либо образом округляется до целого числа, например, до ближайшего меньшего целого). Алгоритм построения гистограммы на неравных интервалах заключается в следующем. По вышеприведенной формуле определяется оптимальное количество таких интервалов. Выборка упорядочивается по возрастанию и определяется количество элементов (N) выборки, которые должны попасть в каждый интервал группировки: N = n/K (конечно, n может не делиться на K нацело, но мы сейчас для простоты будем считать, что получено именно целое значение N; если же деление нацело невозможно, то следует воспользоваться каким-либо простейшим чисто техническим приемом, например, в качестве N взять округленное до ближайшего целого значение n/K, и в каждый интервал группировки поместить N значений выборки, а в последний интервал - N1 оставшихся значений выборки; также можно использовать и более равномерное распределение элементов по интервалам, например, следя за тем, чтобы количество элементов выборки в соседних интервалах группировки не отличалось более чем на 1). Итак, мы имеем K интервалов группировки, в каждом из которых находится fi = N (i = 1, …, K) элементов выборки. Пусть DСi - ширина i-го интервала группировки. Значения DСi выбираются именно таким образом, чтобы в каждом i-м интервале оказалось N элементов выборки. Для этого в упорядоченной по возрастанию выборке, отсчитывается N элементов и они помещаются в первый интервал, затем отсчитывается еще N элементов и они помещаются во второй интервал и т.д. При этом в первом интервале окажутся элементы от у1 до уN включительно, во втором интервале - от уN+1 до у2·N включительно, в третьем - от у2·N+1 до у3·N и т.д. Тогда в качестве границ интервалов группировки можно принять C1 = у1, C2 = (уN + уN+1)/2, C3 = (у2·N + у2·N+1)/2, …, CК+1 = уn. (6) Отсюда DС1 = C2 - C1, DС2 = C3 - C2, …, DСi = Ci+1 - Ci. Поскольку площадь всей гистограммы должна быть равна полной вероятности, т.е. 1 (или 100%), то для определения высот столбиков гистограммы (h1, h2, …, hК) мы имеем систему уравнений: К DС1·h1·α = f1, DС2·h2·α = f2, …, DСК ·hК ·α = fК, где α = Σ fi = n. (№7) i = 1 Отсюда мы вычисляем высоты столбиков h1, h2, …, hК и строим гистограмму. К сожалению, разбиение выборки на неравные интервалы группировки также неоднозначно. В частности, кроме упомянутой выше формулы Уильямса, используется еще и формула Манна-Вальда [Костылев с соавт., 1991: с. 79]. Главный вопрос: А восстанавливается ли плотность распределения? Однако, даже если мы, используя какие-то совершенные ограничения, сможем добиться однозначности выбора интервалов, это не гарантирует нам нахождения истинной плотности распределения. Сложность состоит в том, что правильность восстановления плотности распределения никак не связана с однозначностью количества интервалов. При слишком малом числе интервалов k эмпирическая плотность вероятности (гистограмма) утратит детальность и связь с особенностями теоретического распределения станет малоинформативной. Если же число интервалов будет слишком большим, то некоторые из них окажутся пустыми или слабо заполненными, т.е. гистограмма окажется изрезанной [Костылев с соавт., 1991: с. 78], иначе говоря, эмпирическая плотность вероятности обретет множество экстремумов, которых, возможно, не было в теоретическом распределении. Но если истинная плотность распределения нам не известна, то каким же образом можно решить вопрос о том, восстанавливается ли она при том или ином разбиении интервала выборки? Сделать это можно (и весьма просто) при помощи метода математического моделирования: необходимо задаться некоторым известным распределением, сгенерировать на его основе выборку, а после этого, исходя из данной выборки, попытаться восстановить исходную плотность распределения. Уточнение вопроса: Какова погрешность плотности распределения? Анализ конкретизированного в предыдущем разделе вопроса нельзя проводить, не сделав терминологическое уточнение: а что, вообще говоря, значит «восстановить плотность распределения»? Рассмотрим пример. Пусть какая-то величина характеризуется нормальным распределением с параметрами m = 0 и σ = 1. Если мы, взяв некоторую выборку значений этой величины, установим, что распределение действительно нормальное и его параметры действительно m = 0 и σ = 1, то следует признать, что мы можем в точности восстановить плотность распределения. Между тем, совершенно ясно, что, имея выборку конечного объема, мы сможем абсолютно точно вычислить m и σ лишь случайно, а в подавляющем большинстве случаев полученные значения m и σ будут слегка отличаться от 0 и 1, да и про нормальность распределения мы сможем говорить лишь с той или иной долей уверенности, но полной гарантии у нас никогда не будет. То есть, формально, истинная плотность распределения не может быть восстановлена. Она может быть восстановлена лишь с некоторой погрешностью. Таким образом качественный вопрос «можем ли мы при помощи однозначного разбиения найти истинную плотность распределения» переходит в более глубокий количественный вопрос: с какой погрешностью мы можем найти истинную плотность распределения? Теперь становится ясно, что вскрытая выше проблема неоднозначности восстановления плотности распределения при использовании гистограмм, построенных по различным интервалам группировки, сводится именно к проблеме погрешности. Погрешность построения гистограммы весьма велика и в пределах этой погрешности равнозначны все возможные варианты гистограмм. Оценим точность метода гистограмм. Пусть DС - ширина интервала гистограммы. Если ширина интервала много меньше интервала выборки, а полное количество наблюдений n >> 1, то количество наблюдений в данном интервале имеет пуассоновское распределение и относительная статистическая погрешность числа наблюдений в данном интервале (с центром С) равна[1] δs = [n·p(С)·D]-1/2, где p(x) - плотность вероятность. Кроме статистической погрешности существует еще погрешность аппроксимации [Косарев, 2008: с. 31]. Мы уже приводили формулу для абсолютной погрешности аппроксимации - см. выше формулу (№2). Для относительной погрешности аппроксимации (δа) в частном случае кусочно-постоянной аппроксимации (т.е. для N = 0): δа = RN /p(С) = A·DСN + 1·max| р(N + 1)(C)|/(N + 1)/p(C) = A·DС·| р'(C)|/p(C), здесь мы заменили знак пропорциональности (использованный в №2) на знак равенства, но, естественно, ввели константу пропорциональности A (согласно [Косарев, 2008: с. 31], А = 0.5), а также учли, что на небольшом интервале DС можно принять производную р’ постоянной. Таким образом, суммарная относительная погрешность (δ) равна: δ = (δs2 + δа2)½ = {[n·p(С)·D]-1 + 0.25·[DС·р'(C)/p(C)]2}½. Оценка погрешности метода гистограмм была впервые получена Н.В. Смирновым в 1950 г. При оптимальном интервале группировки Dopt = [2·p(С)·n-1·р'(C)-2]1/3 достигается минимальная относительная погрешность δmin = A1·[n-1·p(С)-2·р'(C)]1/3, где A1 = 3½·4-1/3 » 1.1 [Косарев, 2008: с. 31-32]. Абсолютная погрешность δ·p = δ·p(C) = {p(С)/(n·D) + 0.25·[DС·р'(C)]2}½ (8) Вообще говоря, надо понимать, что полученная формула - это лишь грубая оценка погрешности. Поэтому неудивительно, если окажется, что p - δ·p < 0. Конечно, отрицательная плотность вероятность невозможна, поэтому в таких случаях следует исправить δ·p, приписав ей значение р. Подчеркнем также, что полученные формулы дают только точность гистограммы, т.е. кусочно-постоянной аппроксимации. Для точности полигона частот (т.е. кусочно-линейной аппроксимации) формулы будут уже другие (поскольку тогда в №2 N = 1, следовательно δа ~ DС2). В [Косарев, 2008: с. 32] указано, что для полигона частот δmin ~ n-2/5. Таким образом, мы видим, что метод гистограмм не самый точный - убывание погрешности с ростом n происходит медленнее, чем для метода полигона частот [Косарев, 2008: с. 32]. В заключение этого раздела укажем, что существует и иной подход к вычислению погрешности области определения - построение доверительного интервала плотности распределения. Задача доверительного оценивания восстановленной плотности распределения pn(C) как случайной функции состоит в построении вокруг нее некоторой “трубки”, размер которой определяется вероятностными характеристиками нормированного отклонения pn(С) от истинной плотности распределения p(С), объемом исходной статической выборки и заданным коэффициентом доверия β [Лапко с соавт., 1996: с. 52]. Для простейшего случая гистограммы с равными интервалами группировки доверительная область может быть построена при помощи теоремы Мании, как это изложено, например, в [Лапко с соавт., 1996: с. 52-53]. Особенность задачи восстановления распределений в почвоведении Прежде чем поставить описанный выше вычислительный эксперимент по восстановлению некоторой известной нам плотности распределения, вычленим особенность задачи восстановления, характерную для почвоведения и экологии. Обычно, когда в литературе рассматривается вопрос о построении эмпирической плотности распределения (в виде гистограммы или полигона частот) для некоторой выборки у1, …, уn, то предполагается, что каждый ее элемент точно известен и может быть поэтому отнесен в какой-то конкретный интервал группировки (см., например, [Ллойд и Ледерман, 1989: с. 96-98; Белов с соавт., 1990: с. 221-225; Дмитриев, 1995: с. 54, 55, 172, 173]). Но, с другой стороны, хорошо известно, что данные эксперимента всегда зашумлены некоторой погрешностью, причем на это указывают в том числе и те авторы, (см. [Ллойд и Ледерман, 1989: с. 10; Белов с соавт., 1990: с. 18-23; Дмитриев, 1995: с. 7]), которые при построении гистограмм рассматривают лишь данные без погрешностей. Как видим, эти две взаимоисключающие точки зрения встречаются в одних и тех же литературных источниках. При общем введении в статистику постулируется, что экспериментальные данные всегда заданы с некоторой погрешностью, но при построении гистограммы плотности распределения этих экспериментальных данных об их погрешностях уже не вспоминают. А может быть такой подход имеет право на существование? Действительно, зачем лишний раз вспоминать о погрешности данных, если мы собираемся восстановить полную характеристику случайной величины - ее закон распределения. Ведь «погрешность» - это одна из очень упрощенных характеристик, которыми можно пытаться приближенно охарактеризовать распределение… С этим можно согласиться в том случае, когда распределение измеряемой случайной величины обусловлено только погрешностями измерений. Но при работе с такой сложной средой как почва гораздо чаще встречается иная ситуация. Поясним ее на примере. Рассмотрим эмиссию (из почвы в атмосферу) метана - важного парникового газа. Во многих работах, посвященных математическому моделированию этого процесса (см., например, [Cao et al., 1995; Walter et al., 1996; Glagolev, 1998; Глаголев, 2006]) не без оснований принимается, что интенсивность эмиссии пропорциональна произведению нескольких множителей, каждый из которых отражает влияние одного независимого фактора (логические основания этого вскрыты в [Глаголев с соавт., 2007: с. 203-204], а экспериментальные доказательства независимости действия гидротермических факторов на процесс разложения органики - одну из стадий метаногенеза, приведены в [Van der Linden et al., 1987]). Случайные величины, формирующиеся под воздействием большого числа независимых случайных факторов, взаимодействующих мультипликативно и имеющих один порядок значений, описываются логнормальным распределением [Костылев с соавт., 1991: с. 62]. Действительно, Н.С. Паников [1995: c. 20] опубликовал результаты измерения сотрудниками его лаборатории эмиссии метана на территории Бакчарского болота (Северо-Восточный отрог Большого Васюганского Болота, Бакчарский район Томской обл.), которые удовлетворяли именно логнормальному распределению (подчеркнем, что мы сейчас рассматриваем всего лишь один конкретный пример и вовсе не хотим сказать, что эмиссия удовлетворяет логнормальному распределению всегда, более того, в [Шнырев и Глаголев, 2007] приведены примеры иных распределений эмиссии). Следует отметить, что пока мы не вводили в рассмотрение какие-либо погрешности. Логнормальное распределение потока метана возникает не из-за неточности измерений, а из физической картины влияния факторов среды на образование, транспорт и окисления метана в почве. Таким образом, даже если бы мы могли произвести все измерения абсолютно точно, все равно мы бы не получали всегда какое-то одно значение величины эмиссии, но выборку различных значений у1, …, уn, удовлетворяющих, например, логнормальному распределению. Однако не следует забывать, что мы не можем произвести анализ газа абсолютно точно. Поэтому каждый элемент выборки (уi) будет получен с некоторой абсолютной ошибкой (Dуi). В почвоведении распределение ошибок анализов обычно неплохо аппроксимируется нормальным законом [Дмитриев, 1995: с. 73]. Вообще говоря, предположение о том, что закон распределения ошибок измерения близок к нормальному, принимается достаточно часто и в других предметных областях (см., например, [Белов с соавт., 1990: с. 20; Живописцев и Иванов, 1993: с. 8-9]). Таким образом, для почвоведения и экологии будет характерна задача построения эмпирической плотности распределения некоторой выборки у1, …, уn, при том, что элементы выборки заданы не точно, а имеют погрешности (соответственно, Dу1, …, Dуn), распределенные по некоторому закону, чаще всего - по закону Гаусса. Вычислительные эксперименты и их обсуждение Восстановление распределения потока СН4 Итак, допустим, что для какого-то объекта (например, для конкретного болота) на самом деле эмиссия метана имеет логнормальное распределение с плотностью вероятности æ (2·π·σ2·y 2)-½·exp{-0.5·[ln(y) - ln(m)]2/σ2} при y ≥ 0 W(y) = í è 0 при y < 0 с параметрами m = 5 и σ = 1.9. Задача определения статистических характеристик распределения потока метана для различных объектов в разных природных зонах Западной Сибири впервые решалась в [Шнырев и Глаголев, 2007]. На основании этой работы можно принять, что средний объем выборки для одного объекта составляет порядка 50 значений. Создадим выборку такого объема при помощи генератора случайных чисел (мы использовали функцию lognrnd из «Statistics Toolbox» системы MATLAB 7; результат работы данного генератора приведен в табл. 2, где мы уже упорядочили полученную выборку по возрастанию элементов). Таблица 2 Выборка, полученная при помощи генератора логнормального распределения: lognrnd(1.6,1.9,50,1) i yi i yi i yi i yi i yi 1 0.08 11 1.03 21 4.47 31 8.75 41 25.52 2 0.21 12 1.08 22 4.65 32 9.31 42 37.95 3 0.24 13 1.34 23 4.82 33 13.65 43 47.89 4 0.24 14 1.63 24 5.60 34 14.80 44 48.04 5 0.32 15 2.20 25 6.21 35 17.81 45 48.05 6 0.39 16 2.34 26 6.34 36 18.55 46 54.16 7 0.51 17 3.51 27 6.97 37 19.34 47 58.03 8 0.57 18 3.71 28 7.58 38 19.43 48 73.56 9 0.67 19 3.86 29 8.15 39 19.85 49 109.31 10 0.87 20 4.17 30 8.64 40 23.55 50 316.54 Теперь на каждое измерение нам нужно наложить «погрешности измерений», распределенные в соответствие с нормальным законом - именно результирующая выборка будет представлять собой модель реальных измерений потока метана в природе. Но пока мы этого делать не будем, а зададимся вопросом - насколько точно смогут восстановить исходную плотность распределения (если даже элементы выборки не «зашумлены» никакими дополнительными погрешностями) стандартные методы построения гистограмм. Правда, для большей наглядности (чтобы несоответствие между теоретическим исходным и восстановленным распределениями не «скрадывалось» за счет ширины столбиков гистограммы) мы будем использовать графическое представление не в виде гистограммы, а в виде сглаженного полигона частот, однако подчеркнем, что это - только способ графического представления данных, количество же интервалов группировки, их расположение и ширина совершенно одинаковы что при построении гистограммы (кусочно-постоянной аппроксимации плотности распределения), что при построении обычного полигона частот (кусочно-линейной аппроксимации), что при использовании сглаженного полигона частот (аппроксимации кусочными полиномами порядка выше первого, в частности - при широко используемой в настоящее время аппроксимации сплайнами). Прежде всего сосредоточимся на формуле Старджеса, как наиболее широко применяемой в практической деятельности. Сама по себе эта формула для n = 50 дает k = 1 + log2(50) » 6.6 интервала. Поскольку нецелое количество интервалов невозможно, то рассмотрим несколько вариантов: округление до ближайшего целого, в данном случае совпадающее с округлением до ближайшего большего целого (kR = 7), округление до ближайшего меньшего целого (kL = 6). Оба эти значения удовлетворяют ограничениям Живописцева-Иванова (№1а). Но они не удовлетворяют верхнему ограничению Костылева с соавт. (№1), действительно 1.25·500.4 » 5.98. Следовательно, нам придется проверить еще kLО = 5 интервалов - это количество интервалов группировки удовлетворяет всем ограничениям, в том числе и «рекомендации нечетности». Также обратим внимание на формулу Живописцева-Иванова с обязательным ограничением (№1а), поскольку среди всех формул она дает максимальное количество интервалов - 30. Из табл. 1 видно, что y1 = 0.08, y50 = 316.54. Следовательно при kR = 7 ширина интервала группировки будет составлять DС = 45.2086. При kL = 6 DС = 52.7433, а при kLО = 5 DС = 63.2920. Однако ни одно из этих значений не является удовлетворительным - из рис. 2 видно, что эмпирические плотности распределений, построенные с использованием указанных интервалов группировки не имеют никакого отношения к исходному теоретическому распределению. На первый взгляд удивительно то, что выделение интервалов группировки по формуле Живописцева-Иванова дает наилучший результат. Однако на самом деле ничего удивительного здесь нет. Мы специально выбрали очень «неудобное» распределение с огромным эксцессом (около двух миллионов!!!). Для распределения с таким эксцессом ограничение (№1), выведенное в предположении того, что эксцесс не превышает 6 (!!!), безусловно, не действует. Например, формула (№3a), учитывающая величину эксцесса, рекомендует разбить интервал выборки почти на 4000 равномерных интервалов группировки! Конечно, по сравнению с этим огромным числом кажутся одинаково неправильными разбиения что на 5-7 (по формуле Старджеса), что на 30 интервалов (по формуле Живописцева-Иванова), но, все-таки, 30 интервалов ближе к 4000, чем 5 интервалов! Глаголев%20и%20Сабреков,%20Рис Рис. 2. Восстановление плотности распределения при помощи интервалов группировки равной длины. Неужели действительно необходимо разбивать выборку на такое огромное количество интервалов группировки? Да это и не будет иметь какого-либо смысла для объема выборки n = 50! Сейчас мы покажем, что такой элементарный прием, как использование интервалов неравной длины (но равной вероятности) намного улучшает ситуацию даже при необычайно высоком эксцессе. По формуле (№5) находим, что в случае интервалов равной вероятности их необходимо лишь К = 1.9·500.4 » 9.1. Поскольку дробное количество интервалов невозможно, то опять рассмотрим два варианта: округление до ближайшего целого, в данном случае совпадающее с округлением до ближайшего меньшего целого (КL = 9) и округление до ближайшего большего целого (КR = 10). Последний вариант очень прост - интервалы должны быть таковы, чтобы в каждый из них попало по 5 подряд идущих выборочных значений. Из табл. 2 по формуле (№6) получаем: [0.08, 0.355); [0.355, 0.975); [0.975, 2.27); [2.27, 4.32); [4.32, 6.275); [6.275, 8.695); [8.695, 18.18); [18.18, 24.535); [24.535, 51.105); [51.105, 316.54] Тогда DС1 = 0.275, DС2 = 0.620, DС3 = 1.295, DС4 = 2.050, DС5 = 1.955, DС6 = 2.420, DС7 = 9.485, DС8 = 6.355, DС9 = 26.570, DС10 = 265.435. Учитывая, что для нашего примера fi/α = 5/50 = 0.1, получаем по формуле (7): h1 = f1/(α·DС1) = 0.1/0.275 » 0.364, h2 » 0.161, h3 » 0.077, h4 » 0.049, h5 » 0.051, h6 » 0.041, h7 » 0.011, h8 » 0.016, h9 » 0.004, h10 » 0.000. Для КL = 9 уже не удастся выбрать такие интервалы, чтобы в каждый из них попало равное количество элементов выборки. Как же быть в этом случае? Тут нет однозначного решения. Понятно, что разница между количествами элементов, попавших в соседние интервалы должна быть минимальна, т.е. эти количества могут различаться на 1. Если в 4 интервала попадет по 5 элементов в каждый, а в 5 интервалов попадет по 6 элементов в каждый, то нам как раз удастся разместить все 50 элементов выборки по 9 интервалам группировки. Но как должны быть расположены «пятиричные» и «шестиричные» интервалы один относительно другого? Чередоваться? Или сначала подряд должны идти интервалы одного вида, а потом другого? Здесь опять появляется источник неоднозначности. Впрочем, на практике оказывается, что эта неоднозначность не порождает больших расхождений в результате - независимо от расположения интервалов восстановленная плотность вероятности имеет примерно одинаковый вид. Представляется разумным интервалы с большим количеством элементов размещать на длинных «хвостах» распределений. Раз элементов больше, то, следовательно, больше степень сглаживания (усреднения). Поскольку обычно на «хвостах» значения плотности вероятности достаточно малы, на фоне этих малых значений проявляются «паразитные» высокочастотные колебания, которых в исходном теоретическом распределении не было. Большая степень сглаживания позволяет избавиться от таких осцилляций. Иногда поступают совершенно формально. Если, например, получается, что в каждый интервал группировки должно попасть по 5.5 элемента (что, казалось бы, лишено смысла), то отсчитывают 5 элементов, помещают их в один класс, а 6-ой элемент помещают как в этот класс, так и в следующий. При этом граница между классами проходит именно по этому 6-му элементу. На рис. 3 приведены плотности распределения, восстановленные с использованием 9 и 10 неравномерных интервалов. Прежде всего, следует отметить, что теоретическое распределение восстановлено, вообще говоря, неплохо (обращаем внимание читателя, что по оси абсцисс выбран логарифмический масштаб, поэтому кажущееся несоответствие между исходной теоретической и восстановленной по «экспериментальным данным» плотностями распределения в левой части рисунка на самом деле имеет место для небольшого интервала очень малых эмиссий метана. Надо понимать, что для ограниченной выборки достигнуть полного соответствия теоретическому распределению принципиально невозможно. В той конкретной выборке, которую мы имеем, минимальное значение оказалось равным 0.08 (см. табл. 2). Окажись оно, скажем, 0.04, соответствие было бы лучше. Во-вторых, отметим, что, как и следовало из общетеоретических соображений, при использовании меньшего количества интервалов удалось избежать осцилляций на правом «хвосте» функции. К сожалению, масштаб рис. №3 не позволяет это увидеть (кривые, соответствующие 9 и 10 интервалам, практически сливаются), поэтому мы использовали врезку иного масштаба. Конечно, непосредственно из рис. №3 может показаться, что это не очень существенно. Однако не будем забывать о логарифмическом масштабе! В-третьих, мы видим, что абсолютная погрешность, рассчитанная по формуле №8, действительно очень реалистично отражает возможное отклонение восстановленной нами плотности распределения от истинной. Естественно, рассчитывая абсолютную погрешность, мы не использовали какой-либо информации об истинном распределении, а оценивали значения p(С) и p'(С) по восстановленной плотности. Глаголев%20и%20Сабреков,%20Рис Рис. 3. Восстановление плотности распределения при помощи интервалов группировки равной вероятности (но неравной длины). Завершая данный раздел, заметим, что использование интервалов неравной длины позволило достаточно точно восстановить по «экспериментальным данным» плотность распределения, имеющего чрезвычайно высокий эксцесс. Причем сделано это было, так сказать, совершенно «автоматически» - нам не понадобилась какая-либо информация о свойствах восстанавливаемого распределения, в частности, о величине эксцесса. Обычные методы, основанные на использовании интервалов равной длины оказались совершенно неприменимы - построенные с их помощью распределения чрезвычайно сильно отличались от исходного логнормального распределения. Поскольку априори свойства восстанавливаемого распределения нам не известны, то пользоваться равными интервалами группировки никогда не следует. Некорректность задачи восстановления плотности распределения С математической точки зрения постановка любой задачи включает в себя задание множества допустимых входных данных Ф и множества возможных решений П. Цель вычислительной задачи состоит в нахождении решения p Î П по заданным входным данным F Î Ф. Вычислительная задача называется корректной (по Адамару-Петровскому), если выполнены следующие требования (в том случае когда хотя бы одно из нижеперечисленных требований не выполнено, задача называется некорректной) [Амосов с соавт., 1994: с. 44]: 1) решение задачи существует при любых входных данных F Î П; 2) это решение единственно; 3) решение устойчиво по отношению к малым возмущениям входных данных (отсутствие устойчивости означает, что малым погрешностям входных величин F могут соответствовать сколь угодно большие возмущения идентифицируемого элемента p [Мацевитый и Лушпенко, 1990]). Согласно известному определению, р(С) = dF(C)/dC, где F(z) - функция распределения [Лапко с соавт., 2000: с. 32]. Если вычислять плотность распределения согласно этому определению, то необходимо будет дифференцировать эмпирическую функцию распределения. Существует и другой подход к задаче. Интегрируя вышеприведенное определение можно получить интегральное уравнение, которому должны удовлетворять функции распределения и плотности распределения. Таким образом, функция р(С) является решением интегрального уравнения Фредгольма 1-го рода òθ(z-C)·р(С)dC = F(z), где θ(z-C) - функция Хевисайда [Карандеев и Эйсымонт, 1998]. Некорректность задач численного дифференцирования и решения интегральных уравнений 1-го рода хорошо известна и подробно разъяснялась в литературе (см., например, [Калиткин, 1978: с. 82, 462-464; Тихонов и Арсенин, 1979: с. 9, 11-13, 18]). Для решения некорректных задач применяют специальные методы, использующие так называемую регуляризацию решения. В следующем разделе мы перечислим целый ряд методов восстановления плотности распределения, каждый из которых основан на том или ином способе регуляризации. Какие еще существуют методы восстановления плотности распределения? Нельзя сказать, что задача восстановления плотности распределения относится к числу совсем неизученных, однако по каким-то причинам почти все, достигнутое в этой области, осталось вне поля зрения большинства почвоведов и экологов, сталкивающихся с этой задачей в своей профессиональной деятельности. Не претендуя на полное освещение методов решения указанной задачи, мы, тем не менее, перечислим ряд из них и приведем ссылки, которые, как мы надеемся, позволят заинтересованному читателю подробно познакомиться с каждым методом (рассмотренный выше метод гистограмм мы вторично упоминать не будем): 1. Дескриптивное приближение сплайнами [Воскобойников с соавт., 1984]; 2. Интегральная оценка [Лапко с соавт., 1996: с. 25-29]. 3. Метод Карандеева-Эйсымонта [Карандеев и Эйсымонт, 1998]; 4. Метод корневой оценки [Крянцев и Лукин, 2006: c. 71-73]; 5. Метод Розенблатта-Парзена [Лапко с соавт., 1996: с. 23-25; Крянцев и Лукин, 2006: c. 63-65; Лагутин, 2007: с. 388-392], называемый также методом Е. Парзена (1962 г.) или ядерным методом [Косарев, 2008: с. 32-33]; 6. Метод стохастической регуляризации [Лапко с соавт., 1996: c. 34-36]; 7. Метод структурной минимизации риска [Вапник с соавт., 1984: с. 688-706]; 8. Оценки максимума правдоподобия (в частности, метод решета) [Лапко с соавт., 1996: c. 31-32]; 9. Проекционные методы [Крянцев и Лукин, 2006: c. 66-69], в частности метод Н.Н. Ченцова (1962 г.) [Косарев, 2008: с. 32-33]; 10. Регуляризованный метод гистограмм [Крянцев и Лукин, 2006: c. 69-71]; 11. Рекуррентные ядерные оценки [Лапко с соавт., 1996: c. 30-31]; 12. Ядерная оценка с переменным параметром сглаживания [Лапко с соавт., 1996: с. 29-30]; 13. Ядерные оценки с пониженным смещением [Лапко с соавт., 1996: c. 32-34]. В [Лапко с соавт., 1996: с. 29] проводилось сравнение методов Розенблатта-Парзена и интегральной оценки при решении двух тестовых задач восстановления плотности на основе выборок объемом n = 50, 100, 200: (а) нормального распределения, (б) распределения рЭ(х)=0.5/[π·(x2+1)]. Оказалось, что более точное приближение дает интегральная оценка. В [Карандеев и Эйсымонт, 1998] проводилось сравнение методов Розенблатта-Парзена (с экспоненциальным ядром) и Карандеева-Эйсымонта при решении тестовых задач восстановления плотности: (а) смеси трех нормальных распределений на основе выборки объемом n = 40, (б) распределения Коши по выборке объемом n = 20 и (в) гамма-распределения. Оказалось, что для выборок малого объема более точное приближение дает метод Карандеева-Эйсымонта. Кроме того, для этого метода оценка является более устойчивой относительно выбора константы регуляризации. При оптимальном выборе ширины σ гауссова ядра погрешность метода Парзена в зависимости от числа точек δmin ~ n-2/5, т.е. такая же как и в методе полигона частот [Косарев, 2008: с. 32]. Описание «метода структурной минимизации риска» в [Вапник с соавт., 1984] снабжены текстом программы DENSIT на языке FORTRAN и отладочным примером с подробной распечаткой результатов. В 1978 г. было показано, что не существует никаких других методов оценки плотности, имеющих большую скорость убывания (по n) погрешности оценки плотности, чем метод Ченцова (в этом методе δmin ~ M/n½, где М - число используемых базисных функций) [Косарев, 2008: с. 33]. Преимущество рекуррентных процедур оценивания плотности вероятности заключается в относительной простоте настройки их параметров при получении новых данных. [Лапко с соавт., 1996: c. 30]. Исследование ядерных оценок с переменным параметром сглаживания, представляющих собой сложную сумму зависимых случайных величин, к сожалению, находится пока еще на начальном этапе [Лапко с соавт., 1996: с. 29-30]. В заключение это раздела заметим, что вышеприведенный список (непараметрических методов восстановления плотности распределения) достаточно условен - он не должен рассматриваться как строгая классификация. Например, наиболее популярная непараметрическая оценка - метод Розенблатта-Парзена, которую мы, следуя традиции, поддерживаемой практически всеми исследователями (см., например, [Лапко с соавт., 1996: с. 23-25; Крянцев и Лукин, 2006: c. 63-65; Лагутин, 2007: с. 388-392]), выделяем отдельным пунктом, может рассматриваться как частный случай других методов. Например, в [Лапко с соавт., 1996: с. 26] показано, что при некоторых конкретных значениях параметров интегральная оценка становится тождественной оценке Парзена. Совершенно очевидно также, что оценка Розенблатта-Парзена (получаемая, как известно, при постоянном значении параметра сглаживания) - это частный случай более общей ядерной оценки с переменным параметром сглаживания. Абсолютно то же самое можно сказать и про взаимоотношение оценки Розенблатта-Парзена с рекуррентными ядерными оценками. Наконец, в [Лапко с соавт., 1996: с. 35] показано, как получить одну конкретную оценку Розенблатта-Парзена методом стохастической регуляризации. Вообще для многих конкретных реализации того или иного метода можно обнаружить их принадлежность одновременно к одной и другой группе непараметрических методов из вышеприведенного списка. Так, сравнивая статистику Бримэна, принадлежащую к классу ядерных оценок с переменным параметром сглаживания, и оценку Валветэна-Вагнера, относящуюся к классу рекуррентных ядерных оценок (см., например, [Лапко с соавт., 1996: с. 29-30]), видим их формальную тождественность. Некоторое программное обеспечение Из перечисленных выше методов непараметрического восстановления плотности распределения на сегодняшний день, пожалуй, наиболее популярным является ядерный метод (Розенблатта-Парзена). В частности, именно этот метод реализован в широко распространенной системе MATLAB, где ему соответствует функция ksdensity. Но хороши ли результаты, получаемые с помощью данного метода? Прежде всего заметим, что источниками неоднозначности при использовании функции ksdensity являются: 1) тип ядра (иными словами - вид «окна сглаживания»); 2) ширина «окна сглаживания»; 3) границы отрезка, на котором восстанавливается плотность распределения. Однако при практическом оценивании плотности по заданной реализации выборки важен не столько вид окна, сколько правильное определение ширины «окна сглаживания» [Лагутин, 2007: с. 392]. Если границы отрезка нельзя указать из каких-то теоретических предпосылок, то можно использовать некоторое эмпирическое правило, например, формулу Стефанюка. Эта формула, в частности, используется в программе DENSIT (реализующей «метод структурной минимизации риска»): v = 5·(ymax - ymin)/(n -1), А = ymin - v, B = ymax + v, где ymax, ymin - соответственно, максимальное и минимальное значения элементов выборки; А, B - соответственно, координаты левой и правой границ интервала восстановления плотности вероятности [Вапник с соавт., 1984: с. 692]. В [Лагутин, 2007: с. 390] показано, что главная часть квадратичного риска оценки плотности[2] определяется формулой J(С) = 0.25·β2· р''(C)·hn4 + α·р(C)/(n·hn), (9) где hn - ширина «окна сглаживания»; α, β могут быть выражены через ядро q: α = òq2(y)dy, β = òy2·q(y)dy. Глаголев%20и%20Сабреков,%20Рис Рис. 4. Восстановление плотности методом Розенблатта-Парзена. Для ядра Гаусса q(y) = (2·π)-½·ехр(-у2/2) имеем [Бронштейн и Семендяев, 1980: с. 122]: +¥ +¥ +¥ α = òq2(y)dy = (2·π)-1·òехр(-у2)dy = 2·(2·π)-1·òехр(-у2)dy = π-1·π½/2 » 0.2821. -¥ -¥ 0 +¥ +¥ +¥ β = òy2·q(y)dy = 2·π-½·ò(y/2½)2·ехр[-(у/2½)2]d(y/2½) = 4·π-½·òx2·ехр(-x2)dx = 1. -¥ -¥ 0 Следовательно, можно так подобрать ширину «окна сглаживания» (hn), чтобы минимизировать òJ(C)½dС, т.е. минимизировать площадь полосы доверительного интервала вокруг восстанавливаемой плотности распределения. Естественно, при этом, изменяя hn, необходимо следить (по какому-либо статистическому критерию, например, по критерию Колмогорова-Смирнова), чтобы получающаяся плотность распределения соответствовала исходной выборке. Впрочем, функция ksdensity может автоматически определять величину hn. Но восстанавливается ли при этом истинная плотность распределения? Вновь проведем вычислительный эксперимент. На рис. №4 приведены результаты восстановления плотности распределения по выборке табл. 2 MATLAB-функцией ksdensity. Сделаем несколько замечаний к этому рисунку. 1. Поскольку границы отрезка, на котором должна быть восстановлена плотность распределения, нельзя указать из каких-то теоретических предпосылок, то использовалась формула Стефанюка. Ширина «окна сглаживания» выбиралась функцией ksdensity в автоматическом режиме. Веса элементов выборки задавались обратно пропорциональными квадратам их абсолютных значений (если измерения не являются равноточными, то в статистике обычно вес индивидуального измерения принимается обратно пропорциональным дисперсии этого измерения, т.е. обратно пропорциональным квадрату средней квадратической ошибки данного измерения [Румшиский, 1971: с. 26-27]; в реальности оказывается, что потоки метана из почвы удается измерить с некоторой более или менее постоянной относительной ошибкой; таким образом, абсолютная ошибка или пропорциональная ей средняя квадратическая ошибка оказываются пропорциональны абсолютному значению величины потока, следовательно, дисперсия, равная квадрату средней квадратической ошибки, будет пропорциональна квадрату величины потока, а вес, обратно пропорциональный дисперсии, будет обратно пропорционален квадрату величины потока). 2. На первый взгляд кажется, что плотность распределения восстановлена не очень хорошо, однако обратим внимание на логарифмический масштаб по оси абсцисс. «Не очень хорошо» она восстановилась лишь в небольшом диапазоне значений потока метана от 0.01 до 1. 3. Отложенная на рис. №4 (в виде вертикальных «усов») оценка доверительного интервала при величинах потока порядка сотых-десятых не «дотягивается» до истинной функции плотности распределения. Т.е. мы, вроде бы, недооцениваем погрешность[3] восстановления плотности распределения. Но следует понимать, что J(С)½ - это аналог лишь средней квадратической ошибки, и как индивидуальные измерения вполне могут отклоняться от среднего значения более чем на величину средней квадратической ошибки, так и восстановленная функция местами может отклоняться от истинных значений плотности распределения более чем на J(С)½. Наконец, поставим еще один вычислительный эксперимент, более всего приближенный к реальности. Пока мы пытались восстановить плотность распределения по выборке, каждый элемент которой был точным[4] результатом работы генератора случайных чисел. Но при реальных измерениях точные результаты не могут быть нам известны. Индивидуальные измерения потока метана, имеющего в природе некоторое распределение (пусть, например, логарифмически-нормальное), получаются нами с некоторыми погрешностями. Можем ли мы восстановить исходное распределение, если измерения зашумлены погрешностями? Чтобы смоделировать эту ситуацию, «испортим» каждое значение выборки табл. 2 погрешностью, распределенной по нормальному закону. Для реальных измерений потока погрешность в среднем составляет десятки процентов (см., например, [Panikov et al., 1997; Kotsyurbenko et al., 2004; Глаголев с соавт., 2007; Шнырев и Глаголев, 2007]). На рис. №5 приведены результаты восстановления плотности распределения функцией ksdensity по некоторым использованным в вычислительном эксперименте «испорченным» шумом выборкам. На этом рис. мы не приводим доверительную область, поскольку она практически такая же, что и на рис. №4. Из рис. №5 очевидно, что качество восстановления по зашумленным данным практически такое же, что и для «точных» данных. На первый взгляд совершенно удивительной, почти мистической, кажется возможность восстановления плотности распределения по данным, зашумленным на 160%. Достаточно сказать, что, например, 50-й элемент выборки превратился из «истинного» значения 316.54 в 50.82! Как же можно по настолько сильно «испорченной» выборке восстановить исходную плотность распределения?! Конечно, по значению 50.82 уже никак не восстановишь 316.54! Но не будем забывать о погрешностях и весах!!! Значение 50.82 имеет очень большую погрешность (265.72), а, значит, при построении плотности распределения это значение будет учитываться с пренебрежимо малым весом (1/265.722 » 1.4·10-5), т.е., фактически, 50-й элемент выборки, вообще не будет учитываться при восстановлении плотности распределения. Точно так же автоматически не будут учитываться и другие «совсем испорченные» элементы. Можно сказать, что чем больше шум, тем меньше элементов выборки используется для восстановления плотности распределения. Глаголев%20и%20Сабреков,%20Рис Рис. 5. Восстановление плотности распределения по зашумленным данным методом Розенблатта-Парзена. [1] В [Косарев, 2008: с. 31] принимается, что D = DС. Однако в общем виде это не может быть так, поскольку тогда относительная ошибка (δs) приобретает размерность, обратно пропорциональную квадратному корню из размерности С. Более того, получается, что сама величина относительной ошибки плотности вероятности зависит от того, в каких единицах мы измеряем DС. А ведь это - полный абсурд! Например, если DС измеряется в единицах длины, то достаточно перейти, от метров к сантиметрам, чтобы δs уменьшить в 10 раз!!! Видимо, следует принять D = DС/(Ck+1 - C1) [2] являющаяся аналогом квадрата абсолютной погрешности (№8), аналогом же самой абсолютной погрешности будет, естественно, J(С)½. [3] Кстати, заметим, что MATLAB-функция ksdensity вообще не дает какой-либо оценки погрешности; «усы» рассчитаны с использованием формулы (№9) и представляют собой ± J(С)½. [4] Точнее говоря, «относительно точным». Из табл. 2 видно, что мы везде осуществили округление до 2-го десятичного знака, следовательно погрешность задания наименьшего элемента выборки составляет около 6% (0.005/0.08·100% = 6.25%). Для наибольшего элемента выборки эта погрешность совершенно незначительна: 0.005/316.54·100% = 0.0016%.