Как рассчитать максимальные 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)))
Другие вопросы по тегам