Временные ряды, плотность измерения

Question

Я измеряю ежедневную продолжительность (мин) повторяющегося события (E) в течение 364 дней.

ev1<-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 2.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3.27, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 370.33, 1375.4, 
1394.03, 1423.8, 1360, 1269.77, 1378.8, 1350.37, 1425.97, 1423.6, 
1363.4, 1369.87, 1365.5, 1294.97, 1362.27, 1117.67, 1026.97, 
1077.4, 1356.83, 565.23, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 356.83, 
973.5, 0, 240.43, 1232.07, 1440, 1329.67, 1096.87, 1331.37, 1305.03, 
1328.03, 1246.03, 1182.3, 1054.53, 723.03, 1171.53, 1263.17, 
1200.37, 1054.8, 971.4, 936.4, 968.57, 897.93, 1099.87, 876.43, 
1095.47, 1132, 774.4, 1075.13, 982.57, 947.33, 1096.97, 929.83, 
1246.9, 1398.2, 1063.83, 1223.73, 1174.37, 1248.5, 1171.63, 1280.57, 
1183.33, 1016.23, 1082.1, 795.37, 900.83, 1159.2, 992.5, 967.3, 
1440, 804.13, 418.17, 559.57, 563.87, 562.97, 1113.1, 954.87, 
883.8, 1207.1, 1046.83, 995.77, 803.93, 1036.63, 946.9, 887.33, 
727.97, 733.93, 979.2, 1176.8, 1241.3, 1435.6)

ev2<-c(0, 369.3, 158.2, 347.7, 312.5, 265.47, 334.73, 420.83, 816.9, 
925.6, 926.33, 925.4, 917.57, 675.27, 0, 426.03, 860.03, 1041.43, 
947.8, 1076.83, 709.5, 1014.17, 660.3, 428.2, 718.03, 920.8, 
810, 528.53, 103.83, 300.37, 822.03, 662.13, 393.83, 622.47, 
994.13, 1034.07, 893.8, 643.37, 605.07, 360.97, 158.13, 0, 0, 
678.33, 347.67, 384.87, 495.9, 231.37, 443.23, 638.1, 559.53, 
354, 220.13, 210.4, 425.77, 159.5, 260.13, 1132.9, 77.67, 263.83, 
276.23, 63.6, 1.97, 0, 765.2, 403.03, 214.4, 550.63, 752.47, 
58.7, 475.1, 776.4, 53.87, 106.07, 63.23, 425.5, 461.4, 172.73, 
764.8, 53.27, 20.7, 322.8, 228, 36.07, 27.23, 0, 66.3, 389.77, 
705.23, 9.9, 739.3, 883.73, 0, 0, 347.9, 831.43, 0, 28.2, 4.37, 
596.67, 973.7, 26.33, 0.03, 5.93, 777, 918.43, 0, 54.57, 888.13, 
92.83, 98.13, 808.17, 310.5, 263.57, 248.13, 133.37, 138.37, 
14.73, 55.27, 7.17, 242.6, 206.5, 62.97, 8.67, 670.03, 215.77, 
101, 14.07, 440.33, 603.6, 28.27, 257.07, 64.4, 36.4, 506.17, 
333.3, 121.83, 566, 4.33, 192.83, 77.83, 101.3, 261.67, 15.03, 
298.67, 0.3, 616.4, 90.9, 250.87, 323.17, 36.5, 205.2, 205.3, 
110.67, 33.43, 613.43, 95.27, 3.9, 558.7, 650.83, 0, 179.7, 40.6, 
217.13, 48.23, 423.67, 33.9, 176.3, 139.93, 31.63, 0, 162.77, 
311.47, 22.2, 128.3, 0, 304.9, 281.4, 140.73, 131.8, 393.5, 48.63, 
18.17, 232.7, 294.87, 207.6, 317.13, 51.87, 262.57, 70.73, 9.57, 
480.57, 491.37, 27.03, 625.37, 364.4, 0, 79.93, 723.3, 231.57, 
56.93, 836.43, 713.57, 16.8, 2.23, 56.67, 307.87, 466.77, 270.1, 
143.63, 686.23, 703.77, 0, 167.87, 152.6, 237.97, 278.03, 190.7, 
554.03, 37.5, 177.2, 69.2, 119.13, 225.4, 471.23, 7.43, 273.5, 
75.57, 226.73, 141.17, 40.83, 217.33, 238.2, 15.1, 281.27, 244.03, 
0.83, 186.8, 165.53, 142.1, 121.53, 138.83, 103.5, 42.03, 64.27, 
132.07, 26.73, 150.97, 0, 239.9, 100.47, 95.9, 78.23, 90.73, 
172.7, 9.17, 79.77, 67.67, 2.87, 136.73, 362.1, 78.23, 409.37, 
38.9, 62.73, 459.1, 352.6, 17.43, 241.27, 193.1, 278.4, 124.73, 
256.53, 152.6, 247.03, 229.3, 16.5, 73.9, 0, 545.47, 157.5, 182.2, 
276.57, 76.8, 284.43, 2.83, 1.17, 272.57, 314.77, 98.8, 219.93, 
115.23, 121.77, 453.23, 261.73, 101.83, 381, 118.33, 328.23, 
344, 179.5, 16.7, 99.13, 202.97, 57.57, 83.13, 206.87, 425.27, 
130.97, 113.17, 12.07, 207.4, 77.5, 104.7, 59.77, 59.1, 166.6, 
121.2, 139.77, 96.4, 44.23, 262.6, 61.97, 173.2, 281.03, 27.77, 
91.33, 331.23, 142.73, 103.17, 155.7, 80.47, 52.7, 28.6, 56.67, 
257.23, 90.43, 19.43, 69.43, 358.6, 77.9, 15.07, 592.9, 597.27, 
16.83, 225.53, 176.67, 211.47, 159.83, 211, 187.27, 269.73, 27.1, 
421, 83.1, 11.1, 11.67, 253.1, 326.33, 74.33, 153.93, 12.03, 
70.9, 84.47)

