Гипотеза независимости проверяется по критерию хи-квадрат для стандартного уровня α = 0.05 на сгруппированных данных. Группировка выполняется с помощью медианы, позволяющей разделить первичные данные на две равновероятные группы так, чтобы каждая клетка итоговой таблицы содержала не менее пяти наблюдений в соответствии с требованием критерия хи-квадрат.

Адекватность кубической модели сохраняется при удалении последнего элемента последовательности, что позволяет считать эту модель устойчивой к изменению длины последовательности, и, следовательно, кубическую модель можно использовать для экстраполяции найденной закономерности. Однако экстраполяция кубической модели на один шаг вперед не приводит к содержательному результату, поскольку ожидаемая дата 23-го землетрясения оказывается равной 1835, т.е. меньше реальной даты последнего землетрясения (см. рис. 3). Иначе говоря, кубическая модель является локально адекватным описанием исходных данных, которое не допускает расширения за пределы области наблюдений.

Обычно для улучшения статистической модели используются два способа: повышение степени полинома или увеличение объема выборки. В нашем случае, когда есть всего 22 уникальные даты, повышать степень полинома уже нельзя, но можно попытаться улучшить модель, используя весь массив исходных данных. Поскольку каждой уникальной дате соответствует своя частота появления в массиве исходных данных, то использование реальной частоты наблюдений позволяет увеличить размер выборки с 22 до 37. Однако оказалось, что использование частоты наблюдений в качестве дополнительного параметра не приводит к улучшению модели.

Можно получить еще один вариант модели, поставив в соответствие каждой уникальной дате размер локального массива, который содержит все наблюдения, не превосходящие рассматриваемую дату. Рассматривая число наблюдений в локальном массиве в качестве весовой характеристики, приходим к взвешенной регрессионной модели. Экстраполяция взвешенной модели позволяет получить ожидаемую дату следующего землетрясения в 1860 году, что больше конечной даты последнего землетрясения (1857), но не соответствует реальному состоянию на 2009 год. Повышая, без достаточных на то оснований, степень полинома взвешенной модели, находим, что ожидаемая дата следующего землетрясении начинает превышать 2009 год только при полиноме 5-ой степени. Столь высокая степень полинома свидетельствует о том, что возможности полиномиальной модели практически исчерпаны, и целесообразно перейти к моделям экспоненциального типа.

3.4. Кратные даты

Каждой наблюдаемой дате исходного массива данных соответствует определенная кратность (частота наблюдения), причем средняя кратность равна 37/22 = 1.68 (см. рис. 4).

Распределение кратности заметно отличается от нормального, а потому для проверки статистической значимости этой величины необходимо использовать свободные от формы распределения непараметрические критерии: критерий знаков или ранговый критерий Вилкоксона. Гипотеза о равенстве средней кратности единице отклоняется для стандартного уровня значимости α = 0.05, поскольку p-значение для обоих критериев равно 0.031, и, следовательно, частота наблюдений является существенным статистическим фактором.

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

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

4. Дата землетрясения как случайная величина

4. 1. Статистическое описание

Описание даты землетрясения как случайной величины, включая основные статистики и гистограмму, представлено на рисунке 5.

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

Для получения закона распределения даты землетрясений использовался весь набор стандартных распределений, предлагаемый в процедуре UNIVARIATE системы SAS [4], который включает следующие распределения: бета, экспоненциальное, гамма, логнормальное, нормальное и Вейбулла. Параметры стандартных распределений вычисляются по методу максимального правдоподобия, а качество подгонки оценивается по трем критериям согласия: Колмогорова-Смирнова, Крамера-фон Мизеса и Андерсона-Дарлинга. Проверка показала, что только три распределения согласуются с исходными данными на стандартном уровне значимости α = 0.05, а именно: нормальное, логнормальное и Вейбулла (см. рис. 5), причем эти описания настолько незначительно отличаются друг от друга, что для обоснованного выбора лучшей модели нужны дополнительные исследования.

4.2. Выбор закона распределения

Учитывая, что ожидаемая дата следующего землетрясения находится за пределами интервала наблюдения, для выбора лучшего распределения предлагается использовать устойчивость распределения, т.е. его способность оставаться допустимым описанием исходных данных при изменении объема выборки.

Выборку меньшего объема легко получить удалением максимального значения переменной из исходного массива данных до тех пор, пока это возможно. Для увеличения объема выборки к исходному массиву нужно добавить одно или несколько наблюдений. Одиночные наблюдения встречаются в исходном массиве 16 раз, а кратные только 6. Поскольку соотношение 16/6 не согласуется с гипотезой равных вероятностей 0.5/0.5 (критерий хи-квадрат, p-значение = 0.033), то наиболее вероятным расширением массива исходных данных является одиночное наблюдение, значение которого можно последовательно увеличивать до тех пор, пока это позволяют используемые критерии. Результаты расчетов сведены в таблицу 4, откуда следует, что наибольшей устойчивостью обладает логнормальное распределение (см. отмеченные ячейки), которое согласуется с исходными данными при изменении объема выборки от 35 до 38 в диапазоне от 1555 до 2054 года, т.е. на протяжении почти пятисот лет.

Таким образом, для прогностической модели выбор логнормального распределения становится предпочтительным.

4.3. Логарифмическое преобразование

В общем случае логнормальное распределение определяется тремя параметрами: порога (threshold), масштаба (scale) и формы (shape), которые подбирались методом максимального правдоподобия. Кроме того, известно, что логарифм случайной величины, распределенной по логнормальному закону, имеет нормальное распределение [6], а потому можно перейти к работе с нормальным распределением, используя преобразование x = log(date - theta), где theta - величина порога. Величина порога, найденная по методу максимального правдоподобия равна 423.86, однако выбор новой переменной с таким порогом не выдерживает проверки по критерию Шапиро-Уилка (p-значение = 0.0385), как наиболее мощному среди известных критериев нормальности [5]. Согласно расчетам, p-значение по критерию Шапиро-Уилка растет при уменьшении порога, достигая почти стационарного состояния при значениях порога близких к нулю (см. рис. 6).

Если в качестве эталона выбрать распределение, полученное по методу максимального правдоподобия, то сумма квадратов отклонений от этого эталона позволяет оценить целесообразность уменьшения порога. Как видно из рисунка 6, достаточно хорошим решением можно считать theta = 0, поскольку в этом случае распределение практически не отличается от оптимального, а p-значение по критерию нормальности Шапиро-Уилка равно 0.057.

Кумулятивная функция распределения логарифма даты представлена на рисунке 7, где для наглядности показаны её доверительные границы на основе D-статистики Колмогорова (D = 0.115, p-значение > 0.15).

Итак, массив исходных данных согласуется с нормальным распределением логарифма даты

где cdf - кумулятивная функция распределения [4], mean - среднее значение, а std - стандартное отклонение случайной величины X = {x=log(d_i), 1065 ≤ d_i ≤ 1857; i=1,..., 37}.

5. Прогнозирование

5.1. Расширение массива исходных данных

Проверка устойчивости позволяет установить, что согласие с нормальным распределением логарифма даты по критерию Шапиро-Уилка сохраняется (т.е. p-значение ≥ 0.05) при добавлении 38-го наблюдения с конечной датой вплоть до 2055 года. Это означает, что исходный массив допускает расширение в ограниченных пределах, а именно:

Статистическое описание расширенного массива с неизвестным значением d₃₈ зависит от среднего значения и стандартного отклонения случайной величины X₁, которые можно оценить, опираясь на опыт прошлого. В самом деле, если удалить из массива исходных данных конечную дату, то получим новый массив меньшего размера, у которого есть свои выборочные значения среднего и стандартного отклонения. Продолжая этот процесс до тех пор, пока это возможно, получим ретроспективную картину изменения выборочных характеристик в зависимости от объема выборки (см. рис. 8). Нетрудно заметить устойчивую тенденцию роста обеих характеристик, а высокие значения коэффициентов корреляции: r(mean, N) = 0.986 и r(std, N) = 0.890, свидетельствуют о возможности использования линейных регрессионных моделей для выявления статистических закономерностей.

5.2. Экстраполяция тенденций

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

Согласно проверке, эти модели являются адекватным описанием исходных данных, поскольку распределение остатков (разностей между фактическими и модельными значениями) подчиняются нормальному закону с нулевым средним, разброс остатков не выходит за пределы ±3σ, а величина остатков не зависит от объема выборки. Более того, увеличивая объем выборки до N = 38, можно проверить, что расширенная выборка сохраняет однородность по отношению к массиву исходных данных.

Таким образом, линейные регрессионные модели позволяют оценить неизвестные параметры расширенного массива:

Теперь следует убедиться в том, что существует допустимое распределение с такими параметрами. Последовательно изменяя величину d₃₈ с шагом в один год в пределах от 1858 до 2055, находим допустимые границы для среднего значения и стандартного отклонения логарифма даты, которые представлены в таблице 5.

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

5.3. Оптимальное приближение

Расширение массива в пределах допустимых ограничений не позволяет получить распределение с параметрами регрессионной модели

Однако распределение с регрессионными параметрами можно рассматривать как недостижимый эталон для оценки качества приближенного решения. Тогда среди допустимых распределений найдется такое, у которого сумма квадратов отклонений от эталона будет наименьшей, т.е.

5.4. Прогностическая модель

Рассмотрим расширенный массив данных, содержащий 2040 год в качестве ожидаемой даты следующего землетрясения, который описывается нормальным распределением со следующими параметрами:

Каждой наблюдаемой дате d соответствует квантиль q₁ = log(d) с эмпирической вероятностью p(q₁ = log(d)) = p (см. рис. 10), по которой можно найти гипотетический квантиль q₂ из нормального распределения случайной величины X₂

Затем, используя экспоненту для обратного преобразования, находим гипотетическую дату в годах

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

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

Таким образом, нормальное распределение логарифма даты представляет адекватную прогностическую модель крупных региональных землетрясений в районе разлома Сан-Анреас.

Максимальная относительная ошибка используемого описания достигает 9.1%, и, следовательно, предлагаемая модель относится к моделям с низкой точностью порядка 10%, что вполне естественно для столь малого массива исходных данных.

5.5. Оценка прогноза

Независимость ошибки моделирования от гипотетической даты позволяет перенести центр распределения ошибки в точку, соответствующую прогнозируемой дате. Нормальное распределение ошибки с центром в 2040 году, должно быть ограниченно допустимыми пределами моделирования 1858-2055, и, следовательно, становится усеченным нормальным распределением. Учитывая, что на момент завершения работы (2009) новое крупное землетрясение пока еще не произошло, усеченное нормальное распределение требует коррекции, соответствующей реальному пределу слева. Окончательный вариант распределения ожидаемой даты следующего землетрясения в пределах 2009-2055 представлен на рисунке 12.

Оказалось, что полученное распределение хорошо согласуется с равномерным распределением (критерий хи-квадрат: статистика = 0.2979, число степеней свободы = 8, p-значение = 1.0). Используя квантили равномерного распределения как массив гипотетических дат, находим, что полученный массив статистически однороден массиву исходных данных. Проверка однородности проводилась по трем непараметрическим критериям, поскольку оба массива далеки от нормального распределения (см. табл. 8).

Таким образом, ожидаемая дата следующего землетрясения распределена равномерно на отрезке между 2009 и 2055 годами. Доверительные пределы для стандартного уровня доверия в 95% соответственно равны 2010 и 2054, а доверительный интервал составляет 44 года. В абсолютном выражении эта величина сопоставима с жизнью двух человеческих поколений, а потому может показаться, что подобный прогноз не представляет никакой практической ценности. Однако, если сравнить результаты статистического прогноза с естественным уровнем точности, на котором находятся исходные данные, то получится совсем другая картина. В самом деле, стандартная ошибка среднего для исходного массива данных равна ±28.7 года, а половина доверительного интервала для прогнозируемой даты следующего землетрясения составляет 44/2=22 года. Иначе говоря, точность полученного прогноза вполне соизмерима с точностью описания реальных наблюдений, а потому можно считать, что прогностическая модель не вносит существенных искажений в описание исходных данных. Следовательно, с вероятностью 0.95 ожидаемая дата следующего землетрясения находится между 2010 и 2054 годами, причем любой год является одинаково вероятным.

6. Заключение

По результатам статистического анализа массива дат крупных землетрясений (больше 7 баллов по шкале Рихтера), наблюдаемых в районе разлома Сан-Анреас (Калифорния) можно сделать следующие выводы:
1. Даты крупных землетрясений, наблюдаемых в рассматриваемом регионе, подчиняются логнормальному распределению с нулевым порогом.
2. Логнормальное распределение даты региональных землетрясений сохраняется за пределами интервала наблюдения, допуская расширение массива наблюдений для прогнозирования даты следующего землетрясения.
3. Параметры расширенного распределения можно оценить на основе экстраполяции тенденций изменения среднего значения и стандартного отклонения логарифма даты в зависимости от объема выборки.
4. Ожидаемую дату следующего землетрясения определяет конечная дата расширенного распределения, которое представляет наилучшее приближение к распределению с эталонными параметрами, отражающими наиболее устойчивые тенденции.
5 Точность полученного прогноза вполне соизмерима с точностью отбора наблюдений, но явно недостаточна с практической точки зрения. Остается надеяться, что увеличение массива исходных данных позволит довести точность статистического моделирования до уровня современных потребностей.

Источники: