Фурье / тригонометрическая интерполяция

Фон

В статье Эпштейна (1991): При получении суточных климатологических значений из среднемесячных значений приводятся формулировка и алгоритм расчета интерполяции Фурье для периодических и равномерно распределенных значений.

В статье цель состоит в том, чтобы получить ежедневные значения из месячных средних значений путем интерполяции.

Короче говоря, предполагается, что неизвестные дневные значения могут быть представлены суммой гармонических компонентов: В статье (время) выражается в месяцах.

Y (T) знак равно a_{0} + \underset{J}{Σ} [a_{J} соз (2 π J T / 12) + б_{J} грех (2 π J T / 12)]

$y(t) = a_{0} + \sum_{j}\left[a_{j}\,\cos(2\pi jt/12)+b_{j}\,\sin(2\pi jt/12)\right]$

t

$t$

После некоторого отклонения показано, что члены можно вычислить следующим образом:

\begin{aligned} a_{0} & знак равно \underset{T}{Σ} Y_{T} / 12 \\ a_{J} & знак равно [(π J / 12) / грех (π J / 12)] \times \underset{T}{Σ} [Y_{T} соз (2 π J T / 12) / 6] J знак равно 1, ..., 5 \\ б_{J} & знак равно [(π J / 12) / грех (π J / 12)] \times \underset{T}{Σ} [Y_{T} грех (2 π J T / 12) / 6] J знак равно 1, ..., 5 \\ a_{6} & знак равно [(π J / 12) / грех (π J / 12)] \times \underset{T}{Σ} [Y_{T} соз (π T) / 12] \\ б_{6} & знак равно 0 \end{aligned}

$\begin{align} a_{0} &= \sum_{T}Y_{T}/12 \\ a_{j} &= \left[ (\pi j/12)/\sin(\pi j/12)\right] \times \sum_{T}\left[Y_{T}\,\cos(2\pi jT/12)/6 \right]~~~~~~~j=1,\ldots, 5 \\ b_{j} &= \left[ (\pi j/12)/\sin(\pi j/12)\right] \times \sum_{T}\left[Y_{T}\,\sin(2\pi jT/12)/6 \right]~~~~~~~j=1,\ldots, 5 \\ a_{6} &= \left[ (\pi j/12)/\sin(\pi j/12)\right]\times \sum_{T}\left[Y_{T}\cos(\pi T)/12\right] \\ b_{6} &= 0 \end{align}$ где

Y_{T}

$Y_{T}$ обозначает среднее значение за месяц, а

T

$T$ - месяц.

Harzallah (1995) резюмирует этот подход следующим образом: «Интерполяция выполняется путем добавления нулей к спектральным коэффициентам данных и путем выполнения обратного преобразования Фурье к результирующим расширенным коэффициентам. Этот метод эквивалентен применению прямоугольного фильтра к коэффициентам Фурье «.

Вопросов

Моя цель - использовать вышеуказанную методологию для интерполяции еженедельных средних значений для получения ежедневных данных (см. Мой предыдущий вопрос ). Таким образом, у меня есть 835 еженедельных данных подсчета (см. Пример набора данных внизу вопроса). Есть довольно много вещей, которые я не понимаю, прежде чем я смогу применить подход, изложенный выше:

Как изменить формулы для моей ситуации (еженедельно вместо месячных значений)?
Как можно выразить время ? Я предположил, что (или с точками данных в целом), это правильно? $t$ $t/835$ $t/n$ $n$
Почему автор вычисляет 7 членов (т. )? Сколько терминов я должен рассмотреть? $0\leq j \leq 6$
Я понимаю, что вопрос, вероятно, можно решить с помощью регрессионного подхода и использования прогнозов для интерполяции (спасибо Нику). Тем не менее, некоторые вещи мне неясны: сколько терминов гармоник следует включить в регрессию? И какой период я должен взять? Как можно сделать регрессию, чтобы сохранить средние значения за неделю (поскольку я не хочу точного гармонического соответствия данных)?

Используя регрессионный подход (который также объясняется в этой статье ), мне удалось получить точную гармоническую подгонку к данным ( в моем примере будет проходить через , поэтому я установил 417 терминов). Как можно изменить этот подход - если возможно - чтобы добиться сохранения еженедельных средств? Может быть, применяя поправочные коэффициенты к каждому члену регрессии? $j$ $1, \ldots, 417$ $~$ $~$

