Вычислять среднее значение подмножества вектора только тогда, когда значения подмножества соответствуют условию в R?

Louise

8 августа 2021 в 20:40

У меня есть дневная кривая x, и я пытаюсь аппроксимировать средние пиковые и непиковые значения x: https://ibb.co/Fq1Byzk

Я определил пороговое значение дельты таким образом, что, когда дельта ниже порогового значения, x будет находиться в периоде низкой или максимальной нагрузки. Я хочу получить среднее пиковое значение, где среднее значение относится только к значениям в пределах x, где дельта <порога. Прямо сейчас он также усредняет выбросы.

delta <- matrix(0,24,ncol=1)

for (i in 2:24){
  # i-th element is the i-th hour per day
  delta[i] = x[i,2]-x[i-1,2]
}

# Find hour at which max and min daily values occur
max_threshold = 0.15*max(delta)
min_threshold = 0.15*min(delta)
c <- abs(delta) < max_threshold

t1 <- which(delta>max_threshold)[1]-1 # t1: time index at end of off-peak
t2 <- which.max(delta) + 1 # t2 is time of initial peak
t3 <- which.min(delta)-2 # t3 is time of end peak
t4 <- which.min(delta) # t4 time index of evening off-peak

am <- mean(x[1:t1,2]) # average morning off-peak value
peak <- mean(x[t2:t3,2]) #average peak value
pm <- mean(x[t4:24,2]) # average evening off-peak value

> dput(x)
structure(list(time = structure(c(1451952000, 1451955600, 1451959200, 
1451962800, 1451966400, 1451970000, 1451973600, 1451977200, 1451980800, 
1451984400, 1451988000, 1451991600, 1451995200, 1451998800, 1452002400, 
1452006000, 1452009600, 1452013200, 1452016800, 1452020400, 1452024000, 
1452027600, 1452031200, 1452034800, 1452038400, 1452042000, 1452045600, 
1452049200, 1452052800, 1452056400, 1452060000, 1452063600, 1452067200, 
1452070800, 1452074400, 1452078000, 1452081600, 1452085200, 1452088800, 
1452092400, 1452096000, 1452099600, 1452103200, 1452106800, 1452110400, 
1452114000, 1452117600, 1452121200), class = c("POSIXct", "POSIXt"
), tzone = "UTC"), Crow_education_Omer = c(0.019186330898848, 
0.0192706664192825, 0.0182164724138513, 0.018174304653634, 0.019355001939717, 
0.0197345117816722, 0.023951287803397, 0.0323848398468467, 0.0343245568168401, 
0.0378244809148717, 0.0393003525224754, 0.0403545465279066, 0.0405232175687756, 
0.0393425202826927, 0.0398907011655169, 0.0377401453944372, 0.0344932278577091, 
0.0317101556833707, 0.0304872906370705, 0.0297282709531601, 0.0287584124681633, 
0.0252584883701317, 0.0196080085010205, 0.0197345117816722, 0.0194815052203687, 
0.0196080085010205, 0.0184273112149375, 0.0184694789751548, 0.0191441631386307, 
0.019692344021455, 0.025469327171218, 0.0352522475416196, 0.0376136421137855, 
0.0403967142881239, 0.0435592963044175, 0.0433484575033313, 0.0430532831818105, 
0.042968947661376, 0.043306289743114, 0.044655658070066, 0.0424207667785518, 
0.0416195793344241, 0.0382883262772615, 0.03769797763422, 0.0330173562501054, 
0.0281680638251219, 0.0234452746807901, 0.0225597517162278)), row.names = 97:144, class = "data.frame")

Кроме того, как мне построить новую упрощенную кривую вместе с исходной кривой x на одном графике? Кажется, я не могу расплавить или rbind() новую кривую с уменьшенным количеством точек данных с x, так как мой столбец времени имеет значение POSIXCT.

Спасибо.

Источник

Поделиться Скопировать ссылку Скопировано

Martin Gal

8 августа 2021 в 21:30

Не могли бы вы объяснить, я хочу получить среднее пиковое значение, где среднее значение состоит только из значений в пределах x, где дельта < порог немного дальше? Например: Вы хотите включить в am только значения x[1:t1,2], чьи delta < max_threshold? А как насчет min_threshold?

Louise

8 августа 2021 в 21:37

так, например, я хочу, чтобы мое среднее пиковое значение включало только значения с 10:00 до 14:00, потому что дельта, вычисленная между 10:00 и 14:00, меньше, чем мое значение max_threshold ~0,002. Прямо сейчас мое среднее пиковое значение включает все от t2 до t3. Я вычислил max_threshold и min_threshold, так как они предназначались для определения t2 и t3... но я не был уверен, как это сделать, поэтому я просто установил t2 и t3 на основе добавления 1-2 часов к максимальной и минимальной дельтам. В идеале t2 и t3 должны основываться именно на том, где значения дельты превышают максимальный и минимальный пороги.

