Fourier / trigonometrische Interpolation

Hintergrund

In einer Arbeit von Epstein (1991): Um tägliche klimatologische Werte aus monatlichen Mitteln zu erhalten , werden die Formulierung und ein Algorithmus zur Berechnung der Fourier-Interpolation für periodische und geradzahlige Werte angegeben.

In der Arbeit geht es darum , durch Interpolation Tageswerte aus monatlichen Mitteln zu erhalten .

Kurz gesagt wird angenommen, dass unbekannte Tageswerte durch die Summe der harmonischen Komponenten dargestellt werden können:

y (t) = {ein}_{0} + \sum_{j} [{ein}_{j} \cos (2 π j t /. 12) + b_{j} Sünde (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]$ In der Arbeit wird

t

$t$ (Zeit) in Monaten ausgedrückt.

Nach einiger Abweichung wird gezeigt, dass die Terme berechnet werden können durch: Wobeidas Monatsmittel undden Monat bezeichnet.

\begin{aligned} a_{0} & = \sum_{T} Y_{T} / 12 \\ a_{j} & = [(π j / 12) / \sin (π j / 12)] \times \sum_{T} [Y_{T} \cos (2 π j T / 12) / 6] j = 1, \dots, 5 \\ b_{j} & = [(π j / 12) / \sin (π j / 12)] \times \sum_{T} [Y_{T} \sin (2 π j T / 12) / 6] j = 1, \dots, 5 \\ a_{6} & = [(π j / 12) / \sin (π j / 12)] \times \sum_{T} [Y_{T} \cos (π T) / 12] \\ b_{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) fasst diesen Ansatz wie folgt zusammen: "Die Interpolation wird durchgeführt, indem Nullen zu den Spektralkoeffizienten von Daten hinzugefügt werden und indem eine inverse Fourier-Transformation zu den resultierenden erweiterten Koeffizienten durchgeführt wird. Das Verfahren entspricht der Anwendung eines Rechteckfilters auf Fourier-Koeffizienten . "

Fragen

Mein Ziel ist es, die obige Methode zur Interpolation wöchentlicher Mittelwerte zu verwenden, um tägliche Daten zu erhalten (siehe meine vorherige Frage ). Zusammenfassend habe ich 835 wöchentliche Zählmittel (siehe Beispieldatensatz am Ende der Frage). Es gibt einige Dinge, die ich nicht verstehe, bevor ich den oben beschriebenen Ansatz anwenden kann:

Wie müssten die Formeln für meine Situation geändert werden (wöchentliche statt monatliche Werte)?
Wie könnte die Zeit ausgedrückt werden? Ich habe (oder mit Datenpunkten im Allgemeinen), ist das richtig? $t$ $t/835$ $t/n$ $n$
Warum berechnet der Autor 7 Terme (dh )? Wie viele Begriffe müsste ich berücksichtigen? $0\leq j \leq 6$
Ich verstehe, dass die Frage wahrscheinlich durch einen Regressionsansatz und die Vorhersagen für die Interpolation gelöst werden kann (danke an Nick). Dennoch sind mir einige Dinge unklar: Wie viele Harmonische sollten in die Regression einbezogen werden? Und welche Zeit soll ich nehmen? Wie kann die Regression durchgeführt werden, um sicherzustellen, dass die wöchentlichen Mittelwerte erhalten bleiben (da ich keine exakte harmonische Anpassung an die Daten möchte)?

$j$ $1, \ldots, 417$ $~$ $~$

Die Darstellung der exakten harmonischen Anpassung lautet:

Genaue harmonische Passform

BEARBEITEN

Mit dem Signalpaket und der interp1Funktion habe ich Folgendes mithilfe des Beispieldatensatzes von unten erreicht (vielen Dank an @noumenal). Ich benutze, q=7da wir wöchentliche Daten haben:

# 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

Hier gibt es zwei Probleme:

Die Kurve scheint "zu gut" zu passen: Sie geht durch jeden Punkt
Die wöchentlichen Mittel bleiben nicht erhalten

Beispieldatensatz

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
quelle

Ältere Literatur erfreut sich oft an Formeln von Rokoko-Komplexität für solche Dinge. In der Praxis würde ich mich zurücklehnen, um dies als ein Regressionsproblem zu betrachten, so dass die interpolierten Werte nur die vorhergesagten Werte aus der Regression sind. Siehe z. B. stats.stackexchange.com/questions/60500/…. Das Schlüsselprinzip ist, dass sich der Hauptzyklus einmal pro Jahr wiederholt.

Nick Cox

In der Praxis möchten Sie eine sehr enge Anpassung an die Daten, da Sie eine lokale Glättung wünschen. Das könnte mehrere Fourier-Paare erfordern, aber die sinkenden Renditen setzen sehr schnell ein, so dass jedes neue Begriffspaar bald sehr wenig hinzufügt. Sie müssen es nur saugen und sehen. Alles zu zeichnen macht es klarer.

Nick Cox

Ich habe dies kurz mit Ihren Daten ausprobiert (mit Stata, nicht mit R). Kurz gesagt, obwohl Ihre Daten eine ausgeprägte Saisonalität aufweisen, ist dies nicht regelmäßig genug, damit dieser Ansatz überhaupt gut funktioniert: Beispielsweise variiert nicht nur das Timing der Peaks stark, sondern auch die Anzahl der Fälle am Peak. In einigen Jahren, aber nicht in allen, gibt es spät im Kalenderjahr einen sekundären Höhepunkt. Darüber hinaus wird die Saisonalität durch einen deutlichen langfristigen Trend verstärkt. Ich vermute, um tägliche Fälle zu erhalten, sollten Sie auf eine streng lokale Interpolation oder Glättung der siebenfach erweiterten wöchentlichen Reihe zurückgreifen.

Nick Cox

In der Steuerungssystemtechnik wird das Nyquist-Kriterium als Grundlage für die Abtastrate verwendet. Es heißt Probe mit mehr als der doppelten höchsten Frequenz in Ihren Daten. In der Praxis ist es üblicher, über 5x die höchste Frequenz abzutasten, die Sie auflösen möchten. Wenn Ihre Eingabe wöchentliche Daten sind, schlägt Nyquist vor, dass die höchste auflösbare Häufigkeit in der Größenordnung von 2 Wochen liegt. Es könnte besser sein, wenn Sie andere wöchentliche Statistiken hätten, um die Stichprobe zu informieren und den Mittelwert zu stützen. en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem

EngrStudent

+1 Ausgezeichnete Frage! Wissen Sie zufällig, wie Sie ein Signal im Rauschen erkennen (vorzugsweise in R), vorausgesetzt, Sie haben eine Reihe von Verteilungen, bei denen das Rauschen ein Gaußscher Wert und das Signal ein Gaußscher Wert plus ein weiterer Teil ist? Ich habe mir viele Pakete und Funktionen (Signal, fft () usw.) angesehen und sogar mit Daten gespielt und versucht, Fourier-Transformation und sogar Entropiemaßnahmen anzuwenden, aber bisher ohne Erfolg. Ich habe versucht, eine Frage zu beantworten (nicht meine) und dabei etwas Neues zu lernen, da ich das Thema sehr interessant finde.

Aleksandr Blekh

Ich bin kein Experte für Fourier-Transformationen, aber ...

Der gesamte Stichprobenbereich von Epstein betrug 24 Monate mit einer monatlichen Stichprobenrate von 1/12 Jahren. Ihr Probenumfang beträgt 835 Wochen. Wenn Sie den Durchschnitt für ein Jahr mit Daten aus ~ 16 Jahren basierend auf täglichen Daten schätzen möchten, benötigen Sie eine Abtastrate von 1/365 Jahren. Ersetzen Sie also 12 durch 52, standardisieren Sie jedoch zuerst die Einheiten und erweitern Sie Ihre 835 Wochen auf 835 * 7 = 5845 Tage. Wenn Sie jedoch nur wöchentliche Datenpunkte haben, empfehle ich eine Abtastrate von 52 mit einer Bittiefe von 16 oder 17 für die Peakanalyse, alternativ 32 oder 33 für einen geraden / ungeraden Vergleich. Die Standardeingabeoptionen umfassen also: 1) Verwenden des wöchentlichen Mittels (oder der mittleren absoluten Abweichung, MAD oder etwas in diesem Ausmaß) oder 2) Verwenden der Tageswerte, die eine höhere Auflösung bieten.