График точной гармонической подгонки:

Точная гармоническая посадка

РЕДАКТИРОВАТЬ

Используя пакет сигналов и interp1функцию, вот что мне удалось сделать, используя приведенный ниже пример данных (большое спасибо @noumenal). Я использую, q=7поскольку у нас есть недельные данные:

# Set up the time scale

daily.ts <- seq(from=as.Date("1995-01-01"), to=as.Date("2010-12-31"), by="day")

# Set up data frame 

ts.frame <- data.frame(daily.ts=daily.ts, wdayno=as.POSIXlt(daily.ts)$wday,
                       yearday = 1:5844,
                       no.influ.cases=NA)

# Add the data from the example dataset called "my.dat"

ts.frame$no.influ.cases[ts.frame$wdayno==3] <- my.dat$case

# Interpolation

case.interp1 <- interp1(x=ts.frame$yearday[!is.na(ts.frame$no.influ.case)],y=(ts.frame$no.influ.cases[!is.na(ts.frame$no.influ.case)]),xi=ts.frame$yearday, method = c("cubic"))

# Plot subset for better interpretation
par(bg="white", cex=1.2, las=1)
plot((ts.frame$no.influ.cases)~ts.frame$yearday, pch=20,
     col=grey(0.4),
     cex=1, las=1,xlim=c(0,400), xlab="Day", ylab="Influenza cases")
lines(case.interp1, col="steelblue", lwd=1)

Cubicinterpo

Здесь есть две проблемы:

Кривая кажется «слишком хорошей»: она проходит через каждую точку
Еженедельные средства не сохраняются

Пример набора данных

