CC BY
46
ИССЛЕДОВАНИЕ ГИСТОГРАММ ГЕОЛОГИЧЕСКИХ ПРИЗНАКОВ КОМПЬЮТЕРНЫМ МОДЕЛИРОВАНИЕМ
Д. г.-м. н.
Ю. А. Ткачев
tkachev@geo. котг^с. ги
Гистограмма, как известно, представляет собой столбчатый график распределения частот по интервалам значений признака. Анализ гистограмм традиционно считается среди геологов наглядным и информативным методом решения геологических задач. Даже простое сведение цифровой информации в рисунок приносит большую пользу. Гистограмма позволяет единыш взглядом охватить всю выборку и составить представление о распределении изучаемой величины, чего нельзя сделать даже внимательным изучением таблиц исходных значений. Анализ гистограмм позволяет проверять геологические гипотезы, сформулированные на языке статистики. Одним из наиболее простых действий приемов анализа гистограмм является проверка гипотезы о том, что распределение данных подчиняется предполагаемому, например нормальному (гауссовскому), закону. Техника проверки заключается в сравнении эмпирических (пэ) и теоретических (пт) частот. Расчет теоретических частот производится при условии, что выборка данных получена из генеральной совокупности с предполагаемым, например гауссовским, распределением. Сравнить эти частоты с эмпирическими можно с помощью критерия “согласия” хи-квадрат:
2 (п - п )
х2 =Е^—э—
В зависимости от полученного значения х2 принимают решение. Если значение х2 превышает критическое, х 2,«,
где / — число степеней свободы, /= к - V, где к—число интервалов гистограммы, V — число параметров “подозреваемого” закона, оцененных по эмпирическим данным, а — уровень значимости (так называемая ошибка первого рода — вероятность отклонения верной гипотезы), то гипотезу отклоняют, если значение х2 меньше критического — не отклоняют. Г ипотезу о принадлежности выборки данному за-
кону называют нулевой Н0. В данном случае
Но : /(х) = О, где О — гауссовское распределение. Рассмотренная нулевая гипотеза является простой (независимо от того, сложен или несложен предполагаемый закон), потому что теоретические частоты по этому закону являются вполне определенными, вычислимыми.
Как это ни парадоксально, проверить “менее сильные” гипотезы оказывается значительно сложнее. Под менее сильной гипотезой мы подразумеваем гипотезу более общего, менее притязательного характера, как правило, менее специальную. Например, менее сильной гипотезой является гипотеза, согласно которой распределение симметрично относительно того или иного интервала гистограммы. Менее сильные гипотезы мы называем сложными, так как для них могут существовать много вариантов теоретических частот.
Если верна более сильная гипотеза, то подавно верна и менее сильная, но не наоборот. Например, из гауссовского закона следует, что распределение симметрично, из симметричности гистограммы не следует, что распределение подчиняется гауссовскому закону. Образование планеты Земля из плане-тезималей — гипотеза менее сильная, чем гипотеза образования её из твердых частиц размером до 1 мм.
Менее сильной гипотезой, чем принадлежность распределения гауссовскому закону, является гипотеза, что распределение, представленное на гистограмме, одномодально. Имеющимися средствами математической статистики в общем случае такой вопрос разрешить сложнее, чем установить, подчиняется ли распределение указанному наперед закону. Если бы у нас была только одна определенная альтернатива многомодальному распределению (например, нормальное распределение), задача решалась бы просто: необходимо было бы проверить гипоте-
зу о нормальности и в случае её отклонения считать, что наше распределение би- или полимодальное. Однако существует бесчисленное множество теоретических вариантов одномодальных распределений, не подчиняющихся гауссовскому закону. Наш вариант проверки может незаконно их отклонить.
В наиболее общем виде задача, поставленная в статье, формулируется следующим образом: на заданном уровне значимости проверить гипотезу о том, что представленное на гистограмме распределение является моно-модальным. Поскольку мономодальных распределений бесконечное количество, проверяемая гипотеза оказывается сложной:
Н0 : {(^э = ^) V (^ = ^) V ...
...(^ = Ц)... V (^ = ¥п)},
п ^ го (1)
Точнее, следовало бы в записи (1) писать не ¥э = ¥, а “эмпирическое распределение ¥э есть распределение ¥1 или ^2, или ... ^п”. Еще точнее нулевую гипотезу необходимо было бы сформулировать так: “выборка, представленная гистограммой, является выборкой из одномодального распределения ^ или ^2, или ... Еп”. Именно в этом смысле мы будем понимать запись (1).
Проверить бесконечное число гипотез, образующих в совокупности сложную слабую гипотезу, невозможно, и мы вынуждены сделать ряд упрощений.
1. Распределения: а) проверяемое эмпирическое и б) теоретические одномодальные, будем считать дискретными. Конкретнее, случайная величина в них может принимать к значений — по числу интервалов проверяемой гистограммы.
2. Частоты значений в каждом теоретическом распределении (т.е. фактически сами распределения) определяются числом элементов, размещенных в / -том интервале. Общее число элементов при таком расчете равно объему
и
т
выборки п, т.е. = —, причем п эле-
п
ментов должны быть размещены в к интервалах-ячейках так, чтобы образовалась мономодальная гистограмма. Каждое такое размещение создает одно теоретическое мономодальное распределение.
Общее число размещений п элементов по к ячейкам (в каждой ячейке может быть п/ элементов при условии п > п1 > 0 и X п, = п) рассмотрим на следующих примерах.
Вариант 1. Все элементы сосредоточены в одном интервале, остальные интервалы пусты: М-^ = к, где М-^ — число простых гипотез. Вариант возможен при любом соотношении п и к. Гипотезы из варианта 1 проверены быть не могут, так как число степеней свободы равно нулю: к = 1,/ = к - 1 (одна связь наложена: общая численность теоретических частот равна таковой в эмпирической гистограмме).
Вариант 2. В одном интервале сосредоточено п - 1 элемент, один элемент — в любом другом:
М2 = к ■ (к - 1).
Вариант 3. В одном интервале сосредоточено п - 2 элемента, два элемента сосредоточено в остальных, т. е. 1-й элемент попадает в к - 1 интервал, 2-й элемент — в к - 2 оставшихся, т. е.
М1 = к ■ (к -1),
М3 = к ■ (к-1 )(к-2).
Вариант 4. Аналогично:
М\ = к ■ (к -1),
М4 = к ■ (к-1 )(к-2),
М4 = к ■ (к-1 )(к-2)(к-3).
Вариант 5. По тому же принципу устанавливаем:
м\ = к ■ (к -1),
М52 = к ■ (к-1 )(к-2),
М53 = к ■ (к-1 )(к-2)(к-3),
М54 = к ■ (к-1 )(к-2)(к-3)(к-4),
и т. д.
Общее число таких теоретических дискретных распределений будет равно сумме
М1 + М2 + м1 + м2 + м\ + М 4 + +М 4 + м1 ... + М54.
Эта сумма равна числу разбиений
Rn целого числа n на к целых положительных слагаемых s,-, 0 < < n. Это
число конечно (хотя при больших n и велико), что в принципе позволяет запрограммировать моделирование этих разбиений. В процессе моделирования их можно посчитать, так что при таком подходе отсутствие формулы числа разбиений не препятствует решению задачи. Из разбиений необходимо отобрать мономодальные распределения числом Rn моно по принципу неумень-шения частот к некоторому /-тому интервалу и неувеличения их после этого интервала.
Таким образом получают все дискретные мономодальные распределения, возможные при определенных объемах выборки n и числа интервалов к. В общем случае каждое из этих распределений имеет одинаковую вероятность появления, равную P, = Rn моно^П = const. При некоторых специальных условиях можно рассматривать и неравные вероятности моделируемых распределений.
Специально рассмотрим вопрос о том, как понимать проверку составной гипотезы типа (1) на заданном уровне значимости а. Здесь может быть две трактовки. Первая из них заключается в том, что Н0 не отклоняется, если из всего многообразия мономодальных распределений {F,} найдется хотя бы одно, не отклоняемое на уровне значимости а, или, что то же самое, имеется хотя бы одно F/, принимаемое с доверительной вероятностью P = 1 - а. При таком понимании эмпирическое распределение не противоречит хотя бы одному из предложенных мономодальных распределений. Рассмотрим эту трактовку более подробно.
Если существует такое мономо-дальное распределение, которому не противоречит наше эмпирическое на уровне значимости а, то гипотеза мономодальности не отклоняется. Следовательно, принимать гипотезу би-или полимодальности нет оснований и имеющиеся "провалы" на эмпирической гистограмме следует считать несущественными. На языке математической логики это запишется следующим образом:
если
3(Яог ){P(Hoi = true) > 1- aKp },
/ =1.. M, (2)
где P — вероятность истинности i-той гипотезы Hoi, aKp — принятый уровень
значимости, / — номер теоретического мономодального распределения, М— общее число таких распределений, то распределение ^э мономодально.
По другому условие (2) можно записать следующим образом:
За, {а, > акр},
или
З/Х <х1 }
ыкр
Однако такая трактовка, на наш взгляд, не вполне соответствует существу задачи.
В другом, более приемлемом на наш взгляд варианте, гипотеза мономодальности может быть отклонена также и в том случае, если условие (2) по отдельности не выполняется ни для одного мономодального распределения, но её отклонение (или принятие) делается по совокупности полученных значений уровня значимости {а,}, .
Действительно, для неотклонения (принятия) гипотезы мономодальности необходимо, чтобы доверительная вероятность такого решения Р была бы больше критической, т. е.
р > р,
Kp
1 - акр;
Если проверяемая гипотеза составная, то необходимо, чтобы
Р(А1 + А2 + А3 + ... Ап) > Рк^ где А; - событие, что гипотеза Н0/ может быть принята с доверительной вероятностью Р, = 1 - а,. Для простоты рассмотрения ограничимся случаем, в котором гипотеза мономодальности состоит всего из двух простых гипотез Но1 и Но2. Тогда
Р(А1 + А2) = Р(А1) + Р(А2) - Р(АА) так как события А1 и А 2 — совместимые (наше эмпирическое распределение может не противоречить сразу нескольким мономодальным “теоретическим” распределениям).
Проиллюстрируем вторую трактовку проверки сложной составной гипотезы несколькими примерами, сведенными в табл. 1.
Проверка каждой из трех гистограмм проводилась по критерию х2. Результаты приведены в первых строках для каждой простой гипотезы Но1, Но2 и Ноз. В первом эмпирическом распределении доверительная вероятность принятия гипотез Но1 и Но2 мономодальности не достигает требуемого уровня. Однако по совокупности проверок обеих гипотез гипотеза мономодальности принимается (срабатывает вторая трактовка!).
Таблица 1
Пример последовательной проверки гипотезы мономодальности (критическое значение а= 0.1 (10%), число степеней свободы/= 10)
Номер теоретич. распреде- ления Гистограмма 1* Гистограмма 2 Гистограмма 3
1 х2 а Р х2 а Р х2 а Р
13.1 0.25 0.75 25.2 0.05 0.95 9.34 0.5 0.5
Л = 0.75 <7% = 0.9, гипотеза Н0\ не принимается Р\ = 0.95 > Ркр = 0.90, гипотеза мономодальности принимается на основании проверки Н01 Pi = 0.5 < Ркр = 0.9, гипотеза мономодальности на основании проверки Н0\ не принимается
2 13.4 0.20 0.80 Проверка гипотезы 2 не требуется 8.30 0.6 0.4
Р2 = 0.8 < Ркр = 0.9, гипотеза Н„2 не принимается Рг — 0.4 < Ркр = 0.9, гипотеза мономодальности на основании проверки Но2 не принимается
1+2 По совокупности проверок гипотез Н01 и Но2: Р1;2(Л! +А2) = 0.75 + 0.80- 0.750.80 = 1.55 - 0.60 = 0.95 Л,2 = 0.95 >Лр = 0.90; гипотеза мономодальности принимается по совокупности проверок двух простых гипотез Проверка по совокупности не требуется По совокупности проверок гипотез На\ и Но2: Л,2(^1+Л) = 0.5 + 0.4-0.40.5 = 0.9-0.2 = 0.7 Л;2 = 0.7 <Лр = 0.90; гипотеза мономодальности по совокупности проверок гипотез Н„ 1 и На2 не принимается
3 Проверка гипотезы 3 уже не требуется Проверка гипотезы 3 не требуется 9.30 0.5 0.5
Ръ = 0.5 < Ркр = 0.9, гипотеза мономодальности на основании проверки НоЪ не принимается
1+2 + 3 — — По совокупности проверок гипотез H0i, Н02 и Яо3: Р 1а2,зС^ 1 А2 + A3) = 0.7 + 0.5 — — 0.7-0.5 = 1.2 — 0.35 = 0.85; Л,2,3 = 0.85 <Лр = 0.9; гипотеза мономодальности по совокупности проверок гипотез Н0ъ Н02, Н0з не принимается
* Модельные примеры гистограмм
Во втором эмпирическом распределении мономодальность принимается уже по результатам проверки первой гипотезы Н01, поэтому дальнейшие проверки не требуются: они никаким образом не могут уменьшить доверительную вероятность. Тем не менее проверки других имеющихся простых гипотез, составляющих сложную нулевую гипотезу, если таковые имеются, могут быть продолжены для установления действительного значения доверительной вероятности, что в некоторых случаях бывает полезным или необходимым.
В третьем эмпирическом распределении гипотеза мономодальности не принята ни по результатам проверки отдельных гипотез Н01, Н02, Н03, ни по их совокупности. Если других простых гипотез мономодальности нет, то ре-
зультат проверки сложной гипотезы следует считать окончательным. Если существуют другие теоретические мо-номодальные распределения кроме испытанных Н01, Н02 и НоЪ, проверка должна быть продолжена либо до исчерпания Н0І, либо до достижения критического значения доверительной вероятности.
В связи с изложенными выше процедурами проверки сложных гипотез возникает следующая проблема. Она заключается в решении вопроса о том, следует ли учитывать вероятности появления в данной предметной области “образцовых” для сравнения мономо-дальных распределений, с которыми сравнивается наше эмпирическое. При проверке простой гипотезы, например гипотезы принадлежности эмпирического распределения гауссовскому, так-
же можно задаться вопросом, какова априорная вероятность встречаемости этого распределения в данной предметной области. Если она равна нулю, то независимо от результатов применения критерия согласия необходимо констатировать, что эмпирическое распределение в данной ситуации не может быть гауссовским. И наоборот, если гауссовское распределение здесь единственно возможное, то независимо от проверки принадлежности с необходимостью вытекает, что эмпирическое распределение - гауссовское. В общем случае
Р = Рд ■ РР, (5)
где Р - вероятность истинности принимаемой гипотезы, Рд - доверительная вероятность, полученная в результате проверки гипотезы (или установленная заранее), Рр - априорная вероятность распространения данного “теоретичес-
кого” распределения, испытуемого в качестве нулевой гипотезы.
В обычной практике проверки гипотез о законах распределений вероятность Рр не устанавливается и не учитывается, что соответствует условиям полной неизвестности априорных вероятностей распространенности видов распределений, когда для всех подбираемых для сравнения законов принимается Рр = const. Мы можем поступить так же, но в нашем случае объектов для сравнения с ними эмпирического распределения может быть много тысяч, и вопрос о том, какие из них вероятны в природе, не тривиален. Без учета распространенности тех или иных распределений описанная выше процедура должна квалифицироваться как подгонка. Ввод в такие процедуры априорных вероятностей представляется нам совершенно необходимым элементом.
Технические трудности применения критерия согласия х2 к дискретным распределениям выражаются в том, что во многих теоретических мономо-дальных распределениях часть интервалов будет иметь нулевые частоты, что недопустимо, так как они в формуле для расчета х2 появляются как в числителе, так и в знаменателе. Для преодоления этой трудности “нулевые” интервалы объединяются с соседними, чтобы ни в одном не оказалось нулевых частот. При этом аналогичным образом объединяются с суммированием частот те же интервалы эмпирической гистограммы. После такой фильтрации в эмпирической гистограмме могут оставаться интервалы с частотами меньше пяти, что при проверке допускать не рекомендуется (по другим источникам частоты должны быть не менее трех). Снова применяем фильтр с объединением интервалов эмпирической гистограммы, а также соответствующих интервалов теоретического распределения. После такой подготовки используем критерий х2.
Статистические таблицы, содержащие значения хи-квадрат распределения для различных степеней свободы и критические значения х2 для разных уровней значимости а, оказались непригодными для решения нашей задачи. Таблицы сделаны для дискретных значений уровня значимости, а нам необходимы точные значения а для получаемых в результате применения формулы (1) значений х2. В програм-
мах, реализующих описываемую методику, необходимо было предусмотреть способ расчета а для заданных степеней свободы и получаемых значений X2- Для этого в свою очередь необходимо было разработать ряд процедур вычисления специальных математических функций, и в их числе гамма-функцию — аналога факториала для дробного аргумента.
Разработанный нами программный комплекс на языке Паскаль-7 содержит ... программных единиц общим объемом ... операторов. Он способен обрабатывать гистограммы, содержащие до 25 интервалов с объемом выборки до 300. Обработка больших объемов выборки сильно удлиняет время счета, и без того немалое.
Методика испытана на небольших искусственных примерах, а также на реальном материале, любезно предоставленном Н. П. Юшкиным, которому автор обязан также постановкой задачи именно в широком плане, без предположений о виде теоретического распределения, как эталона для сравнения с эмпирическим. Исходные данные были заданы в виде гистограммы (см. рисунок). Необходимо было установить, какие из “провалов” гистограммы существенны без априорных предположений о виде распределения случайных величин в целой гистограмме и в отдельных её частях.
При анализе сложных гистограмм, содержащих несколько мод и разделяющих их “провалов”, предусмотрена
такая последовательность действий. Сначала рассчитывается значение X2 и а в целом для всей гистограммы. Если гипотеза мономодальности не отклоняется, анализ можно закончить. Если гистограмма разделяется на две части по интервалу самого существенного провала, то анализ продолжается аналогичным образом для каждой части. На практике “для надежности”, независимо от результатов проверки общей гистограммы или её части с “провалом”, анализ ведется до конца.
Для ускорения анализа можно воспользоваться вариантом программы, в котором не предусматривается моделирование мономодальных распределений, а используется ручной ввод частот мономодальной гистограммы, на глаз наименее отличающейся от эмпирической. Если при этом не превы-
сит критического значения а, эмпирическая гистограмма считается мономо-дальной. В противоположном случае частотами "теоретической" мономо-дальной гистограммы немного варьируют, но в пределах, оставляющих её мономодальной, анализ повторяется и ведется проверка по описанным схемам первой и второй трактовок проверки сложных гипотез.
В реальном примере исходные данные гистограммы (см. рисунок) представляют собой определения возраста цирконов из долины р. Оби изотопным методом (или по трекам осколков деления?). В этой гистограмме насчитыва-
S
и
(D
ч
(D
ч
(D
30
25
20
15
10
5
5
41
7
МІМІ
8
ib-
9
900 1000
78
0
100 200 300 400 500 600 700 800
1 2 3 4 5 6
Возраст, млн лет
Исходная гистограмма распределения цирконов из долины р. Оби по возрасту, млн лет. Цифры сверху - номера модальных интервалов, цифры снизу - номера провальных интервалов. Результаты анализа приведены в табл. 2
2
4
3
6
1
Таблица2
Результаты проверки на мономодальность исходной гистограммы и ее частей (критический уровень значимости акр = 0.1)
Анализируемая часть гистограммы Число интервалов к Число данных п Число степеней свободы / Значение X Уровень значимости а Вывод
Левая от провала 4 23 112 13 55.2 0.4 • 10-6 Провал 2 существен
Левая от начала до провала 2 13 77 7 1.03 0.996 Провал 1 несуществен
От провала 2 до провала 4 9 35 5 2.41 0.75 Провал 3 несуществен
От провала 4 до провала 6 6 8 1 0.67 0.42 Провал 5 несуществен
От провала 4 до правого конца гистограммы 21 17 2 10.3 0.005 Провал 6 существен
От провала 6 до провала 8 14 8 3 10.0 0.002 Провал 7 существен
От провала 7 до правого конца гистограммы 7 4 0 — — О существенности провала 8 судить невозможно ввиду отсутствия степеней свободы
ется девять мод. Визуально наиболее контрастным представляется провал между модами 2 и 3. Расчеты это подтверждают. Поэтому первым шагом является расчленение гистограммы по
провалу между модами 4 и 5 с целью анализа левой части гистограммы. Результат анализа этой части гистограммы на мономодальность приведен в первой строке табл. 2: она не моно-
модальна, и вероятность ошибки этого вывода не превысит 0.4 • 10-6, что ничтожно мало. Дальнейший ход анализа и его результаты приведены в этой же таблице.
Гравюра О. Велегжанинова “Бабочка на льду”