Оба индивида (ev1, ev2) имеют приблизительно одинаковую общую продолжительность событий, однако временной «спред» больше в ev2 и более «сфокусирован» в ev1

plot(1:364, ev1, type="l", xlab="Days", ylab="Daily Event duration", main="ev1")
plot(1:364, ev2, type="l", xlab="Days", ylab="Daily Event duration", main="ev2")

Я хотел бы описать или дать количественную оценку этого временного спреда или агрегации продолжительности ежедневного события во времени. Есть ли стандартизированный способ сделать это?

Я думал примерно так: какое минимальное количество дней будет составлять x процентов от общей продолжительности мероприятия. Для приведенного выше примера это минимальное количество дней будет больше для ev2, чем ev1. Есть ли способ рассчитать это?

Любые мысли или ссылки будут полезны.

Большое спасибо

Answer 1

Если вы хотите, чтобы «время» было индексом, то, возможно, было бы проще работать с представлением, которое явно признает, что:

dfrm <- data.frame(tm <-1:364, ev1=ev1, ev2=ev2)

Поскольку вас действительно интересует «плотность»значения индекса ("tm"), используйте весовой аргумент для плотности:

 ev1dens <-  density(dfrm$tm, weights=dfrm$ev1/sum(dfrm$ev1), from=0, to=364, n=364)
 plot( ev1dens, lwd=5)
 which.max(ev1dens$y)
#[1] 326
 abline(v=326)  #

Теперь (поскольку, к счастью, плотности являются унимодальными), необходимо отсортировать нормализованные значения плотности в порядке убывания и найтиИндекс, при котором совокупность становится> целевой пропорцией:

 which(cumsum(ev1dens$y[ order(ev1dens$y, decreasing =TRUE) ])/sum(ev1dens$y) > 0.9)[1]
#[1] 124
 ev1dens$x[order(ev1dens$y, decreasing =TRUE) ][124]
#1] 240.6612

Я приложил усилия, чтобы определить, где будет установлена ​​точка отсечения для включения 90%, но после рассмотрения ответа Томми на ваши последующиена вопрос мои опасения по поводу точности этого метода усиливаются.124 будет индексом для точки отсечения, захватывающей 90%, а 240 будет значением х.Посмотрите на нанесенную на график последовательность нисходящих значений ev1dens $ y, пройденных во время кумулятивного процесса, нанесенного пунктирной красной линией, и зеленую вертикальную линию, где уровень 90% в конечном итоге накоплен:

Вы можете проверить совместное распределение тм и двух векторов.

  require(hexbin)
 hexev1 <- with(dfrm,  hexbin(tm, ev1))
 plot(hexev1)
 hexev2 <- with(dfrm,  hexbin(tm, ev2))
 plot(hexev2)
 plot(hexev1)

Индекс для получения x% от общего числа (который, я думаю, совершенно отличается от кластеризации, приведенной выше, просто:

> min(which(cumsum(ev1) >= sum(ev1)*(x/100) ) )
[1] 317
> min(which(cumsum(ev2) >= sum(ev2)*(x/100) ) )
[1] 112

Answer 2

Звучит как то, что вы хотите, это временное стандартное отклонение.Вы можете вычислить среднее время, взвешенное наблюдением (или медианой, используя код в ответе DWin с x = 50), а затем взять квадратный корень из среднего квадрата разницы, взвешенной наблюдением.(Использование кадра данных и кода DWin)

t.med <- min(which(cumsum(ev1) >= sum(ev1)*(.5)))
t.sd <- with(dfrm, sqrt(mean(((tm - t.med) * ev1)^2)))

Это дает 174031 для ev1 и 41416 для ev2.

В качестве альтернативной статистики, как доля ненулевых измерений, так иДлина самого длинного последовательного прогона с ненулевыми измерениями, поделенными на общую длину, оба очень просты, понятны и, кажется, отражают вашу точку зрения.