structure(list(date = structure(c(9134, 9141, 9148, 9155, 9162, 
9169, 9176, 9183, 9190, 9197, 9204, 9211, 9218, 9225, 9232, 9239, 
9246, 9253, 9260, 9267, 9274, 9281, 9288, 9295, 9302, 9309, 9316, 
9323, 9330, 9337, 9344, 9351, 9358, 9365, 9372, 9379, 9386, 9393, 
9400, 9407, 9414, 9421, 9428, 9435, 9442, 9449, 9456, 9463, 9470, 
9477, 9484, 9491, 9498, 9505, 9512, 9519, 9526, 9533, 9540, 9547, 
9554, 9561, 9568, 9575, 9582, 9589, 9596, 9603, 9610, 9617, 9624, 
9631, 9638, 9645, 9652, 9659, 9666, 9673, 9680, 9687, 9694, 9701, 
9708, 9715, 9722, 9729, 9736, 9743, 9750, 9757, 9764, 9771, 9778, 
9785, 9792, 9799, 9806, 9813, 9820, 9827, 9834, 9841, 9848, 9855, 
9862, 9869, 9876, 9883, 9890, 9897, 9904, 9911, 9918, 9925, 9932, 
9939, 9946, 9953, 9960, 9967, 9974, 9981, 9988, 9995, 10002, 
10009, 10016, 10023, 10030, 10037, 10044, 10051, 10058, 10065, 
10072, 10079, 10086, 10093, 10100, 10107, 10114, 10121, 10128, 
10135, 10142, 10149, 10156, 10163, 10170, 10177, 10184, 10191, 
10198, 10205, 10212, 10219, 10226, 10233, 10240, 10247, 10254, 
10261, 10268, 10275, 10282, 10289, 10296, 10303, 10310, 10317, 
10324, 10331, 10338, 10345, 10352, 10359, 10366, 10373, 10380, 
10387, 10394, 10401, 10408, 10415, 10422, 10429, 10436, 10443, 
10450, 10457, 10464, 10471, 10478, 10485, 10492, 10499, 10506, 
10513, 10520, 10527, 10534, 10541, 10548, 10555, 10562, 10569, 
10576, 10583, 10590, 10597, 10604, 10611, 10618, 10625, 10632, 
10639, 10646, 10653, 10660, 10667, 10674, 10681, 10688, 10695, 
10702, 10709, 10716, 10723, 10730, 10737, 10744, 10751, 10758, 
10765, 10772, 10779, 10786, 10793, 10800, 10807, 10814, 10821, 
10828, 10835, 10842, 10849, 10856, 10863, 10870, 10877, 10884, 
10891, 10898, 10905, 10912, 10919, 10926, 10933, 10940, 10947, 
10954, 10961, 10968, 10975, 10982, 10989, 10996, 11003, 11010, 
11017, 11024, 11031, 11038, 11045, 11052, 11059, 11066, 11073, 
11080, 11087, 11094, 11101, 11108, 11115, 11122, 11129, 11136, 
11143, 11150, 11157, 11164, 11171, 11178, 11185, 11192, 11199, 
11206, 11213, 11220, 11227, 11234, 11241, 11248, 11255, 11262, 
11269, 11276, 11283, 11290, 11297, 11304, 11311, 11318, 11325, 
11332, 11339, 11346, 11353, 11360, 11367, 11374, 11381, 11388, 
11395, 11402, 11409, 11416, 11423, 11430, 11437, 11444, 11451, 
11458, 11465, 11472, 11479, 11486, 11493, 11500, 11507, 11514, 
11521, 11528, 11535, 11542, 11549, 11556, 11563, 11570, 11577, 
11584, 11591, 11598, 11605, 11612, 11619, 11626, 11633, 11640, 
11647, 11654, 11661, 11668, 11675, 11682, 11689, 11696, 11703, 
11710, 11717, 11724, 11731, 11738, 11745, 11752, 11759, 11766, 
11773, 11780, 11787, 11794, 11801, 11808, 11815, 11822, 11829, 
11836, 11843, 11850, 11857, 11864, 11871, 11878, 11885, 11892, 
11899, 11906, 11913, 11920, 11927, 11934, 11941, 11948, 11955, 
11962, 11969, 11976, 11983, 11990, 11997, 12004, 12011, 12018, 
12025, 12032, 12039, 12046, 12053, 12060, 12067, 12074, 12081, 
12088, 12095, 12102, 12109, 12116, 12123, 12130, 12137, 12144, 
12151, 12158, 12165, 12172, 12179, 12186, 12193, 12200, 12207, 
12214, 12221, 12228, 12235, 12242, 12249, 12256, 12263, 12270, 
12277, 12284, 12291, 12298, 12305, 12312, 12319, 12326, 12333, 
12340, 12347, 12354, 12361, 12368, 12375, 12382, 12389, 12396, 
12403, 12410, 12417, 12424, 12431, 12438, 12445, 12452, 12459, 
12466, 12473, 12480, 12487, 12494, 12501, 12508, 12515, 12522, 
12529, 12536, 12543, 12550, 12557, 12564, 12571, 12578, 12585, 
12592, 12599, 12606, 12613, 12620, 12627, 12634, 12641, 12648, 
12655, 12662, 12669, 12676, 12683, 12690, 12697, 12704, 12711, 
12718, 12725, 12732, 12739, 12746, 12753, 12760, 12767, 12774, 
12781, 12788, 12795, 12802, 12809, 12816, 12823, 12830, 12837, 
12844, 12851, 12858, 12865, 12872, 12879, 12886, 12893, 12900, 
12907, 12914, 12921, 12928, 12935, 12942, 12949, 12956, 12963, 
12970, 12977, 12984, 12991, 12998, 13005, 13012, 13019, 13026, 
13033, 13040, 13047, 13054, 13061, 13068, 13075, 13082, 13089, 
13096, 13103, 13110, 13117, 13124, 13131, 13138, 13145, 13152, 
13159, 13166, 13173, 13180, 13187, 13194, 13201, 13208, 13215, 
13222, 13229, 13236, 13243, 13250, 13257, 13264, 13271, 13278, 
13285, 13292, 13299, 13306, 13313, 13320, 13327, 13334, 13341, 
13348, 13355, 13362, 13369, 13376, 13383, 13390, 13397, 13404, 
13411, 13418, 13425, 13432, 13439, 13446, 13453, 13460, 13467, 
13474, 13481, 13488, 13495, 13502, 13509, 13516, 13523, 13530, 
13537, 13544, 13551, 13558, 13565, 13572, 13579, 13586, 13593, 
13600, 13607, 13614, 13621, 13628, 13635, 13642, 13649, 13656, 
13663, 13670, 13677, 13684, 13691, 13698, 13705, 13712, 13719, 
13726, 13733, 13740, 13747, 13754, 13761, 13768, 13775, 13782, 
13789, 13796, 13803, 13810, 13817, 13824, 13831, 13838, 13845, 
13852, 13859, 13866, 13873, 13880, 13887, 13894, 13901, 13908, 
13915, 13922, 13929, 13936, 13943, 13950, 13957, 13964, 13971, 
13978, 13985, 13992, 13999, 14006, 14013, 14020, 14027, 14034, 
14041, 14048, 14055, 14062, 14069, 14076, 14083, 14090, 14097, 
14104, 14111, 14118, 14125, 14132, 14139, 14146, 14153, 14160, 
14167, 14174, 14181, 14188, 14195, 14202, 14209, 14216, 14223, 
14230, 14237, 14244, 14251, 14258, 14265, 14272, 14279, 14286, 
14293, 14300, 14307, 14314, 14321, 14328, 14335, 14342, 14349, 
14356, 14363, 14370, 14377, 14384, 14391, 14398, 14405, 14412, 
14419, 14426, 14433, 14440, 14447, 14454, 14461, 14468, 14475, 
14482, 14489, 14496, 14503, 14510, 14517, 14524, 14531, 14538, 
14545, 14552, 14559, 14566, 14573, 14580, 14587, 14594, 14601, 
14608, 14615, 14622, 14629, 14636, 14643, 14650, 14657, 14664, 
14671, 14678, 14685, 14692, 14699, 14706, 14713, 14720, 14727, 
14734, 14741, 14748, 14755, 14762, 14769, 14776, 14783, 14790, 
14797, 14804, 14811, 14818, 14825, 14832, 14839, 14846, 14853, 
14860, 14867, 14874, 14881, 14888, 14895, 14902, 14909, 14916, 
14923, 14930, 14937, 14944, 14951, 14958, 14965, 14972), class = "Date"), 
    cases = c(168L, 199L, 214L, 230L, 267L, 373L, 387L, 443L, 
    579L, 821L, 1229L, 1014L, 831L, 648L, 257L, 203L, 137L, 78L, 
    82L, 69L, 45L, 51L, 45L, 63L, 55L, 54L, 52L, 27L, 24L, 12L, 
    10L, 22L, 42L, 32L, 52L, 82L, 95L, 91L, 104L, 143L, 114L, 
    100L, 83L, 113L, 145L, 175L, 222L, 258L, 384L, 755L, 976L, 
    879L, 846L, 1004L, 801L, 799L, 680L, 530L, 410L, 302L, 288L, 
    234L, 269L, 245L, 240L, 176L, 188L, 128L, 96L, 59L, 63L, 
    44L, 52L, 39L, 50L, 36L, 40L, 48L, 32L, 39L, 28L, 29L, 16L, 
    20L, 25L, 25L, 48L, 57L, 76L, 117L, 107L, 91L, 90L, 83L, 
    76L, 86L, 104L, 101L, 116L, 120L, 185L, 290L, 537L, 485L, 
    561L, 1142L, 1213L, 1235L, 1085L, 1052L, 987L, 918L, 746L, 
    620L, 396L, 280L, 214L, 148L, 148L, 94L, 107L, 69L, 55L, 
    69L, 47L, 43L, 49L, 30L, 42L, 51L, 41L, 39L, 40L, 38L, 22L, 
    37L, 26L, 40L, 56L, 54L, 74L, 99L, 114L, 114L, 120L, 114L, 
    123L, 131L, 170L, 147L, 163L, 163L, 160L, 158L, 163L, 124L, 
    115L, 176L, 171L, 214L, 320L, 507L, 902L, 1190L, 1272L, 1282L, 
    1146L, 896L, 597L, 434L, 216L, 141L, 101L, 86L, 65L, 55L, 
    35L, 49L, 29L, 55L, 53L, 57L, 34L, 43L, 42L, 13L, 17L, 20L, 
    27L, 36L, 47L, 64L, 77L, 82L, 82L, 95L, 107L, 96L, 106L, 
    93L, 114L, 102L, 116L, 128L, 123L, 212L, 203L, 165L, 267L, 
    550L, 761L, 998L, 1308L, 1613L, 1704L, 1669L, 1296L, 975L, 
    600L, 337L, 259L, 145L, 91L, 70L, 79L, 63L, 58L, 51L, 53L, 
    39L, 49L, 33L, 47L, 56L, 32L, 43L, 47L, 19L, 32L, 18L, 34L, 
    39L, 63L, 57L, 55L, 69L, 76L, 103L, 99L, 108L, 131L, 113L, 
    106L, 122L, 138L, 136L, 175L, 207L, 324L, 499L, 985L, 1674L, 
    1753L, 1419L, 1105L, 821L, 466L, 274L, 180L, 143L, 82L, 101L, 
    72L, 55L, 71L, 50L, 33L, 26L, 25L, 27L, 21L, 24L, 24L, 20L, 
    18L, 18L, 25L, 23L, 13L, 10L, 16L, 9L, 12L, 16L, 25L, 31L, 
    36L, 40L, 36L, 47L, 32L, 46L, 75L, 63L, 49L, 90L, 83L, 101L, 
    78L, 79L, 98L, 131L, 83L, 122L, 179L, 334L, 544L, 656L, 718L, 
    570L, 323L, 220L, 194L, 125L, 95L, 77L, 46L, 42L, 29L, 35L, 
    21L, 29L, 16L, 14L, 19L, 15L, 19L, 18L, 21L, 10L, 14L, 7L, 
    7L, 5L, 9L, 14L, 11L, 18L, 22L, 39L, 36L, 46L, 44L, 37L, 
    30L, 39L, 37L, 45L, 71L, 59L, 57L, 80L, 68L, 88L, 72L, 74L, 
    208L, 357L, 621L, 839L, 964L, 835L, 735L, 651L, 400L, 292L, 
    198L, 85L, 64L, 41L, 40L, 23L, 18L, 14L, 22L, 9L, 19L, 8L, 
    14L, 12L, 15L, 14L, 4L, 6L, 7L, 7L, 8L, 13L, 10L, 19L, 17L, 
    20L, 22L, 40L, 37L, 45L, 34L, 26L, 35L, 67L, 49L, 77L, 82L, 
    80L, 104L, 88L, 49L, 73L, 113L, 142L, 152L, 206L, 293L, 513L, 
    657L, 919L, 930L, 793L, 603L, 323L, 202L, 112L, 55L, 31L, 
    27L, 15L, 15L, 6L, 13L, 21L, 10L, 11L, 9L, 8L, 11L, 7L, 5L, 
    1L, 4L, 7L, 2L, 6L, 12L, 14L, 21L, 29L, 32L, 26L, 22L, 44L, 
    39L, 47L, 44L, 93L, 145L, 289L, 456L, 685L, 548L, 687L, 773L, 
    575L, 355L, 248L, 179L, 129L, 122L, 103L, 72L, 72L, 36L, 
    26L, 31L, 12L, 14L, 14L, 14L, 7L, 8L, 2L, 7L, 8L, 9L, 26L, 
    10L, 13L, 13L, 5L, 5L, 3L, 6L, 1L, 10L, 6L, 7L, 17L, 12L, 
    21L, 32L, 29L, 18L, 22L, 24L, 38L, 52L, 53L, 73L, 49L, 52L, 
    70L, 77L, 95L, 135L, 163L, 303L, 473L, 823L, 1126L, 1052L, 
    794L, 459L, 314L, 252L, 111L, 55L, 35L, 14L, 30L, 21L, 16L, 
    9L, 11L, 6L, 6L, 8L, 9L, 9L, 10L, 15L, 15L, 11L, 6L, 3L, 
    8L, 4L, 7L, 7L, 13L, 10L, 23L, 24L, 36L, 25L, 34L, 37L, 46L, 
    39L, 37L, 55L, 65L, 54L, 60L, 82L, 55L, 53L, 61L, 52L, 75L, 
    92L, 121L, 170L, 199L, 231L, 259L, 331L, 357L, 262L, 154L, 
    77L, 34L, 41L, 21L, 17L, 16L, 7L, 15L, 11L, 7L, 5L, 6L, 13L, 
    7L, 6L, 8L, 7L, 1L, 11L, 9L, 3L, 9L, 9L, 8L, 15L, 19L, 16L, 
    10L, 12L, 26L, 35L, 35L, 41L, 34L, 30L, 36L, 43L, 23L, 55L, 
    107L, 141L, 217L, 381L, 736L, 782L, 663L, 398L, 182L, 137L, 
    79L, 28L, 26L, 16L, 14L, 8L, 4L, 4L, 6L, 6L, 11L, 4L, 5L, 
    7L, 7L, 6L, 8L, 2L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 8L, 8L, 11L, 
    10L, 11L, 8L, 24L, 25L, 25L, 33L, 36L, 51L, 61L, 74L, 92L, 
    89L, 123L, 402L, 602L, 524L, 494L, 406L, 344L, 329L, 225L, 
    136L, 136L, 84L, 55L, 55L, 42L, 19L, 28L, 8L, 7L, 2L, 7L, 
    6L, 4L, 3L, 5L, 3L, 3L, 0L, 1L, 2L, 3L, 2L, 1L, 2L, 2L, 9L, 
    4L, 9L, 10L, 18L, 15L, 13L, 12L, 10L, 19L, 15L, 22L, 23L, 
    34L, 43L, 53L, 47L, 57L, 328L, 552L, 787L, 736L, 578L, 374L, 
    228L, 161L, 121L, 96L, 58L, 50L, 37L, 14L, 9L, 6L, 15L, 12L, 
    9L, 1L, 6L, 4L, 7L, 7L, 3L, 6L, 9L, 15L, 22L, 28L, 34L, 62L, 
    54L, 75L, 65L, 58L, 57L, 60L, 37L, 47L, 60L, 89L, 90L, 193L, 
    364L, 553L, 543L, 676L, 550L, 403L, 252L, 140L, 125L, 99L, 
    63L, 63L, 76L, 85L, 68L, 67L, 38L, 25L, 24L, 11L, 9L, 9L, 
    4L, 8L, 4L, 6L, 5L, 2L, 6L, 4L, 4L, 1L, 5L, 4L, 1L, 2L, 2L, 
    2L, 2L, 3L, 4L, 4L, 7L, 5L, 2L, 10L, 11L, 17L, 11L, 16L, 
    15L, 11L, 12L, 21L, 20L, 25L, 46L, 51L, 90L, 123L)), .Names = c("date", 
"cases"), row.names = c(NA, -835L), class = "data.frame")

r interpolation fourier-transform COOLSerdash
источник

Старая литература часто восхищается формулами сложности рококо для такого рода вещей. На практике я бы предпочел думать об этом как о проблеме регрессии, так что интерполированные значения - это просто предсказанные значения из регрессии. См., Например, stats.stackexchange.com/questions/60500/… Основной принцип заключается в том, что основной цикл повторяется один раз в год.

Ник Кокс

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

Ник Кокс

Я вкратце опробовал это на ваших данных (используя Stata, а не R). Вкратце, хотя в ваших данных отмечена выраженная сезонность, этот подход не совсем подходит для того, чтобы этот подход вообще работал хорошо: например, не только сильно варьируется время пиков, но также и количество случаев на пике; в некоторые годы, но не во всех, наблюдается вторичный пик в конце календарного года. Кроме того, сезонность усугубляется заметной долгосрочной тенденцией. Я предполагаю, что для получения ежедневных случаев вам следует прибегнуть к некоторой строго локальной интерполяции или сглаживанию еженедельных рядов, расширенных в семь раз.

Ник Кокс

В разработке систем управления критерий Найквиста используется как основа для частоты дискретизации. Это говорит выборка на более чем в два раза большей частоте в ваших данных. На практике более обычным является сэмплирование выше 5х максимальной частоты, которую вы хотите разрешить. Если вы вводите еженедельные данные, то Найквист предполагает, что самая высокая разрешаемая частота составляет порядка 2 недель. Возможно, было бы лучше, если бы у вас была другая еженедельная статистика, чтобы сообщить выборку и подкрепить среднее значение. ru.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem

EngrStudent,

+1 Отличный вопрос! Знаете ли вы случайно, как обнаружить сигнал в шуме (предпочтительно в R), если у вас есть несколько распределений, где шум - это гауссов, а сигнал - это гауссов плюс другая часть? Я просмотрел множество пакетов и функций (signal, fft () и т. Д.) И даже поиграл с данными, пытаясь применить преобразование Фурье и даже меры энтропии, но пока безрезультатно. Я пытался ответить на вопрос (не мой) и узнать что-то новое по пути, поскольку я нахожу эту тему довольно интересной.

Александр Блех

Я не эксперт по преобразованиям Фурье, но ...

Общий диапазон выборки Эпштейна составлял 24 месяца с месячной частотой выборки: 1/12 года. Ваш выборочный диапазон составляет 835 недель. Если ваша цель - оценить среднее значение за один год с данными за ~ 16 лет на основе ежедневных данных, вам нужна частота выборки 1/365 лет. Поэтому замените 52 на 12, но сначала стандартизируйте единицы и увеличьте свои 835 недель до 835 * 7 = 5845 дней. Однако, если у вас есть только еженедельные точки данных, я предлагаю частоту выборки 52 с битовой глубиной 16 или 17 для пикового анализа, или 32 или 33 для четного / нечетного сравнения. Таким образом, параметры ввода по умолчанию включают: 1) использовать еженедельные средние значения (или срединное абсолютное отклонение, MAD или что-то в этом роде) или 2) использовать дневные значения, которые обеспечивают более высокое разрешение.

Liebman et al. выбрал точку отсечки jmax = 2. Следовательно, рис. 3. содержит меньше частичных значений и, таким образом, является более симметричным в верхней части синуса по сравнению с рис. 2. (Один частичный сигнал на базовой частоте приведет к чистой синусоидальной волне. Если бы Эпштейн выбрал более высокое разрешение (например, jmax = 12), преобразование предположительно дало бы лишь незначительные колебания с дополнительными компонентами, или, возможно, ему не хватало вычислительной мощности.

При визуальном осмотре ваших данных у вас появляется 16-17 пиков. Я бы посоветовал вам установить jmax или «битовую глубину» на 6, 11, 16 или 17 (см. Рисунок) и сравнить результаты. Чем выше пики, тем больше они вносят вклад в исходную сложную форму волны. Таким образом, предполагая 17-полосное разрешение или битовую глубину, 17-е частичное вносит минимальный вклад в исходную форму сигнала по сравнению с 6-м пиком. Тем не менее, с разрешением 34 полосы вы можете обнаружить разницу между четными и нечетными пиками в соответствии с довольно постоянными долинами. Глубина в битах зависит от вашего исследовательского вопроса, интересуетесь ли вы только пиками или пиками и долинами, а также тем, как именно вы хотите приблизиться к исходной серии.

Анализ Фурье уменьшает ваши данные. Если бы вам пришлось инвертировать функцию на определенной глубине в битах с помощью преобразования Фурье, вы, вероятно, могли бы перепроверить, соответствуют ли новые средние оценки вашим первоначальным средним значениям. Итак, чтобы ответить на ваш четвертый вопрос: упомянутые вами параметры регрессии зависят от требуемой чувствительности и разрешения. Если вы не хотите точной подгонки, то просто введите еженедельные средние в преобразовании. Однако следует помнить, что более низкая битовая глубина также уменьшает данные. Например, обратите внимание, как гармоническое наложение Эпштейна на анализ Либермана и его коллег пропускает среднюю точку ступенчатой функции с искаженной кривой немного вправо (т.е. температура слишком высокая), в декабре на рисунке 3.

Параметры Либмана и коллег:

Битовая глубина: 2

Параметры Эпштейна:

Частота дискретизации: 12 [каждый месяц]
Диапазон образцов: 24 месяца
Битовая глубина: 6

Ваши параметры:

Частота дискретизации: 365 [каждый день]
Диапазон выборки: 5845 дней

Точный разрядный подход

Точная посадка на основе визуального осмотра. (Если у вас есть сила, просто посмотрите, что происходит по сравнению с более низкой битовой глубиной.)

Полный спектр (пики): 17
Полный спектр (четный / нечетный): 34

Подход с переменной битовой глубиной

Это, вероятно, то, что вы хотите сделать:

Сравните только пики: 6, 11, 16, 17
Сравните четные / нечетные: 12, 22, 32, 34
Ресинтезировать и сравнить средства

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

ОБНОВЛЕНИЕ 1:

Бит глубина

Бит - сокращение от двоичной цифры - 0 или 1. Биты 010101 описывают прямоугольную волну. Глубина в битах составляет 1 бит. Чтобы описать пильную волну, вам понадобится больше битов: 0123210. Чем сложнее волна, тем больше бит вам нужно:

Это несколько упрощенное объяснение, но чем сложнее временной ряд, тем больше битов требуется для его моделирования. На самом деле, «1» является синусоидальной составляющей, а не прямоугольной волной (прямоугольная волна больше похожа на 3 2 1 0 - см. Прилагаемый рисунок). 0 битов будут плоской линией. Информация теряется при уменьшении битовой глубины. Например, звук CD-качества обычно 16-битный, но качество звука стационарного телефона часто составляет около 8 бит.

Пожалуйста, прочтите это изображение слева направо, сосредоточив внимание на графиках:

FFT

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

ОБНОВЛЕНИЕ 2

Фильтровать или не фильтровать

Я не совсем уверен, как вы могли бы ввести ограничение в регрессию, поскольку я знаком только с интервальными ограничениями, но, возможно, DSP - это ваше решение. Это то, что я понял до сих пор:

Шаг 1. Разбейте ряд на синусовые компоненты с помощью функции Фурье по всему набору данных (в днях)
Шаг 2. Воссоздайте временные ряды с помощью обратного преобразования Фурье, с дополнительным средним ограничением, связанным с исходными данными: отклонения интерполяций от исходных средних должны компенсировать друг друга (Harzallah, 1995).

Мое лучшее предположение состоит в том, что вам придется ввести авторегрессию, если я правильно понимаю Harzallah (1995, Рис. 2). Так что, вероятно, будет соответствовать фильтру бесконечной реакции (IIR)?

IIR http://paulbourke.net/miscellaneous/ar/

В итоге:

Получите средства из необработанных данных
Необработанные данные преобразования Фурье
Преобразованные Фурье данные, преобразованные.
Отфильтруйте результат, используя IIR

Возможно, вы могли бы использовать БИХ-фильтр, не проходя анализ Фурье? Единственное преимущество анализа Фурье, как я вижу, состоит в том, чтобы изолировать и определить, какие шаблоны являются влиятельными и как часто они повторяются (то есть колеблются). Затем вы можете отфильтровать те, которые вносят меньший вклад, например, используя узкий режекторный фильтр с наименьшим вносящим пиковым значением (или фильтр на основе ваших собственных критериев). Для начала вы можете отфильтровать менее значимые нечетные долины, которые больше похожи на шум в «сигнале». Шум характеризуется очень малым количеством случаев и отсутствием паттерна. Гребенчатый фильтр с нечетными частотными составляющими может уменьшить шум - если только вы не найдете там образец.

Вот некоторый произвольный биннинг - только для пояснительных целей: Ты видишь шум в долинах?

К сожалению, есть функция R для этого !?

При поиске IIR-фильтра я обнаружил, что функции R интерполируются в пакете Signal. Забудь все, что я сказал до этого момента. Интерполяции должны работать так же, как у Гарзаллы: http://cran.r-project.org/web/packages/signal/signal.pdf

Поиграйте с функциями. Должен сделать свое дело.

ОБНОВЛЕНИЕ 3

interp1 не interp

case.interp1 <- interp1(x=(ts.frame$no.influ.cases[!is.na(ts.frame$no.influ.case)]),y=ts.frame$yearday[!is.na(ts.frame$no.influ.case)],xi=mean(WEEKLYMEANSTABLE),method = c("cubic"))

Установите XI в исходное еженедельное средство.

ноуменальное
источник

Большое спасибо за этот ответ! Моя цель исследования проста: у меня есть еженедельные средние значения и я хочу получать ежедневные оценки, и среднее значение (интерполированных) ежедневных оценок за одну неделю должно равняться среднему недельному значению (то есть исходной точке данных). Как вы думаете, это возможно? Кроме того, я не понимаю, что означают «битовая глубина» и «пиковый анализ» (у меня нет никакого опыта с преобразованиями Фурье вообще).

COOLSerdash

@COOLSerdash посмотреть мое обновление. Да, это возможно! Но вам нужно будет найти лучший способ сравнить оценочные средние из повторно синтезированного временного ряда с исходными средними из исходного временного ряда.

Ноуменал

(Кстати: +1 завтра, я не могу больше голосовать сегодня). Большое спасибо за обновление, теперь оно стало понятнее. Я подумал о следующей процедуре: 1) подгонка функции Фурье к еженедельным средним значениям посредством регрессии, 2) использование прогнозов регрессии, чтобы «заполнить» промежутки между недельными значениями (то есть, чтобы получить дневные значения) 3) для каждую неделю рассчитывайте среднее значение всех дневных значений, и это среднее значение должно равняться исходному значению. В статье Эпштейн использовал некоторый «поправочный» коэффициент, чтобы заставить функцию иметь желаемые свойства, но я все еще не уверен, как это сделать с регрессией.

COOLSerdash

@COOLSerdash Смотрите обновление 2! Перейти к последнему абзацу.

Ноуменал

Абсолютно фантастично! Большое спасибо за исследования. Обратите внимание, что мне удалось реализовать подход Гарзаллы, используя сплайны (линейные и кубические). Так что я думаю, что мне нужно interp. Я отредактировал свой вопрос. Еще раз большое спасибо.

COOLSerdash

Фурье / тригонометрическая интерполяция

Фон

Вопросов

Ответы: