Как рассчитать максимальные 8-часовые скользящие (скользящие) средние в R?
Я использую R. Я знаю, что вычисление скользящего среднего - это тема с несколькими ответами на этом сайте, но у меня есть некоторые проблемы, которые делают мой вопрос уникальным.
У меня есть фрейм данных, включающий 8784 почасовых концентрации (366 дней * 24 часа) загрязнителя воздуха (озон). Этот фрейм данных содержит некоторые значения NaN (пропущенные значения). Процедура содержит следующие шаги:
1- вычисление скользящих (скользящих) средних значений за 8 часов для часовых концентраций: т.е. каждые 8 концентраций должны быть усреднены следующим образом: среднее значение от 1 до 8, среднее от 2 до 9, среднее от 3 до 10 и т. Д. Это приводит к получению 24 скользящих средних за каждый день (каждые 24 часа).
2 - для каждого дня я хочу максимум 8-часовых скользящих средних: то есть среди 24 скользящих средних должно быть выбрано наибольшее число. Наконец, будет выбрано 366 скользящих средних (366 дней).
3- Необходимо создать новый фрейм данных, содержащий 366 скользящих средних.
Я знаю, что есть некоторые пакеты (openair, zoo, TTR), которые делают что-то подобное, но есть ли способы написать коды без пакетов?
An Exmaple of my data
ColName
1 18.76
2 12.92
3 8.12
4 NaN
5 12.92
6 3.77
7 18.76
8 9.52
9 94.09
10 18.76
11 14.13
12 8.12
13 2.04
14 12.92
15 9.17
.
.
.
8783 34.58
8784 64.23
Название основного фрейма данных - "Почасово". Я попробовал эти коды:
Hourly1 <- c(0, cumsum(ifelse(is.nan(Hourly), 0, Hourly)))
rsum <- (Hourly1[(Hourly1+1):length(Hourly1)] - Hourly1[1:(length(Hourly1) - 8)]) / 8
Но когда я пробую первую строку, возникает следующая ошибка:
Error in is.nan(Hourly) : default method not implemented for type 'list'
ОБНОВЛЕНИЕ: я использовал следующие коды, но максимум 8-часовых средних не рассчитан правильно:
Hourly2<-as.numeric(Hourly$Average)
names(Hourly2) <- rep(seq.Date(as.Date("2017-01-01"), by=1, length.out=366), each=24)
x<-Hourly2
#use cumsum to get the moving average, keep NaNs
cx <- c(0, cumsum(ifelse(is.nan(x), 0, x))) + c(0,x)*0
n <- 8
rsum <- (cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]) / n
res <- data.frame(mov_avg=rsum, days=names(rsum))
#select max from each day, ignoring NaN's
mx <- aggregate(mov_avg~days, data=res, max)
Я сравнил окончательные результаты (максимум 366 из 8-часовых средних, каждый за 1 день года) со стандартным предварительно утвержденным набором данных. В некоторые дни коды рассчитывают средние значения правильно, а в другие дни нет! Я не понял его логику.
Вы можете найти мой набор необработанных данных здесь!
ОБНОВЛЕНИЕ 2:
Вот ссылка для скачивания окончательных результатов, полученных разными методами!
Update3:
Разница между результатами была обусловлена различными методами расчета скользящих средних. Существует три метода расчета скользящих средних: слева, справа и по центру. Коды, предложенные ребятами здесь, следуют "правильному" методу.
2 ответа
Вот пример того, как это сделать с cumsum
когда у вас есть пропущенные значения. Я бы внимательно рассмотрел, как они распределяются в ваших данных и как вы хотите с ними справляться.
#create some sample data
set.seed(1)
x <- rnorm(24*366)
names(x) <- rep(seq.Date(as.Date("2017-01-01"), by=1, length.out=366), each=24)
x[sample(100, 1:length(x))] <- NaN #add some missing values
#use cumsum to get the moving average, keep NaNs
cx <- c(0, cumsum(ifelse(is.nan(x), 0, x))) + c(0,x)*0
n <- 8
rsum <- (cx[(n+1):length(cx)] - cx[1:(length(cx) - n)]) / n
res <- data.frame(mov_avg=rsum, days=names(rsum))
#select max from each day, ignoring NaN's
mx <- aggregate(mov_avg~days, data=res, max)
days mov_avg
1 2017-01-01 0.6404849
2 2017-01-02 0.3456389
3 2017-01-03 0.5998888
4 2017-01-04 0.6635502
5 2017-01-05 0.7244289
6 2017-01-06 0.1715349
Я работал именно над этим, и нашел решение, которое использует map2()
# create a day of ozone data
o3day <- data.frame(o3hrly = runif(24, 0.04, 0.1))
# 8hr average function
avg_8hr <- function(.x, .y, o3) {
# print(.x)
# print(.y)
# print(o3)
o3 %>% slice(.x:.y) %>% summarize(o38hr = mean(o3hrly))
}
max(unlist(map2(.x = 1:17, .y = 8:24, .f = avg_8hr, o3 = o3day)))