Liebman et al. wählte den Grenzpunkt jmax = 2. Daher enthält Fig. 3 weniger Teiltöne und ist daher am oberen Ende des Sinus symmetrischer als Fig. 2. (Ein einzelner Teil bei der Grundfrequenz würde zu einer reinen Sinuswelle führen. ) Wenn Epstein eine höhere Auflösung gewählt hätte (z. B. jmax = 12), würde die Transformation vermutlich nur geringfügige Schwankungen mit den zusätzlichen Komponenten ergeben, oder ihm fehlte möglicherweise die Rechenleistung.

Durch visuelle Überprüfung Ihrer Daten scheinen Sie 16-17 Peaks zu haben. Ich würde vorschlagen, dass Sie jmax oder die "Bittiefe" auf 6, 11, 16 oder 17 setzen (siehe Abbildung) und die Ausgänge vergleichen. Je höher die Peaks sind, desto mehr tragen sie zur ursprünglichen komplexen Wellenform bei. Unter der Annahme einer 17-Band-Auflösung oder Bittiefe trägt der 17. Teil im Vergleich zum 6. Peak minimal zum ursprünglichen Wellenformmuster bei. Bei einer Auflösung von 34 Bändern würden Sie jedoch einen Unterschied zwischen geraden und ungeraden Spitzen feststellen, wie dies durch die ziemlich konstanten Täler nahegelegt wird. Die Bittiefe hängt von Ihrer Forschungsfrage ab, ob Sie sich nur für die Peaks oder sowohl für Peaks als auch für Täler interessieren, aber auch, wie genau Sie die ursprüngliche Serie approximieren möchten.

Die Fourier-Analyse reduziert Ihre Datenpunkte. Wenn Sie die Funktion bei einer bestimmten Bittiefe mithilfe einer Fourier-Transformation invertieren, können Sie wahrscheinlich überprüfen, ob die neuen Mittelwertschätzungen Ihren ursprünglichen Mitteln entsprechen. Um Ihre vierte Frage zu beantworten: Die von Ihnen genannten Regressionsparameter hängen von der Empfindlichkeit und Auflösung ab, die Sie benötigen. Wenn Sie keine genaue Anpassung wünschen, geben Sie auf jeden Fall einfach die wöchentlichen Mittelwerte in die Transformation ein. Beachten Sie jedoch, dass eine geringere Bittiefe auch die Daten reduziert. Beachten Sie beispielsweise, dass Epsteins harmonische Überlagerung der Analyse von Lieberman und Kollegen im Dezember in Abbildung 3 den Mittelpunkt der Schrittfunktion mit einer leicht nach rechts geneigten Kurve (dh eine zu hohe Temperatur) verfehlt.

Parameter von Liebman und Kollegen:

Bittiefe: 2

Epsteins Parameter:

Abtastrate: 12 [jeden Monat]
Probenbereich: 24 Monate
Bittiefe: 6

Ihre Parameter:

Abtastrate: 365 [jeden Tag]
Probenbereich: 5845 Tage

Genaue Bittiefenannäherung

Genaue Passform basierend auf Sichtprüfung. (Wenn Sie die Leistung haben, sehen Sie einfach, was im Vergleich zu niedrigeren Bittiefen passiert.)