Martin Gal

8 августа 2021 в 21:51

Ваша дельта в 13:00 составляет около -0.00118. Это < max_threshold, но также и < min_threshold. Почему этот должен быть включен? (Извините за глупые вопросы, но я еще не понял алгоритм...). В 14:00 дельта равна < max_thr и снова > min_thr.

Louise

8 августа 2021 в 22:14

вы абсолютно правы. я пересмотрел вектор c выше, чтобы вместо этого показать абсолютное значение дельты. В основном я хочу усреднить x между t2 и t3 по индексам, где abs(delta) < max_thr

Ответы (1)

Martin Gal

8 августа 2021 в 23:24

Это всего лишь частичное решение, так как ломается уже второй день. Я назвал data.frame df вместо x.

library(ggplot2)
library(dplyr)
library(lubridate)

df_obj <- df %>% 
  group_by(day = day(time)) %>%  # group by days
  filter(day == 5) %>%           # filter for day 5
  mutate(
    delta_rev = Crow_education_Omer - lag(
      Crow_education_Omer, 
      default = first(Crow_education_Omer)
      ), # delta between day n and n-1
    delta_for = lead(
      Crow_education_Omer, 
      default = last(Crow_education_Omer)
      ) - Crow_education_Omer, # delta between day n-1 and n
    max_tresh  = 0.15 * max(delta_rev)
    )  %>%
  group_by(grp   = 1 - (abs(delta_rev) < 0.15 * max(delta_rev) | abs(delta_for) < 0.15 * max(delta_for)),
           grp2 = cumsum(grp != lag(grp, default = 0))
  ) %>% 
  mutate(
    average = mean(Crow_education_Omer) * 
      (1 - grp) * 
      (abs(first(Crow_education_Omer) - last(Crow_education_Omer)) < max_tresh)
  )

Сначала нам нужно изменить ваш существующий data.frame, чтобы построить ваши средние значения. На основе этого расчета мы используем ggplot2 для построения графика:

df_obj %>% 
  ggplot(aes(x = time, y = Crow_education_Omer)) +
  geom_point() +
  geom_line(aes(color = "sample")) +
  geom_line(data = df_obj[df_obj$average != 0, ], aes(x = time, y = average, color = "average")) +
  xlab("Time") +
  ylab("Value")

возвращает

Но для day 6 это не работает должным образом: изменение на filter(day == 6) и повторное построение графика возвращает

, что не является ожидаемым результатом. Изменение порогового значения на 0.33 * max(delta) и повторное построение графика создает

Итак, возможно, вы сможете использовать этот код, чтобы создать правильное и работающее решение. Удачи!

Несколько пояснений:

Наращиваем delta_rev и delta_for. delta_rev равно вашему delta, поэтому для данной строки/точки данных i мы вычисляем df[i,2] - df[i-1,2].
delta_for меняет это, теперь мы вычисляем df[i + 1,2] - df[i,2] для заданного i. Моя идея здесь такова: использование обоих, delta_rev и delta_for позволяет нам посмотреть на предыдущие и последующие точки. Это дает нам больше информации о соседях данной точки и полезно, чтобы определить, принадлежит ли точка к группе (am,peak,pm).
Функция group_by пытается построить группы на основе порогового значения. grp проверяет, является ли точка данных < 0.15 max(delta), grp2 создает уникальный групповой номер.

Есть несколько проблем:

На основе этого алгоритма может быть более трех групп.
group_by находит другую группу между 15:00 и 20:00, мы отфильтровываем ее (это часть (abs(first(Crow_education_Omer) - last(Crow_education_Omer)) < max_tresh)). Я не уверен, что это хорошее решение.
Как указано выше, это не возвращает разумный график для дня 6. Возможно, это вызвано df_obj[df_obj$average != 0, ]-частью geom_point.

Поделиться Скопировать ссылку Скопировано

Louise

9 августа 2021 в 02:27

Привет, не могли бы вы объяснить, что происходит в коде df_obj? Я не понимаю, что происходит после того, как x фильтруется только до 5-го дня. Что делают delta_for и delta_rev?

Martin Gal

9 августа 2021 в 06:31

Я добавлю объяснение позже.

Martin Gal

9 августа 2021 в 10:08

@Луиза добавила несколько пояснений. Я упустил одну вещь: очевидно, я использовал множество dplyr-функций/синтаксисов. Надеюсь, вы тверды с этим. Особенно легко создать delta с помощью dplyr, и мы избегаем использования цикла for.