Практическая работа
по курсу «Математические модели окружающей среды»
Задано временное изменение уровня воды в некоторых пунктах запериод примерно в 170 лет.
Применить методы математической статистики для оценкихарактеристик и качества имеющихся данных наблюдений. Выполнить прогноз подъемауровня воды на будущее и проверить качество прогноза на уже имеющихся данных.
1. Рассчитатьмоменты ряда (среднее и среднеквадратичное значение), построить функциюраспределения и плотность функции распределения. Выполнить ее аппроксимациютеоретическими зависимостями.
/>
Рис. 1.1. Изменение уровня воды за период в 102 года
Минимальный уровень воды = 0.06328, максимальное значение уровня =0.6792
Заменим простой статистический ряд на статистический ряд сменьшим числом слагаемых, равным 100. И для такого ряда рассчитаем частотусобытия (в качестве события берем средний уровень воды).
Таким образом, имеем 100 интервалов, для каждого вычисляется частотасобытия (число событий в статистическом ряде, когда X = x, к общему числусобытий)
/>.
В нашем случае имеем N=1024 события, а m – число уровней,попавших i-ыйинтервал Очевидны свойства этой частоты
/>
Частоту различных уровней воды можно изобразить графически
/>
Рис. 1.2. График зависимости частоты от среднего уровняводы
Статистическая функция распределения есть «частота» события Хx в данном статистическом интервале
/>.
/>/>
/>Рис. 1.3. Функцияраспределения
/>Эта функция F*(x) являетсянеубывающей со следующими пределами:
/> F*(x ® –¥) = 0, F*(x ® + ¥) = 1.
С функцией распределения F(x) связана плотность функциираспределения f(x)
/>.
которая удовлетворяет следующим соотношениям:
f(x) ³ 0, ò f(x) dx = 1,
/>
Рис. 1.4. Плотность функции распределения
Была выполнена аппроксимация плотности функции распределениятеоретическими зависимостями: полиномами 6-ой, 9-ой, 15-ой степени,тригонометрическими многочленами. Оптимальным приближением оказался полином 9-ойстепени.
В качества критерия оптимальной аппроксимации использоваликритерий Пирсона
/>
/>
Рис. 1.5. Аппроксимация плотность функции распределенияполиномом 9-ой степени
Для нового ряда по имеющимся данным можно рассчитать математическоеожидание, характеризующее среднее значение уровня воды
/>,
и среднеквадратичное отклонение, характеризующее среднийразброс этих значений:
s*=/>.
где /> — дисперсия:
/>
xi – среднее значение случайной величинывнутри разряда.
В нашем случае, средний уровень воды равен 0.41, асреднеквадратичное отклонение – 0.119
2. В какойстепени данный ряд является стационарным? На каких временах данный ряд можносчитать стационарным? Дать оценки моментов для «кусков» ряда и построитьгистограммы оценок
Для того чтобы ряд был стационарным, должны быть выполненыусловия
- корреляционнаяфункция не зависит от времени
/>
математическое ожидание />
- дисперсия/>
-
Для проверки стационарности делим исходный ряд на />кусков, и для каждоготакого куска проверяем выполнение трех условий.
– Корреляционная функция.
Фиксируем />, где N – количество точек.
Считаем автокорреляционную функцию для первого отрезка, азатем – корреляционную функцию для каждых двух соседних кусков. Получаемзначение корреляционной функции при фиксированном />для каждого куска ряда.
Если процесс стационарный, то все значения должны совпадатьсо значением автокорреляционной функции.
/>
Рис. 2.1. Графики зависимости корреляционной функции отномера отрезка при различных />.
В качестве оценки корреляционной функции вычислилисреднеквадратичное отклонение от значения автокорреляционной функции.
/>
Рис. 2.2. Зависимость среднеквадратичного отклонения от />
– Математическое ожидание
Для каждого «куска» ряда вычисляется математическое ожидание />(/>; />). Затем находим среднее от/> и среднеквадратичноеотклонение s*=/>.
- Дисперсия
Для каждого «куска» ряда вычисляется дисперсия />(/>; />). Затем считаем среднее от/> и среднеквадратичноеотклонение s*=/>.
/>
Рис. 2.3. Зависимости среднеквадратичного отклонения D(M) и D(D) от />.
В нашем случае, критерием стационарности является минимумсреднеквадратичного отклонения от значения автокорреляционной функции, минимум /> и />.
Этому условию удовлетворяет «кусок» ряда, длиной />.
Таким образом, исходный ряд стационарен на периоде T=21 год.
1. Оценка математического ожидания
Проверяем
Ø состоятельность оценки для каждогостационарного «куска» ряда.
/>,
где /> – математическоеожидание на стационарном периоде
/> – среднее значение в зависимости от числа данных
/>/>
/>/>
/>/>
/>/>
Ø несмещенностьоценки
M [a*] = a,
Номер интервала
1
2
3
4
5
6
7
8
M [a*] 0,40938 0,41218 0,41058 0,41152 0,40758 0,41118 0,41259 0,40985
a 0,40714 0,40661 0,40437 0,4080 0,40492 0,40906 0,41206 0,41018
2. Оценка дисперсии
Проверяем
Ø состоятельность оценки для каждогостационарного «куска» ряда.
/>,
где /> – среднеквадратичноеотклонение на стационарном периоде
/> — среднеквадратичное отклонение взависимости от числа данных
/>/>
/>/>
/>/>
/>/>
Ø несмещенностьоценки
M [a*] = a,
Номер интервала
1
2
3
4
5
6
7
8
M [a*] 0,11862 0,11507 0,11944 0,1235 0,12227 0,11891 0,11709 0,1185
a 0,11477 0,11391 0,12122 0,11959 0,11959 0,11674 0,12163 0,11842
3. Вычислить куммулятивную частоту превышения уровня и сделатьоценки ее стационарности
/> – куммулятивная частота превышения уровняа.
/>
Рис. 3.1.Зависимость куммулятивной частоты от уровня превышения.
Разбиваем исходный ряд на на /> отрезков,для каждого «куска» ряда строим функции зависимости куммулятивной частоты отуровня превышения, и оцениваем стационарность полученных
n-зависимостей критериям, полученным в п. 2.
Анализируя оценки корреляционной функции, математическогоожидания и дисперсии, находим период стационарности куммулятивной частотыпревышения уровня T=17.85 лет.
/>
Рис. 3.2. Оценка корреляционной функции
/>
Рис. 3.3. Оценки математического ожидания и дисперсии.
4. С помощью пуассоновской статистики дать долгосрочный прогнозпревышения уровня. На какой срок Вы гарантируете такой прогноз?
/>
– искомая формула вероятности того, что за периодвремени T произойдет ровно m превышений уровня a. Вероятностьвозникновения таких ситуаций определяется средней частотой превышения уровня a и временем прогноза. Изанализа данных за 176 лет на стационарном периоде находим функцию зависимости среднейчастоты от уровня a/>.
/>
Рис. 4.1. Средняя частота превышения уровня
Нас интересуют редкие события, например, превышение уровня />. Соответственно, среднюючастоту превышения такого уровня можно определить из графика />, и она равна/>0.0635
Зная среднюю частоту, теперь можно вычислить вероятность того, чтоза период, например, равный 220 лет, произойдет ровно 1,2,3,4.
/>
Рис. 4.2. Вероятность возникновения ровно m – аварий
Зависимости вероятности от прогнозируемого времени дляразного числа превышения уровня (m=1,2,3,4) являются немонотонными, и их максимумприходится на моменты времени />
Как видно изграфика, сначала более вероятным является только одно превышение уровня a=0.6, затем два, три…
Например, вероятность того, что за 62 года произойдет трипревышения уровня a=0.6, самая высокая и равна p= 0,19775.
5. Рассчитайте среднюю частоту появления выбросов и среднее времявыброса
Будем рассматривать выбросы (превышение уровня, например, a=0.6) за период, накотором исходный ряд стационарен.
Выброс характеризуется следующим условием
/>
Тогда можем найти среднюю частоту выброса за уровень а=0.6/>=0.063
Средняя продолжительность выброса для стационарных случайныхпроцессов
/>
Отсюда, зная функцию распределения случайной величины, порог,начиная с которого процесс ведет к катастрофе, и период прогноза, можнорассчитать среднее время катастрофической ситуации, в нашем случае – превышениезаданного уровня.
Среднее число выбросов за период T для стационарногопроцесса определяется:
/>
где P – среднюю частоту выброса
Таким образом, среднее число катастроф пропорциональнопродолжительности времени прогноза и падает с увеличением порогового значения,определяющего возникновение катастрофической ситуации.
Средняя продолжительность выброса может быть вычислена поформуле
/>
и она не зависит от прогнозируемого времени (для стационарныхпроцессов
Средняя продолжительность выброса за уровень а=0.6равна /> лет.
6. Попытайтесь сделать краткосрочный прогноз уровня воды,используя линейный и корреляционный анализ. Проверьте его на уже имеющихсяданных
1. Линейный прогноз.
Рассматриваем наш ряд на стационарном периоде T. Выбираем n=127 точек. Тогда
/>
Нам известны значения />.Решая систему уравнений, находим коэффициенты />.
Тогда можем спрогнозировать любую точку ряда, например, />:
/>
В качестве проверки найденных коэффициентов сделали прогноз с22-ого по 33-ий год (рис. 6.1) и c 33-ого по 44-ий год (рис. 6.2).
/>
Рис. 6.1.Линейный прогноз с 22-ого по 33-ий год
/>
Рис. 6.2.Линейный прогноз с 33-ого по 44-ый год
2.Корреляционный прогноз
Рассматриваем наш ряд на стационарном периоде T. Тогда можем представитьв виде
/>
где
/>
Зная коэффициенты /> и /> можно спрогнозироватьуровень воды, например, на 33 года.
В качестве проверки найденных коэффициентов сделали прогнозна первые 22 года (рис. 6.3) и с 22-ого по 33-ий год (рис. 6.4). Каквидно из графиков. Результаты совпадают в пределах погрешности.
/>
Рис. 6.3. Корреляционный прогноз на первые 22 года
/>
Рис. 6.4. Корреляционный прогноз с 22-ого по 33-ий год
Анализируя полученные результаты и используя уже известныеданные, можно сказать, что корреляционный прогноз более точный.