Vollspektrum (Peaks): 17
Volles Spektrum (gerade / ungerade): 34

Ansatz mit variabler Bittiefe

Dies ist wahrscheinlich das, was Sie tun möchten:

Nur Peaks vergleichen: 6, 11, 16, 17
Vergleiche gerade / ungerade: 12, 22, 32, 34
Mittel neu synthetisieren und vergleichen

Dieser Ansatz würde etwas Ähnliches wie der Vergleich der Figuren in Epstein ergeben, wenn Sie die Transformation erneut umkehren, dh die Teiltöne in eine Annäherung an die ursprüngliche Zeitreihe synthetisieren. Sie können auch die diskreten Punkte der resynthetisierten Kurven mit den Mittelwerten vergleichen und möglicherweise sogar auf signifikante Unterschiede testen, um die Empfindlichkeit Ihrer Bittiefenauswahl anzuzeigen.

UPDATE 1:

Bittiefe

Ein Bit - kurz für Binärziffer - ist entweder 0 oder 1. Die Bits 010101 würden eine Rechteckwelle beschreiben. Die Bittiefe beträgt 1 Bit. Um eine Sägewelle zu beschreiben, benötigen Sie mehr Bits: 0123210. Je komplexer eine Welle wird, desto mehr Bits benötigen Sie:

Dies ist eine etwas vereinfachte Erklärung, aber je komplexer eine Zeitreihe ist, desto mehr Bits sind erforderlich, um sie zu modellieren. Tatsächlich ist "1" eine Sinuswellenkomponente und keine Rechteckwelle (eine Rechteckwelle ähnelt eher 3 2 1 0 - siehe beigefügte Abbildung). 0 Bits wären eine flache Linie. Informationen gehen mit der Verringerung der Bittiefe verloren. Beispielsweise beträgt Audio in CD-Qualität normalerweise 16 Bit, Audio in Festnetz-Telefonqualität jedoch häufig etwa 8 Bit.

Bitte lesen Sie dieses Bild von links nach rechts und konzentrieren Sie sich dabei auf die Grafiken:

FFT

Sie haben gerade eine Leistungsspektrumanalyse abgeschlossen (obwohl in Ihrer Abbildung mit hoher Auflösung). Ihr nächstes Ziel wäre es herauszufinden: Wie viele Komponenten benötige ich im Leistungsspektrum, um die Mittelwerte der Zeitreihen genau zu erfassen?

UPDATE 2

Filtern oder nicht filtern

Ich bin nicht ganz sicher, wie Sie die Einschränkung in die Regression einführen würden, da ich nur mit Intervallbeschränkungen vertraut bin, aber vielleicht ist DSP Ihre Lösung. Das habe ich mir bisher gedacht:

Schritt 1. Zerlegen Sie die Reihe durch die Fourier-Funktion in Sinuskomponenten für den gesamten Datensatz (in Tagen).
Schritt 2. Erstellen Sie die Zeitreihen durch eine inverse Fourier-Transformation neu, wobei die zusätzliche Mittelwertbeschränkung an die ursprünglichen Daten gekoppelt ist: Die Abweichungen der Interpolationen von den ursprünglichen Mitteln sollten sich gegenseitig aufheben (Harzallah, 1995).

Meine beste Vermutung ist, dass Sie eine Autoregression einführen müssten, wenn ich Harzallah (1995, Abb. 2) richtig verstehe. Das würde also wahrscheinlich einem Infinite-Response-Filter (IIR) entsprechen?

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

In Summe:

Mittel aus Rohdaten ableiten
Fourier-Transformations-Rohdaten
Fourier Inverse Transform transformierte Daten.
Filtern Sie das Ergebnis mit IIR

