Если нас поразит столкновение, приведшее к вымиранию человечества, оно, скорее всего, произойдет весной или осенью.
Введение
Если бы вы спросили 100 человек, что, по их мнению, является самым большим риском для человеческой цивилизации, я бы поспорил, что тремя главными ответами были бы ядерная война, глобальная пандемия и глобальное потепление/изменение климата. Однако менее 10 лет назад метеор диаметром около 20 метров и массой 10 000 тонн взорвался в 30 км над городом Челябинск в России. Хотя погибших не было, по оценкам, в результате взрыва был нанесен ущерб на сумму 30 миллионов долларов и ранено 1500 человек. Около 100 лет назад, в 1908 году, над Сибирью взорвался метеор размером 50–60 метров мощностью 12 мегатонн, уничтоживший около 2200 квадратных километров леса. 66 миллионов лет назад астероид диаметром 12 км врезался в Землю со скоростью 43 000 км/ч, оставив кратер диаметром 200 км и убив по меньшей мере 80% всего живого на Земле. Возможно, это не так очевидно, как ядерная война или глобальная пандемия, учитывая текущую ситуацию на Земле, но я считаю, что разрушительное столкновение с цивилизацией продолжает представлять серьезную угрозу для выживания человечества в долгосрочной перспективе.
В этом анализе я буду использовать модели машинного обучения (ML) ETS, ARIMA и LSTM, чтобы попытаться понять тенденции сезонности и предсказать близкое сближение с Землей околоземных объектов (NEO).
Данные
Данные об ОСЗ были получены из Центра изучения околоземных объектов Лаборатории реактивного движения НАСА: https://cneos.jpl.nasa.gov/ca/. Данные были собраны с 1 января 2010 г. по 31 декабря 2021 г. и учитывали только ОСЗ с расстоянием между их центром и центром Земли ≤ 0,05 астрономических единиц (7 479 894 км).
Подготовка данных и исследование
# Load Libraries #
library(forecast)
library(dplyr)
library(quantmod)
library(tidyverse)
library(glue)
library(forcats)
library(timetk)
library(tidyquant)
library(tibbletime)
library(cowplot)
library(recipes)
library(rsample)
library(yardstick)
library(keras)
# Read and prep Data #
data <- read.csv("...NEO.csv")
tail(data)
dim(data)
# Convert Date column to date type # df <- mutate(data, Date= as.Date(Date, format= "%Y-%m-%d")) # Visualize # ggplot(data=df, aes(x=Date, y=Count, group=1)) + geom_line()+ geom_point()
Мы установим набор данных для обучения с 1 января 2010 г. по 31 декабря 2020 г. и набор данных для тестирования с 1 января 2021 г. по 31 декабря 2021 г.:
# Set up Training and Test Sets # train <- subset(df, Date < "2021-01-01") # 2010 to 2020 data train_count <- train$Count train.ts <- ts(train_count, frequency=12, start=c(2010,1)) test <- subset(df, Date >= "2021-01-01") # 2021 data test_count <- test$Count test.ts <- ts(test_count, frequency=12, start=c(2021,1)) seasonplot(train.ts)
СТВ
Модели ETS обозначают:
- Ошибкаошибка: может быть Aаддитивной или Ммультипликативной
- T rend: может быть Nодин, Aдобавочный или Mмногократный
- S easonality: может быть None, Additiv или Mmultiplicative.
Так, например, модель с «AAA» означает модель с аддитивной ошибкой, аддитивной тенденцией и аддитивной сезонностью. В этом примере я установлю модель «ZZZ», которая выберет для меня оптимальную модель.
# ETS # ets_fit <- ets(train.ts, model="ZZZ") autoplot(ets_fit)
Модель ETS для данных NEO — «MAM», что означает наличие Mмножественной Eошибки, Aддитивной T ренд и Mмногократная Sсезонность. Насколько хорошо эта модель работала на тестовых данных?:
ets_future <- as.data.frame(forecast(ets_fit, h=12)) df_ets <- cbind(test, ets_future) MAPE <- mean(abs(df_ets$Count - df_ets$`Point Forecast`)/df_ets$Count) MAE <- mean(abs(df_ets$Count - df_ets$`Point Forecast`)) MSE<-mean((df_ets$Count - df_ets$`Point Forecast`)^2) RMSE<-sqrt(MSE)
Модель ETS была чрезвычайно точной при прогнозировании ОСЗ на данных испытаний, отличаясь только средним абсолютным средним значением 17 (средняя абсолютная ошибка в процентах 12%).
АРИМА
ARIMA означает авторегрессивные интегрированные скользящие средние. Подобно моделям LSTM, которые мы обсудим ниже, модели ARIMA работают на основе автокорреляций, обнаруженных в данных. Другими словами, они используют корреляцию между значениями, которые предшествуют и следуют за ними, чтобы прогнозировать будущее.
# ARIMA # arima_fit <- auto.arima(train.ts) summary(arima_fit)
Модель ARIMA, сгенерированная из данных NEO, представляет собой ARIMA(0, 1, 1)(0, 1, 1), которая является наиболее часто используемой сезонной моделью ARIMA. Я не буду вдаваться во все подробности, но хотел показать уравнение:
Ŷₜ = Yₜ-₁₂ + Yₜ-₁ — Yₜ-₁₃ — (-0.8794)eₜ-₁ — (-0.5597)eₜ-₁₂ + ((-0.8794)(-0.5597)eₜ-₁₃)
На английском языке это уравнение говорит, что прогнозируемое количество ОСЗ будет зависеть от количества ОСЗ за 12 месяцев до этого, плюс число в предыдущем месяце минус число 13 месяцев назад, а затем добавляется кратность ошибок при лагах 1, 12 и 13. Стоит отметить, что коэффициент ошибки lag-13 является произведением коэффициентов MA(1) и SMA(1), показанных в сводке выше. Теперь, когда мы понимаем наше уравнение, насколько хорошо работает эта модель?:
arima_future <- as.data.frame(forecast(arima_fit, h=12)) df_arima <- cbind(test, arima_future) MAPE <- mean(abs(df_ets$Count - df_arima$`Point Forecast`)/df_ets$Count) MAE <- mean(abs(df_ets$Count - df_arima$`Point Forecast`)) MSE<-mean((df_ets$Count - df_arima$`Point Forecast`)^2) RMSE<-sqrt(MSE)
Хотя модель ARIMA была очень точной, отставая от среднего абсолютного среднего только на 20 (средняя абсолютная ошибка в процентах — 18%), модель ETS по-прежнему работала лучше.
ЛСТМ
Я не буду слишком углубляться в LSTM, но я написал статью, в которой содержится более подробная информация, которую вы можете прочитать здесь: https://medium.datadriveninvestor.com/utilizing-long-short-term-memory- сети, чтобы спрогнозировать цену биткойна в r-8a77a6512ddc
Наш первый шаг с LSTM — проверка автокорреляции:
# LSTM #
tidy_acf <- function(df, Count, lags = 0:20) {
value_expr <- enquo(Count)
acf_values <- df %>%
pull(Count) %>%
acf(lag.max = tail(lags, 1), plot = FALSE) %>%
.$acf %>%
.[,,1]
ret <- tibble(acf = acf_values) %>%
rowid_to_column(var = "lag") %>%
mutate(lag = lag - 1) %>%
filter(lag %in% lags)
na.action = na.omit
return(ret)
}
max_lag <- 150
# Find NEO's ACF Values #
df %>%
tidy_acf(Count, lags = 0:max_lag)
df %>%
tidy_acf(Count, lags = 0:max_lag) %>%
ggplot(aes(lag, acf)) +
geom_segment(aes(xend = lag, yend = 0), color = palette_light()[[1]]) +
geom_vline(xintercept = 120, size = 3, color = palette_light()[[2]]) +
annotate("text", label = "10 Year Mark", x = 75, y = 0.8,
color = palette_light()[[2]], size = 6, hjust = 0) +
theme_tq() +
labs(title = "Near Earth Objects")
# Train and Test Data #
periods_train <- 132
periods_test <- 12
rolling_origin_resamples <- rolling_origin(
df,
initial = periods_train,
assess = periods_test,
cumulative = FALSE,
skip = FALSE
)
rolling_origin_resamples
# Create a function that plots the Split #
plot_split <- function(split, expand_y_axis = TRUE, alpha = 1, size = 1, base_size = 14) {
# Manipulate data
train_tbl <- training(split) %>%
add_column(key = "training")
test_tbl <- testing(split) %>%
add_column(key = "testing")
data_manipulated <- bind_rows(train_tbl, test_tbl) %>%
as_tbl_time(index = Date) %>%
mutate(key = fct_relevel(key, "training", "testing"))
# Collect attributes
train_time_summary <- train_tbl %>%
tk_index() %>%
tk_get_timeseries_summary()
test_time_summary <- test_tbl %>%
tk_index() %>%
tk_get_timeseries_summary()
# Visualize
g <- data_manipulated %>%
ggplot(aes(x = Date, y = Count, color = key)) +
geom_line(size = size, alpha = alpha) +
theme_tq(base_size = base_size) +
scale_color_tq() +
labs(
title = glue("Split: {split$id}"),
subtitle = glue("{train_time_summary$start} to {test_time_summary$end}"),
y = "", x = ""
) +
theme(legend.position = "none")
if (expand_y_axis) {
NEO_time_summary <- df %>%
tk_index() %>%
tk_get_timeseries_summary()
g <- g +
scale_x_date(limits = c(NEO_time_summary$start,
NEO_time_summary$end))
}
return(g)
}
rolling_origin_resamples$splits[[1]] %>%
plot_split(expand_y_axis = TRUE) +
theme(legend.position = "bottom")
split <- rolling_origin_resamples$splits[[1]]
split_id <- rolling_origin_resamples$id[[1]]
df_trn <- training(split)
df_tst <- testing(split)
df <- bind_rows(
df_trn %>% add_column(key = "training"),
df_tst %>% add_column(key = "testing")
) %>%
as_tbl_time(index = Date)
df
# Standardize Data #
rec_obj <- recipe(Count ~ ., df) %>%
step_sqrt(Count) %>%
step_center(Count) %>%
step_scale(Count) %>%
prep()
df_processed_tbl <- bake(rec_obj, df)
df_processed_tbl
center_history <- rec_obj$steps[[2]]$means["Count"]
scale_history <- rec_obj$steps[[3]]$sds["Count"]
c("center" = center_history, "scale" = scale_history)
nrow(df_tst)
dim(df_trn)
Теперь, когда у нас настроены и стандартизированы данные для обучения и тестирования, мы можем приступить к созданию нашего LSTM. Сначала мы должны определить оптимальное количество эпох и оптимальное количество используемых узлов. Мы будем использовать двухслойный LSTM, так как этого должно быть достаточно для нашей проблемы регрессии.
# Model inputs
lag_setting <- 12 # = nrow(df_tst)
batch_size <- 12
train_length <- 132
tsteps <- 1
epochs <- 50
# Training Set
lag_train_tbl <- df_processed_tbl %>%
mutate(Count_lag = lag(Count, n = lag_setting)) %>%
filter(!is.na(Count_lag)) %>%
filter(key == "training") %>%
tail(train_length)
x_train_vec <- lag_train_tbl$Count_lag
x_train_arr <- array(data = x_train_vec, dim = c(length(x_train_vec), 1, 1))
y_train_vec <- lag_train_tbl$Count
y_train_arr <- array(data = y_train_vec, dim = c(length(y_train_vec), 1))
# Testing Set
lag_test_tbl <- df_processed_tbl %>%
mutate(
Count_lag = lag(Count, n = lag_setting)
) %>%
filter(!is.na(Count_lag)) %>%
filter(key == "testing")
x_test_vec <- lag_test_tbl$Count_lag
x_test_arr <- array(data = x_test_vec, dim = c(length(x_test_vec), 1, 1))
y_test_vec <- lag_test_tbl$Count
y_test_arr <- array(data = y_test_vec, dim = c(length(y_test_vec), 1))
model <- keras_model_sequential()
model %>%
layer_lstm(units = 88,
input_shape = c(tsteps, 1),
batch_size = batch_size,
return_sequences = TRUE,
stateful = TRUE) %>%
layer_lstm(units = 88,
return_sequences = FALSE,
stateful = TRUE) %>%
layer_dense(units = 1)
model %>%
compile(loss = 'mae', optimizer = 'adam')
model
for (i in 1:epochs) {
model %>% fit(x = x_train_arr,
y = y_train_arr,
batch_size = batch_size,
epochs = 1,
verbose = 1,
shuffle = FALSE)
model %>% reset_states()
cat("Epoch: ", i)
}
# Make Predictions
pred_out <- model %>%
predict(x_test_arr, batch_size = batch_size) %>%
.[,1]
# Retransform values
pred_tbl <- tibble(
Date = lag_test_tbl$Date,
Count = (pred_out * scale_history + center_history)^2
)
# Combine actual data with predictions
tbl_1 <- df_trn %>%
add_column(key = "actual")
tbl_2 <- df_tst %>%
add_column(key = "actual")
tbl_3 <- pred_tbl %>%
add_column(key = "predict")
# Create time_bind_rows() to solve dplyr issue
time_bind_rows <- function(data_1, data_2, Date) {
index_expr <- enquo(Date)
bind_rows(data_1, data_2) %>%
as_tbl_time(index = !! index_expr)
}
ret <- list(tbl_1, tbl_2, tbl_3) %>%
reduce(time_bind_rows, Date = Date) %>%
arrange(key, Date) %>%
mutate(key = as_factor(key))
ret
Очевидно, что для создания LSTM требуется много работы, но обеспечивают ли они гораздо более точные прогнозы?
# Determining Model Performance # MSE<-mean((tbl_2$Count - tbl_3$Count)^2) RMSE<-sqrt(MSE) MAPE <- mean(abs(tbl_2$Count - tbl_3$Count)/tbl_2$Count) MAE <- mean(abs(tbl_2$Count - tbl_3$Count))
Несмотря на то, что модель LSTM работала довольно хорошо, модель ETS по-прежнему была наиболее точной в предсказании ОСЗ 2021 года.
Предсказание будущего
Поскольку модель ETS была наиболее точной, мы будем использовать ее для прогнозирования ОСЗ в 2022 году:
full_count <- df$Count full.ts <- ts(full_count, frequency=12) ets_full_fit <- ets(full.ts, model="ZZZ") ets_full_future <- as.data.frame(forecast(ets_full_fit, h=12))
Прогноз снова предсказывает большое количество ОСЗ весной и осенью 2022 года, причем наибольшее количество, 203, ожидается в октябре, затем 186 в ноябре и 174 в марте. Модель ETS прогнозирует в общей сложности 1666 ОСЗ в 2022 году, что на 72 больше, чем в 2021 году, и на 222 больше, чем в 2020 году.
Выводы
В этом эксперименте я показал, что методы прогнозирования временных рядов могут использоваться для точного прогнозирования близких сближений ОСЗ с Землей. С помощью автокорреляции я показал, что ОСЗ происходят сезонно на основе 12-месячной модели, которая, конечно же, совпадает с орбитой Земли вокруг Солнца. С помощью графиков исторических данных и прогнозов я продемонстрировал, что ожидается увеличение числа ОСЗ, что является тенденцией, которая наблюдается с 1980-х годов, когда исследователи стали лучше осознавать опасности, которые событие столкновения представляет для людей. Я также показал с помощью графиков исторических данных, графиков сезонности и с помощью прогнозов, что большинство ОСЗ происходят весной и осенью, что указывает на то, что если столкновение действительно произойдет, оно, скорее всего, произойдет в октябре/ноябре, а затем в марте/апреле. Однако я не смог найти никакой информации о том, почему существует тенденция увеличения количества NEO весной и осенью. Наконец, я показал, что глубокое обучение, в данном случае моделирование LSTM, не всегда может давать наилучший прогноз и не всегда может стоить дополнительной работы по реализации.
Запланируйте сеанс DDIChat в Data Science / AI / ML / DL:
Подайте заявку на участие в программе DDIChat Expert здесь.
Работайте с DDI: https://datadriveninvestor.com/collaborate
Подпишитесь на DDIntel здесь.
