Известна ли принадлежность объектов к группам?
Нужно ли искать диагностические значения признаков?
Деревья классификации да
SVM,
Random Forest,
Дискриминантый анализ нет да
Предполагается ли иерархия групп?
Иерархический кластерный анализ да
Данные —
это матрица расстояний?
Многомерное шкалирование да
Анализ главных компонент нет нет нет
Рис. 51. Как выбрать подходящий метод многомерного анализа ных матрицу расстояний), и тогда станут доступны и другие методы анализа.
* * *
Ответ к задаче про иерархическую классификацию пива
Чтобы построить иерархическую классификацию, надо вначале вычис- лить матрицу расстояний, а уже потом сгруппировать объекты и изоб- разить это при помощи дендрограммы. Сначала посмотрим на данные:
> pivo <- read.table("data/pivo-s.txt", h=TRUE)
> head(pivo)
C.1 C.2 C.3 C.4 C.5 C.6 C.7 C.8 C.9 C.10
baltika.3 0
1 0
0 0
0 1
0 1
1
baltika.6 0
1 1
1 1
0 1
0 1
1
baltika.9 0
1 1
1 1
0 1
0 1
1
ochak.klassicheskoe
0 1
0 0
0 0
0 0
1 1
ochak.temnoe
0 1
1 1
1 0
0 0
1 1
sib.kor.prazdnichnoe
1 1
1 1
1 1
0 0
0 0
Классификация с обучением,
или Дискриминантный анализ
173
Данные, как видите, бинарные, и, значит, надо применять специфи- ческий метод вычисления расстояний. (Кстати говоря, бинарные дан- ные часто используются там, где нужно классифицировать объекты,
описанные качественными характеристиками). Воспользуемся функци- ей vegdist() из пакета vegan, которая умеет использовать популярный для бинарных данных метод Жаккара (method="jaccard"). Наверное,
до этого момента вы не знали про такой метод и могли решить зада- чу несколько иначе, например, используя метод binary из стандартной функции dist(). Вот так мы получили дендрограмму, изображенную на рис. 52: 0
> library(vegan)
> pivo.d <- vegdist(pivo, "jaccard")
> plot(hclust(pivo.d, "ward.D"), main="", xlab="", sub="")
ve na
.p or te r
si ne br yu xo v
so ko l.b ez al ko go ljn oe so ko l.s ve tl oe bu dw ei se r
si b.
ko r.
or ig in al jn oe ef es he in ek en bo ch ka re v.
sv et lo e
kl in sk oe
.z ol ot oe kl in sk oe
.t em no e
ve na
.p et er go f
st ar yj
.m el jn ik
.s ve tl oe zh ig ul ev sk oe xa m
ov ni ch es ko e
tu ljs ko e.
kr ep ko e
tu ljs ko e.
ar se na ljn oe tu ljs ko e.
or ig in al jn oe ba lt ik a.
3
oc ha k.
kl as si ch es ko e
af an as ij.
sv et lo e
oc ha k.
sp ec ia ljn oe ba lt ik a.
6
ba lt ik a.
9
oc ha k.
te m
no e
si b.
ko r.
pr az dn ic hn oe
0.
0 0.
5 1.
0 1.
5
Рис. 52. Иерархическая классификация сортов пива
Получились две большие группы, одна с «Жигулевским» и «Балти- кой», другая с «Клинским» и «Старым мельником».
174
Анализ структуры: data mining
Ответ к задаче про два вида растений
. Попробуем построить дерево классификации (рис. 53):
> eq <- read.table("data/eq.txt", h=TRUE, stringsAsFactors=TRUE)
> eq.tree <- tree(eq[,1] ., eq[,-1])
> plot(eq.tree); text(eq.tree)
|
N.REB < 15.5
DL.R < 631.5
arvense arvense fluviatile
Рис. 53. Дерево классификации показывает, что два вида хвощей раз- личаются по признаку N.REB (это количество ребер стебля)
Глава 7
Узнаем будущее: анализ временных рядов
В этой главе мы рассмотрим только самые основные принципы ра- боты с временн
0
ыми рядами. За более подробными сведениями рекомен- дуем обратиться к специальной литературе.
7.1. Что такое временные ряды
Во многих областях деятельности людей замеры показателей прово- дятся не один раз, а повторяются через некоторые интервалы времени.
Иногда
этот интервал равен многим годам, как при переписи населения страны, иногда — дням, часам, минутам и даже секундам, но интервал между измерениями во временном ряду есть всегда. Его называют ин- тервалом выборки
(sampling interval). А образующийся в результате выборки ряд данных называют временн
0
ым рядом
(time series).
В любом временном ряду можно выделить две компоненты:
1) неслучайную (детерминированную) компоненту;
2) случайную компоненту.
Неслучайная компонента обычно наиболее интересна, так как она дает возможность проверить гипотезы о производящем временной ряд явлении. Математическая модель неслучайной компоненты может быть использована для прогноза поведения временного ряда в будущем.
7.2. Тренд и период колебаний
Если явление, результатом которого является изучаемый ряд, за- висит от времени года (или времени суток, или дня недели, или ино- го фиксированного периода календаря), то из неслучайной компоненты может быть выделена еще одна компонента — сезонные колебания явле- ния. Ее следует отличать от циклической компоненты, не привязанной к какому-либо естественному календарному циклу.
176
Узнаем будущее: анализ временных рядов
Под трендом (тенденцией) понимают неслучайную и непериодиче- скую компоненту ряда. Первый вопрос, с которым сталкивается иссле- дователь, анализирующий временной ряд,— существует ли в нем тренд?
При наличии во временн
0
ом ряде тренда и периодической компонен- ты значения любого последующего значения ряда зависят от преды- дущих. Силу и знак этой связи можно измерить уже знакомым вам инструментом — коэффициентом корреляции. Корреляционная зави- симость между последовательными значениями временного ряда назы- вается автокорреляцией.
Коэффициент автокорреляции первого порядка определяет зависи- мость между соседними значениями ряда t n
и t n−1
, б
0
ольших порядков —
между более отдаленными значениями. Лаг (сдвиг) автокорреляции —
это количество периодов временного ряда, между которыми опреде- ляется коэффициент автокорреляции. Последовательность коэффици- ентов автокорреляции первого, второго и других порядков называется автокорреляционной функцией временного ряда.
Анализ автокорреляционной
функции позволяет найти лаг, при ко- тором автокорреляция наиболее высокая, а следовательно, связь между текущим и предыдущими уровнями временного ряда наиболее тесная.
Если значимым оказался только первый коэффициент автокорреляции
(коэффициент автокорреляции первого порядка), временной ряд, ско- рее всего, содержит только тенденцию (тренд). Если значимым ока- зался коэффициент автокорреляции, соответствующий лагу n, то ряд содержит циклические колебания с периодичностью в n моментов вре- мени. Если ни один из коэффициентов автокорреляции не является зна- чимым, то можно сказать, что либо ряд не содержит тенденции (тренда)
и циклических колебаний, либо ряд содержит нелинейную тенденцию,
которую линейный коэффициент корреляции выявить не способен.
Взаимная корреляция (кросс-корреляция) отражает, есть ли связь между рядами. В этом случае расчет происходит так же, как и для ав- токорреляции, только коэффициент корреляции рассчитывается между основным рядом и рядом, связь с которым основного ряда нужно опре- делить. Лаг (сдвиг) при этом может быть и отрицательной величиной,
поскольку цель расчета взаимной корреляции — выяснение того, какой из двух рядов «ведущий».
7.3. Построение временного ряда
Анализ временн
0
ого ряда часто строится вокруг объяснения выяв- ленного тренда и циклических колебаний значений ряда в рамках неко- торой статистической модели.
Построение временного ряда
177
Найденная модель позволяет: прогнозировать будущие значения ря- да (forecasting), генерировать искусственный временной ряд, все стати- стические характеристики которого эквивалентны исходному (simula- tion), и заполнять пробелы в исходном временном ряду наиболее веро- ятными значениями.
Нужно отличать экстраполяцию временного ряда (прогноз буду- щих значений ряда) от интерполяции (заполнение пробелов между имеющимися данными ряда). Не всегда модели ряда, пригодные для интерполяции, можно использовать для прогноза. Например, полином
(степенн
0
ое уравнение с коэффициентами) очень хорошо сглаживает ис- ходный ряд значений и позволяет получить оценку показателя, который описывает данный ряд в промежутках между значениями ряда. Но ес- ли мы попытаемся продлить полином за стартовое значение ряда, то получим совершенно случайный результат.
Вместе с тем обычный ли- нейный тренд, хотя и не столь изощренно следует изгибам внутри ряда,
дает устойчивый прогноз развития ряда в будущем.
Разные участки ряда могут описывать различные статистические модели, в этом случае говорят, что ряд нестационарный. Нестационар- ный временной ряд во многих случаях удается превратить в стационар- ный путем преобразования данных.
В базовые возможности R входят средства для представления и ана- лиза временн
0
ых рядов. Основным типом временн
0
ых данных является
«ts», который представляет собой временной ряд, состоящий из значе- ний, разделенных одинаковыми интервалами времени. Временн
0
ые ряды могут быть образованы и неравномерно отстоящими друг от друга зна- чениями. В этом случае следует воспользоваться специальными типами данных — zoo и its, которые становятся доступными после загрузки пакетов с теми же именами.
Часто необходимо обрабатывать календарные даты. По умолчанию read.table()
считывает все нечисловые даты (например, «12/15/04»
или «2004-12-15») как текстовые строки или факторы. Поэтому после загрузки таких данных при помощи read.table() нужно обязатель- но применить функцию as.Date(). Она «понимает» описание шаблона даты и преобразует строки символов в тип данных Date. В последних версиях R она работает и с факторами:
> dates.df <- data.frame(dates=c("2011-01-01","2011-01-02",
+ "2011-01-03","2011-01-04","2011-01-05"))
> str(dates.df$dates)
chr [1:5] "2011-01-01" "2011-01-02" "2011-01-03" ...
> dates.1 <- as.Date(dates.df$dates, "%Y-%m-%d")
> str(dates.1)
Date[1:5], format: "2011-01-01" "2011-01-02" "2011-01-03"
178
Узнаем будущее: анализ временных рядов "2011-01-04" ...
А вот как создаются временные ряды типа ts:
> ts(1:10,
# ряд данных
+ frequency = 4,
# поквартально
+ start = c(1959, 2)) # начинаем во втором квартале 1959 года
Qtr1 Qtr2 Qtr3 Qtr4 1959 1
2 3
1960 4
5 6
7 1961 8
9 10
Можно конвертировать сразу матрицу, тогда каждая колонка мат- рицы станет отдельным временн
0
ым рядом:
# матрица данных из трех столбцов
> z <- ts(matrix(rnorm(300), 100, 3),
+ start=c(1961, 1), # начинаем в 1-й месяц 1961 года
+ frequency=12)
# помесячно class(z)
[1] "mts" "ts"
# тип данных -- многомерный временной ряд
Ряды отображаются графически с помощью стандартной функции plot()
(рис. 54):
> plot(z,
+ plot.type="single", # поместить все ряды на одном графике
+ lty=1:3)
# типы линий временных рядов
Методы для анализа временн
0
ых рядов и их моделирования включа- ют ARIMA-модели, реализованные в функциях arima(), AR() и VAR(),
структурные модели в StructTS(), функции автокорреляции и частной автокорреляции в acf() и pacf(), классическую декомпозицию вре- менного ряда в decompose(), STL-декомпозицию в stl(), скользящее среднее и авторегрессивный фильтр в filter().
Покажем на примере, как применить некоторые из этих функций.
В текстовом файле leaf2-4.txt в директории data записаны результа- ты длившихся трое суток непрерывных наблюдений над хищным рас- тением росянкой. Листья этого растения постоянно открываются и за- крываются «в надежде» поймать и затем переварить мелкое насекомое.
Файл содержит результаты наблюдений над четвертым листом второго растения, поэтому он так называется. Состояние листа отмечали каж-
Построение временного ряда
179 1962 1964 1966 1968
−
2
−
1 0
1 2
Рис. 54.
График трех временн0
ых рядов с общими осями времени дые 40 минут, всего в сутки делали 36 наблюдений. Попробуем сделать временной ряд из колонки FORM, в которой закодированы изменения формы пластинки листа (1 — практически плоская, 2 — вогнутая); это шкальные данные, поскольку можно представить себе форму = 1.5.
Сначала посмотрим, как устроен файл данных, при помощи команды file.show("data/leaf2-4.txt")
:
K.UVL;FORM;ZAGN;OTOGN;SVEZH;POLUPER;OSTATKI
2;1;2;2;1;0;1 1;1;2;1;0;1;1 1;1;2;1;1;0;1 2;2;3;1;1;0;0 1;2;3;1;1;0;0
Теперь можно загружать его:
> leaf <- read.table("data/leaf2-4.txt", head=TRUE,
180
Узнаем будущее: анализ временных рядов
+ as.is=TRUE, sep=";")
Сразу ст
0
оит посмотреть, все ли в порядке:
> str(leaf)
’data.frame’: 80 obs. of
7 variables:
$ K.UVL
: int
2 1 1 2 1 1 1 1 1 1 ...
$ FORM
: int
1 1 1 2 2 2 2 2 2 2 ...
$ ZAGN
: int
2 2 2 3 3 3 2 2 2 2 ...
$ OTOGN
: int
2 1 1 1 1 1 1 1 1 1 ...
$ SVEZH
: int
1 0 1 1 1 1 1 1 1 1 ...
$ POLUPER: int
0 1 0 0 0 0 0 0 0 0 ...
$ OSTATKI: int
1 1 1 0 0 0 1 1 1 1 ...
> summary(leaf)
K.UVL
FORM
ZAGN
Min.
:1.000
Min.
:1.0
Min.
:1.0 1st Qu.:1.000 1st Qu.:1.0 1st Qu.:2.0
Median :1.000
Median :2.0
Median :2.0
Mean
:1.325
Mean
:1.7
Mean
:2.5 3rd Qu.:2.000 3rd Qu.:2.0 3rd Qu.:3.0
Max.
:2.000
Max.
:2.0
Max.
:4.0
OTOGN
SVEZH
POLUPER
Min.
:0.00
Min.
:0.0000
Min.
:0.000 1st Qu.:0.00 1st Qu.:0.0000 1st Qu.:1.000
Median :1.00
Median :0.0000
Median :1.000
Mean
:0.75
Mean
:0.1625
Mean
:0.925 3rd Qu.:1.00 3rd Qu.:0.0000 3rd Qu.:1.000
Max.
:2.00
Max.
:1.0000
Max.
:2.000
OSTATKI
Min.
:0.0000 1st Qu.:0.0000
Median :1.0000
Mean
:0.6125 3rd Qu.:1.0000
Max.
:2.0000
Все загрузилось правильно, видимых ошибок и выбросов нет. Теперь преобразуем колонку FORM во временной ряд:
> forma <- ts(leaf$FORM, frequency=36)
Проверим:
Построение временного ряда
181
> str(forma)
Time-Series [1:80] from 1 to 3.19: 1 1 1 2 2 2 2 2 2 2 ...
Все правильно, наблюдения велись чуть больше трех (3.19) суток.
Ну вот, а теперь попробуем понять, насколько наши данные периодичны и есть ли в них тренд (рис. 55):
> (acf(forma, main=""))
Autocorrelations of series ‘forma’, by lag
0.0000 0.0278 0.0556 0.0833 0.1111 0.1389 0.1667 0.1944 0.2222 1.000 0.614 0.287 0.079 0.074 0.068 -0.038 -0.085 -0.132 0.2500 0.2778
-0.137 -0.024 0.3056 0.3333 0.3611 0.3889 0.4167 0.4444 0.4722 0.5000 0.5278 0.030 0.025 -0.082 -0.087 0.027 0.021 -0.043 -0.090 -0.137 0.0 0.1 0.2 0.3 0.4 0.5
−
0.
2 0.
0 0.
2 0.
4 0.
6 0.
8 1.
0
Lag
A
C
F
Рис. 55. График автокорреляций для состояния формы листа росянки
182
Узнаем будущее: анализ временных рядов
Эта команда («auto-correlation function», ACF) выводит коэффици- енты автокорреляции и рисует график автокорреляции, на котором в нашем случае можно увидеть, что значимой периодичности нет — все пики лежат внутри обозначенного пунктиром доверительного интерва- ла, за исключением самых первых пиков, которые соответствуют авто- корреляции без лага или с очень маленьким лагом. По сути, это пока- зывает, что в пределах 0.05 суток (у нас период наблюдений — сутки),
то есть около 1 часа, следующее состояние листа будет таким же, как и текущее, а вот на б
0
ольших интервалах таких предсказаний сделать нельзя. То, что волнообразный график пиков как бы затухает, говорит о том, что в наших данных возможен тренд. Проверим это (рис. 56):
> plot(stl(forma, s.window="periodic")$time.series, main="")
−
0.
4
−
0.
2 0.
0 0.
2
se as on al
1.
7 1.
8 1.
9 2.
0
tr en d
−
0.
6
−
0.
2 0.
2 0.
6 1.0 1.5 2.0 2.5 3.0
re m
ai nd er
Time
Рис. 56. График сезонной декомпозиции для состояния формы листа росянки. Возможный тренд изображен на среднем графике
Действительно, наблюдается тенденция к уменьшению значения фор- мы с течением времени. Мы выяснили это при помощи функции stl()
Прогноз
183
(названа по имени метода, STL — «Seasonal Decomposition of Time Se- ries by Loess»), которая вычленяет из временного ряда три компоненты:
сезонную (в данном случае суточную), тренд и случайную, при помощи сглаживания данных методом LOESS.
Задача
. Попробуйте понять, имеет ли другой признак того же листа
(K.UVL, коэффициент увлажнения, отражающий степень «мокрости»
листа) такую же периодичность и тренд.
7.4. Прогноз
Научившись
базовым манипуляциям с временн0
ыми рядами, мы мо- жем попробовать решить задачу построения модели временного ряда.
Построенная модель позволит нам проследить развитие наблюдаемого процесса в будущем. Кроме того, мы сможем рассмотреть применение более сложных функций анализа временн
0
ых рядов, а заодно познако- мить читателя с общими принципами сравнительного анализа статис- тических моделей.
Наш пример будет заключаться в прогнозе числа абонентов интернет- провайдера. Исходные данные состоят из:
1) данных о подключениях за декабрь 2004 года;
2) помесячных данных о подключениях в 2005–2008 годах.
> polzovateli <- ts(read.table("data/data.txt")$V3,
+ start=c(2004,12), frequency=12)
Общее количество абонентов на каждый месяц отчетного периода составило:
> cum.polzovateli <- ts(cumsum(polzovateli),
+ start=c(2004,12), frequency=12)
Отобразим данные помесячного подключения и изменения коли- чества пользователей графически. Оба временн
0
ых ряда показывают экспоненциальный рост, поэтому выведем их в полулогарифмических координатах (рис. 57):
> oldpar <- par(mfrow=c(2,1))
> plot(polzovateli, type="b", log="y", xlab="")
> plot(cum.polzovateli, type="b", ylim=c(1,3000), log="y")
> par(oldpar)
184
Узнаем будущее: анализ временных рядов
2005 2006 2007 2008 2009 10 20 50 10 0
20 0
2005 2006 2007 2008 2009 1
5 50 50 0
Рис. 57. Изменения помесячного и общего количества пользователей
Линейный рост временн
0
ых рядов в полулогарифмических коорди- натах подтверждает предположение об экспоненциальном росте во вре- мени и количества подключений в месяц, и общего количества пользо- вателей.
Попробуем теперь построить модель временного ряда общего числа подключений распространенным методом ARIMA («Autoregressive In- tegrated Moving Average», авторегрессия интегрированного скользяще- го среднего). Нам надо выбрать подходящее значение параметра order,
для этого мы будем перебирать различные значения каждого из трех его компонентов:
> model01 <- arima(cum.polzovateli, order=c(0,0,1))
> model02 <- arima(cum.polzovateli, order=c(0,0,2))
> model03 <- arima(cum.polzovateli, order=c(0,0,3))
> model04 <- arima(cum.polzovateli, order=c(0,0,4))
> model05 <- arima(cum.polzovateli, order=c(0,0,5))
> model06 <- arima(cum.polzovateli, order=c(0,0,6))
> model07 <- arima(cum.polzovateli, order=c(0,0,7))
Прогноз
185
> model08 <- arima(cum.polzovateli, order=c(0,0,8))
> model09 <- arima(cum.polzovateli, order=c(0,0,9))
> model010 <- arima(cum.polzovateli, order=c(0,0,10))
> model011 <- arima(cum.polzovateli, order=c(0,0,11))
> model012 <- arima(cum.polzovateli, order=c(0,0,12))
> model013 <- arima(cum.polzovateli, order=c(0,0,13))
> model014 <- arima(cum.polzovateli, order=c(0,0,14))
(Можно было бы написать здесь цикл с оператором for, как мы де- лали для данных об отравлении в главе про двумерные данные, но нам было лень. Так что вот вам, дорогой читатель, задача: напишите такой цикл. Важная подсказка: нужно использовать функцию assign().)
Теперь сравним модели. «Лучшая» модель будет соответствовать минимуму AIC (рис. 58):
> plot(AIC(model01, model02, model03, model04, model05, model06,
+ model07, model08, model09, model010, model011, model012,
+ model013, model014), type="b")
На графике можно увидеть, что по ходу кривой первый минимум наблюдается в районе компонента, соответствующего коэффициенту =
12, а дальше ход вычислений становится нестабильным.
Теперь выберем лаг авторегресии (первый элемент order):
> model012 <- arima(cum.polzovateli, order=c(0,0,12))
> model112 <- arima(cum.polzovateli, order=c(1,0,12))
> model212 <- arima(cum.polzovateli, order=c(2,0,12))
> model312 <- arima(cum.polzovateli, order=c(3,0,12))
> model412 <- arima(cum.polzovateli, order=c(4,0,12))
Здесь тоже можно применить AIC:
> AIC(model012, model112, model212, model312, model412)
df
AIC
model012 14 477.4036
model112 15 438.3171
model212 16 435.7380
model312 17 438.1186
model412 18 439.0120
AIC минимален, когда лаг авторегресии равен 2, поэтому принимаем это значение для дальнейшего анализа.
Аналогично выберем второй компонент:
186
Узнаем будущее: анализ временных рядов
4 6
8 10 12 14 16 50 0
55 0
60 0
65 0
70 0
75 0
df
A
IC
Рис. 58. Поиск наилучшей модели по минимуму значения AIC
> model2120 <- arima(cum.polzovateli, order=c(2,0,12))
> model2121 <- arima(cum.polzovateli, order=c(2,1,12))
> model2122 <- arima(cum.polzovateli, order=c(2,2,12))
> model2123 <- arima(cum.polzovateli, order=c(2,3,12))
> model2124 <- arima(cum.polzovateli, order=c(2,4,12))
> model2125 <- arima(cum.polzovateli, order=c(2,5,12))
> AIC(model2120, model2121 ,model2122, model2123, model2124)
df
AIC
model2120 16 435.7380
model2121 15 421.6246
model2122 15 405.5416
model2123 15 399.1918
model2124 15 407.6942
Видно, что оптимальной моделью является model2123.
Ну а теперь, зная оптимальную модель, построим прогноз изме- нения общего числа абонентов на 2009 год (рис. 59):
Прогноз
187
> plot(cum.polzovateli, # Накопленное число пользователей
+ xlim=c(2004.7,2010), ylim=c(0,6500))
# Границы включают и данные прогноза. Линия прогноза:
> lines(predict(model2123, n.ahead=12, se.fit = TRUE)$pred,
+ col="green")
# Верхняя граница прогноза:
> lines(predict(model2123, n.ahead=12, se.fit = TRUE)$se +
+ predict(model2123, n.ahead=12, se.fit = TRUE)$pred,
+ col="red")
# Нижняя граница прогноза:
> lines(-predict(model2123, n.ahead=12, se.fit = TRUE)$se +
+ predict(model2123, n.ahead=12, se.fit = TRUE)$pred,
+ col="red")
2005 2006 2007 2008 2009 2010 0
10 00 20 00 30 00 40 00 50 00 60 00
Рис. 59. Прогноз изменения общего числа абонентов на 2009 год
Максимальное и минимальное ожидаемое количество абонентов по месяцам 2009 года составит:
> round(predict(model2123,
188
Узнаем будущее: анализ временных рядов
+ n.ahead=12,
# период прогноза
+ se.fit = TRUE)$se +
# берем из объекта только ошибку
+ predict(model2123,
# складываем ошибку и данные прогноза
+ n.ahead=12,
+ se.fit = TRUE)$pred) #
берем только данные самого прогнозаJan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec
2009 3158 3337 3533 3767 3992 4239 4494 4765 5050 5349 5663 5991
> round(-predict(model2123, n.ahead=12, se.fit = TRUE)$se +
+ predict(model2123, n.ahead=12, se.fit = TRUE)$pred)
Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec
2009 3134 3283 3437 3628 3817 4031 4253 4490 4728 4969 5213 5463
Что же, теперь можно заняться обещанным в предисловии «пред- сказанием курса доллара на следующую неделю». Итак, задача: в фай- ле dollar.txt содержатся значения курса доллара Центрального Бан- ка с 1 июля по 9 августа 2011 года, всего за 11 недель. Попробуйте предсказать курс доллара на две недели вперед. Чтобы проверить эф- фективность предсказания, возьмите для модели данные по 26 июля, а предскажите последние две недели.
* * *
Ответ к задаче про лист росянки
. Применим тот же подход,
который мы использовали для признака формы листа:
> uvl <- ts(leaf$K.UVL, frequency=36)
> str(uvl)
Time-Series [1:80] from 1 to 3.19: 2 1 1 2 1 1 1 1 1 1 ...
> acf(uvl)
> plot(stl(uvl, s.window="periodic")$time.series)
(Графики мы предлагаем вам построить самостоятельно.)
Видно, что на этот раз присутствует некая периодичность, с интер- валом около 0.2 суток (примерно 5 часов). А вот выраженного тренда нет.
Ответ к задаче про цикл
. Вот как можно его сделать:
> for (m in 1:14)
+ {
+ assign(paste("model0",m,sep=""), arima(cum.polzovateli,
+ order=c(0,0,m)))
+ }
Прогноз
189
Функция assign() заменяет стрелочку (<-) и при этом позволяет программно менять название объекта — адреса присвоения.
Ответ к задаче о курсе доллара
. Поступим ровно так же, как в примере о пользователях провайдера, то есть сначала превратим дан- ные во временной ряд:
> dollar1 <- read.table("data/dollar.txt", dec=",")$V3
> dollar <- ts(dollar1[1:56], frequency=7)
(Десятичным разделителем была запятая! Еще мы учли, что курс доллара имеет недельную периодичность, и организовали данные поне- дельно.)
Затем проверим разные коэффициенты модели ARIMA:
> for (m in 1:7)
+ {
+ mm <- paste("model0", m, sep="")
+ assign(mm, arima(dollar, order=c(0,0,m)))
+ cat(paste(mm, ": ", AIC(get(mm)), "\n", sep=""))
+ }
model01: -67.8129312683548
model02: -80.9252194390681
model03: -82.7498433251648
model04: -82.4022042909309
model05: -84.5913013237983
model06: -83.0836200480987
model07: -82.2056122336345
> for (m in 0:5)
+ {
+ mm <- paste("model", m, "05", sep="")
+ assign(mm, arima(dollar, order=c(m,0,5)))
+ cat(paste(mm, ": ", AIC(get(mm)), "\n", sep=""))
+ }
model005: -84.5913013237983
model105: -82.772885907062
model205: -81.8467316861654
model305: -79.9942752423287
model405: -83.3710277055304
model505: -80.7835758224462
> for (m in 0:5)
+ {
+ mm <- paste("model0", m, "5", sep="")
+ assign(mm, arima(dollar, order=c(0,m,5)))
190
Узнаем будущее: анализ временных рядов
+ cat(paste(mm, ": ", AIC(get(mm)), "\n", sep=""))
+ }
model005: -84.5913013237983
model015: -78.0533071948488
model025: -69.4383206785473
model035: -58.0629434523778
model045: -36.1736715036462
model055: -18.2117126978254
Для ускорения процесса был придуман цикл, который не только меняет значения коэффициентов и имена моделей, но еще и вычисляет и выводит AIC для каждой из них. Обратите внимание на функцию get()
, ее использование в чем-то противоположно assign(): если име- ется имя name, то get(name) будет искать объект с именем name и пере- давать его наружу именно как объект, а не как строку текста. Функция cat()
использовалась для печати «наружу», без нее цикл бы ничего не вывел.
Итого, модель model005 имеет минимальный AIC, и, стало быть,
нам лучше выбрать для предсказания именно ее. Построим график предсказания, добавив туда еще и реальные значения курса доллара за последние две недели, с 27 июля по 9 августа (рис. 60):
> plot(dollar, xlim=c(1,11), ylim=c(27.3,28.5))
> lines(predict(model005, n.ahead=14, se.fit = TRUE)$pred,
+ lty=2)
> lines(ts(dollar1[56:70], start=9, frequency=7))
Итак, наша модель не смогла предсказать резкое падение курса 29
июля, но тем не менее указала на некоторое увеличение курса в следу- ющие две недели. На самом деле результаты даже лучше, чем
можно было ожидать, поскольку график автокорреляций (проверьте это са- мостоятельно) показывает, что значимых автокорреляций в нашем ряду мало.
Прогноз
191 2
4 6
8 10 27
.4 27
.6 27
.8 28
.0 28
.2 28
.4
Рис. 60. Предсказание курса доллара (предсказанные значения обозна- чены пунктирной линией). По оси x отложены недели, прошедшие с 1
июня 2011 года
Глава 8
Статистическая разведка
Вот и окончилось наше изложение основных методов статистики и их использования в R. Теперь настало время применить знания на практике. Но перед тем, как сделать этот шаг, неплохо представить по- лученные знания в «концентрированном» виде. Итак, как нужно ана- лизировать данные?
8.1. Первичная обработка данных
Вначале надо выяснить несколько вопросов:
Скачать 3.04 Mb.