Vielleicht könnten Sie einen IIR-Filter verwenden, ohne die Fourier-Analyse zu durchlaufen? Der einzige Vorteil der Fourier-Analyse besteht meines Erachtens darin, zu isolieren und zu bestimmen, welche Muster Einfluss haben und wie oft sie erneut auftreten (dh oszillieren). Sie können dann entscheiden, diejenigen herauszufiltern, die weniger beitragen, z. B. einen Filter mit schmaler Kerbe am wenigsten beitragenden Peak (oder Filter nach Ihren eigenen Kriterien). Für den Anfang könnten Sie die weniger beitragenden ungeraden Täler herausfiltern, die im "Signal" eher wie Rauschen erscheinen. Lärm ist durch sehr wenige Fälle und kein Muster gekennzeichnet. Ein Kammfilter bei ungeraden Frequenzkomponenten kann das Rauschen reduzieren - es sei denn, Sie finden dort ein Muster.

Hier ist ein beliebiges Binning - nur zu Erklärungszwecken: Kannst du den Lärm in den Tälern sehen?

Ups - dafür gibt es eine R-Funktion!?

Bei der Suche nach einem IIR-Filter habe ich zufällig festgestellt, dass die R-Funktionen im Signalpaket interpolieren . Vergiss alles, was ich bis jetzt gesagt habe. Die Interpolationen sollten wie bei Harzallah funktionieren: http://cran.r-project.org/web/packages/signal/signal.pdf

Spielen Sie mit den Funktionen herum. Sollte den Trick machen.

UPDATE 3

interp1 nicht 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"))

Stellen Sie xi auf das ursprüngliche wöchentliche Mittel ein.

noumenal
quelle

Vielen Dank für diese Antwort! Mein Forschungsziel ist einfach: Ich habe wöchentliche Mittelwerte und möchte tägliche Schätzungen erhalten, und der Mittelwert der (interpolierten) täglichen Schätzungen über eine Woche sollte dem wöchentlichen Mittelwert (dh dem ursprünglichen Datenpunkt) entsprechen. Glaubst du, das ist möglich? Außerdem verstehe ich nicht, was "Bittiefe" und "Peakanalyse" bedeuten (ich habe keinerlei Erfahrung mit Fourier-Transformationen).

COOLSerdash

@COOLSerdash siehe mein Update. Ja, es ist möglich! Sie müssten jedoch herausfinden, wie Sie die geschätzten Mittelwerte aus der resynthetisierten Zeitreihe am besten mit den ursprünglichen Mittelwerten aus der ursprünglichen Zeitreihe vergleichen können.

noumenal

(Übrigens: +1 morgen, ich kann heute nicht mehr wählen). Vielen Dank für das Update, es ist jetzt klarer. Ich dachte über das folgende Verfahren nach: 1) eine Fourier-Funktion durch Regression an die wöchentlichen Mittelwerte anpassen, 2) die Vorhersagen der Regression verwenden, um die Lücken zwischen den wöchentlichen Werten zu "füllen" (dh um die täglichen Werte zu erhalten) 3) für Berechnen Sie jede Woche den Mittelwert aller Tageswerte. Dieser Mittelwert sollte dem ursprünglichen Wert entsprechen. In der Arbeit verwendete Epstein eine Art "Korrektur" -Faktor, um die Funktion zu zwingen, die gewünschten Eigenschaften zu haben, aber ich bin mir immer noch nicht sicher, wie ich das mit Regression machen soll.

COOLSerdash

@COOLSerdash Siehe Update 2! Fahren Sie mit dem letzten Absatz fort.

noumenal

Absolut fantastisch! Vielen Dank für die Recherche. Beachten Sie, dass ich es geschafft habe , Harzallahs Ansatz mithilfe von Splines (linear und kubisch) zu implementieren. Also würde ich wohl brauchen interp. Ich habe meine Frage bearbeitet. Nochmals vielen Dank.

COOLSerdash

Fourier / trigonometrische Interpolation

Hintergrund

Fragen

Antworten: