Interpolation de Fourier / trigonométrique

9

Contexte

Dans un article d'Epstein (1991): Sur l'obtention de valeurs climatologiques quotidiennes à partir de moyennes mensuelles , la formulation et un algorithme de calcul de l'interpolation de Fourier pour les valeurs périodiques et à intervalles pairs sont donnés.

Dans cet article, l'objectif est d' obtenir des valeurs quotidiennes à partir de moyennes mensuelles par interpolation.

En bref, on suppose que des valeurs quotidiennes inconnues peuvent être représentées par la somme des composantes harmoniques: Dans l'article, le (temps) est exprimé en mois. t

y(t)=a0+j[ajcos(2πjt/12)+bjpéché(2πjt/12)]
t

Après une certaine déviation, il est montré que les termes peuvent être calculés par: Où désigne les moyennes mensuelles et le mois.

une0=TOuiT/12unej=[(πj/12)/péché(πj/12)]×T[OuiTcos(2πjT/12)/6]       j=1,,5bj=[(πj/12)/péché(πj/12)]×T[OuiTpéché(2πjT/12)/6]       j=1,,5une6=[(πj/12)/péché(πj/12)]×T[OuiTcos(πT)/12]b6=0
OuiTT

Harzallah (1995) résume cette approche comme suit: "L'interpolation est effectuée en ajoutant des zéros aux coefficients spectraux des données et en effectuant une transformée de Fourier inverse aux coefficients étendus résultants. La méthode équivaut à appliquer un filtre rectangulaire aux coefficients de Fourier. . "


Des questions

Mon objectif est d'utiliser la méthodologie ci-dessus pour l'interpolation des moyennes hebdomadaires pour obtenir des données quotidiennes (voir ma question précédente ). En résumé, j'ai 835 moyennes hebdomadaires de données de comptage (voir l'exemple de jeu de données au bas de la question). Il y a pas mal de choses que je ne comprends pas avant de pouvoir appliquer l'approche décrite ci-dessus:

  1. Comment les formules devraient-elles être modifiées pour ma situation (valeurs hebdomadaires au lieu de mensuelles)?
  2. Comment le temps pourrait être exprimé? J'ai supposé (ou avec points de données en général), est-ce correct?tt/835t/nn
  3. Pourquoi l'auteur calcule-t-il 7 termes (ie )? Combien de termes devrais-je considérer?0j6
  4. Je comprends que la question peut probablement être résolue en utilisant une approche de régression et en utilisant les prédictions d'interpolation (merci à Nick). Cependant, certaines choses ne sont pas claires pour moi: combien de termes d'harmoniques devraient être inclus dans la régression? Et quelle période dois-je prendre? Comment faire la régression pour s'assurer que les moyennes hebdomadaires sont préservées (car je ne veux pas un ajustement harmonique exact aux données)?

En utilisant l' approche de régression (qui est également expliquée dans cet article ), j'ai réussi à obtenir un ajustement harmonique exact des données (le dans mon exemple parcourrait , j'ai donc ajusté 417 termes). Comment modifier cette approche - si possible - pour assurer la conservation des moyens hebdomadaires? Peut-être en appliquant des facteurs de correction à chaque terme de régression?j1,,417  

Le tracé de l'ajustement harmonique exact est:

Ajustement harmonique exact

ÉDITER

En utilisant le paquet de signaux et la interp1fonction, voici ce que j'ai réussi à faire en utilisant l'exemple de jeu de données ci-dessous (merci beaucoup à @noumenal). J'utilise q=7car nous avons des données hebdomadaires:

# 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

Il y a deux problèmes ici:

  1. La courbe semble correspondre "trop ​​bien": elle passe par tous les points
  2. Les moyennes hebdomadaires ne sont pas conservées

Exemple de jeu de données

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")
COOLSerdash
la source
2
La littérature plus ancienne se délecte souvent de formules de complexité rococo pour ce genre de choses. En pratique, je pencherais en arrière pour penser à cela comme un problème de régression, de sorte que les valeurs interpolées ne sont que les valeurs prédites de la régression. Voir par exemple stats.stackexchange.com/questions/60500/… Le principe clé est que le cycle principal se répète une fois par an.
Nick Cox
1
En pratique, vous voulez un ajustement très proche des données, car vous voulez un lissage local. Cela peut nécessiter plusieurs paires de Fourier, mais les rendements décroissants s'installent très rapidement, de sorte que chaque nouvelle paire de termes ajoute bientôt très peu. Il suffit de le sucer et de voir. Tout tracer le rend plus clair.
Nick Cox
1
J'ai essayé brièvement avec vos données (en utilisant Stata, pas R). En bref, bien qu'il y ait une saisonnalité marquée dans vos données, elle n'est pas assez régulière pour que cette approche fonctionne bien du tout: par exemple, non seulement le moment des pics varie beaucoup, mais aussi le nombre de cas au sommet; certaines années, mais pas toutes, il y a un pic secondaire tard dans l'année civile. En outre, la saisonnalité est aggravée par une tendance marquée à long terme. Je suppose que pour obtenir des cas quotidiens, vous devez recourir à une interpolation strictement locale ou à un lissage de la série hebdomadaire multipliée par sept.
Nick Cox
1
Dans l'ingénierie des systèmes de contrôle, le critère de Nyquist est utilisé comme plancher pour le taux d'échantillonnage. Il indique que l'échantillon est plus de deux fois la fréquence la plus élevée de vos données. En pratique, il est plus conventionnel d'échantillonner au-dessus de 5 fois la fréquence la plus élevée que vous espérez résoudre. Si vos données sont des données hebdomadaires, Nyquist suggère que la fréquence résolvable la plus élevée est de l'ordre de 2 semaines. Il serait préférable que vous disposiez d'autres statistiques hebdomadaires pour informer l'échantillonnage et étayer la moyenne. en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem
EngrStudent
+1 Excellente question! Savez-vous par hasard comment détecter le signal dans le bruit (de préférence dans R), à condition que vous ayez un certain nombre de distributions, où le bruit est un gaussien et le signal est un gaussien plus une autre partie? J'ai regardé de nombreux packages et fonctions (signal, fft (), etc.) et j'ai même joué avec des données, en essayant d'appliquer la transformation de Fourier et même des mesures d'entropie, mais en vain jusqu'à présent. J'essayais de répondre à une question (pas la mienne) et d'apprendre quelque chose de nouveau en cours de route, car je trouve le sujet très intéressant.
Aleksandr Blekh

Réponses:

4

Je ne suis pas un expert des transformées de Fourier, mais ...

La plage d'échantillonnage totale d'Epstein était de 24 mois avec un taux d'échantillonnage mensuel: 1/12 ans. Votre plage d'échantillonnage est de 835 semaines. Si votre objectif est d'estimer la moyenne sur un an avec des données de ~ 16 ans basées sur des données quotidiennes, vous avez besoin d'un taux d'échantillonnage de 1/365 ans. Remplacez donc 52 par 12, mais standardisez d'abord les unités et augmentez vos 835 semaines à 835 * 7 = 5845 jours. Cependant, si vous ne disposez que de points de données hebdomadaires, je suggère un taux d'échantillonnage de 52 avec une profondeur de bits de 16 ou 17 pour l'analyse des pics, ou 32 ou 33 pour la comparaison paire / impaire. Ainsi, les options d'entrée par défaut incluent: 1) pour utiliser les moyennes hebdomadaires (ou l'écart absolu médian, MAD, ou quelque chose dans cette mesure) ou 2) pour utiliser les valeurs quotidiennes, qui fournissent une résolution plus élevée.

Liebman et al. a choisi le point de coupure jmax = 2. Par conséquent, la figure 3. contient moins de partiels et est donc plus symétrique en haut du sinus par rapport à la figure 2. (Un seul partiel à la fréquence de base entraînerait une onde sinusoïdale pure. ) Si Epstein avait choisi une résolution plus élevée (par exemple jmax = 12), la transformée ne produirait vraisemblablement que des fluctuations mineures avec les composants supplémentaires, ou peut-être manquait-il de la puissance de calcul.

Grâce à l'inspection visuelle de vos données, vous semblez avoir 16-17 pics. Je vous suggère de régler jmax ou la "profondeur de bits" sur 6, 11, 16 ou 17 (voir figure) et de comparer les sorties. Plus les pics sont élevés, plus ils contribuent à la forme d'onde complexe d'origine. Ainsi, en supposant une résolution de 17 bandes ou une profondeur de bits, le 17e partiel contribue de façon minimale au modèle de forme d'onde d'origine par rapport au 6e pic. Cependant, avec une résolution de 34 bandes, vous détecteriez une différence entre les pics pairs et impairs comme le suggèrent les vallées assez constantes. La profondeur de bits dépend de votre question de recherche, si vous êtes intéressé par les pics uniquement ou à la fois par les pics et les vallées, mais aussi par la façon dont vous souhaitez approximer exactement la série d'origine.

L'analyse de Fourier réduit vos points de données. Si vous deviez inverser la fonction à une certaine profondeur de bits en utilisant une transformée de Fourier, vous pourriez probablement vérifier si les nouvelles estimations moyennes correspondent à vos moyennes originales. Donc, pour répondre à votre quatrième question: les paramètres de régression que vous avez mentionnés dépendent de la sensibilité et de la résolution dont vous avez besoin. Si vous ne souhaitez pas un ajustement exact, entrez simplement les moyens hebdomadaires dans la transformation. Cependant, sachez qu'une profondeur de bits plus faible réduit également les données. Par exemple, notez comment la superposition harmonique d'Epstein sur l'analyse de Lieberman et ses collègues manque le point médian de la fonction de pas, avec une courbe asymétrique légèrement à droite (c.-à-d. Température est trop élevée), en décembre sur la figure 3.


Paramètres de Liebman et de ses collègues:

  • Profondeur de bits: 2

Paramètres d'Epstein:

  • Taux d'échantillonnage: 12 [chaque mois]
  • Gamme d'échantillons: 24 mois
  • Profondeur de bits: 6

Vos paramètres:

  • Taux d'échantillonnage: 365 [tous les jours]

  • Plage d'échantillonnage: 5845 jours

Approche exacte de la profondeur de bits

Ajustement exact basé sur une inspection visuelle. (Si vous avez le pouvoir, voyez ce qui se passe par rapport à des profondeurs de bits plus faibles.)

  • Spectre complet (pics): 17
  • Spectre complet (pair / impair): 34

Approche à profondeur de bits variable

C'est probablement ce que vous souhaitez faire:

  • Comparer les pics uniquement: 6, 11, 16, 17
  • Comparer les nombres pairs / impairs: 12, 22, 32, 34
  • Resynthétiser et comparer les moyens

Cette approche produirait quelque chose de similaire à la comparaison des figures dans Epstein si vous inversez à nouveau la transformation, c'est-à-dire synthétisez les partiels en une approximation de la série temporelle d'origine. Vous pouvez également comparer les points discrets des courbes resynthétisées aux valeurs moyennes, peut-être même tester des différences significatives pour indiquer la sensibilité de votre choix de profondeur de bits.


MISE À JOUR 1:

Peu profond

Un bit - abréviation de chiffre binaire - est soit 0 soit 1. Les bits 010101 décriraient une onde carrée. La profondeur de bits est de 1 bit. Pour décrire une onde de scie, vous auriez besoin de plus de bits: 0123210. Plus une onde est complexe, plus vous avez besoin de bits:

Il s'agit d'une explication quelque peu simplifiée, mais plus une série chronologique est complexe, plus il faut de bits pour la modéliser. En fait, "1" est une composante d'onde sinusoïdale et non une onde carrée (une onde carrée ressemble plus à 3 2 1 0 - voir figure ci-jointe). 0 bits serait une ligne plate. Les informations sont perdues avec la réduction de la profondeur de bits. Par exemple, l'audio de qualité CD est généralement de 16 bits, mais l'audio de qualité de téléphone fixe est souvent d'environ 8 bits.

Veuillez lire cette image de gauche à droite, en vous concentrant sur les graphiques:

FFT

Vous venez de terminer une analyse du spectre de puissance (bien qu'à haute résolution dans votre figure). Votre prochain objectif serait de déterminer: de combien de composants ai-je besoin dans le spectre de puissance pour capturer avec précision les moyennes des séries temporelles?

MISE À JOUR 2


Filtrer ou ne pas filtrer

Je ne suis pas tout à fait sûr de la manière dont vous introduiriez la contrainte dans la régression car je ne connais que les contraintes d'intervalle, mais peut-être que DSP est votre solution. Voici ce que je pensais jusqu'à présent:

  • Étape 1. Décomposer la série en composantes sinusales via la fonction de Fourier sur l'ensemble de données complet (en jours)

  • Étape 2. Recréez la série chronologique par une transformée de Fourier inverse, avec la contrainte moyenne supplémentaire couplée aux données originales: les écarts des interpolations par rapport aux moyennes originales devraient s'annuler (Harzallah, 1995).

Ma meilleure supposition est que vous devrez introduire une autorégression si je comprends bien Harzallah (1995, Fig 2). Cela correspondrait donc probablement à un filtre à réponse infinie (IIR)?

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

En résumé:

  1. Dériver les moyennes des données brutes
  2. Transformation de Fourier Données brutes
  3. Transformation inverse de Fourier transformée.
  4. Filtrer le résultat avec IIR

Peut-être pourriez-vous utiliser un filtre IIR sans passer par l'analyse de Fourier? Le seul avantage de l'analyse de Fourier telle que je la vois est d'isoler et de déterminer quels modèles sont influents et à quelle fréquence ils se reproduisent (c'est-à-dire oscillent). Vous pouvez alors décider de filtrer ceux qui contribuent le moins, par exemple en utilisant un filtre coupe-bande étroit au pic le moins contributeur (ou un filtre basé sur vos propres critères). Pour commencer, vous pouvez filtrer les vallées impaires moins contributives qui ressemblent davantage à du bruit dans le "signal". Le bruit est caractérisé par très peu de cas et aucun motif. Un filtre en peigne à composantes de fréquence impaires pourrait réduire le bruit - à moins que vous ne trouviez un motif à cet endroit.

Voici quelques regroupements arbitraires, à des fins explicatives uniquement: Pouvez-vous voir le bruit dans les vallées?

Oups - il y a une fonction R pour ça!?

Lors de la recherche d'un filtre IIR, j'ai découvert les fonctions R interpolées dans le paquet Signal. Oubliez tout ce que j'ai dit jusqu'ici. Les interpolations devraient fonctionner comme celles de Harzallah: http://cran.r-project.org/web/packages/signal/signal.pdf

Jouez avec les fonctions. Devrait faire l'affaire.


MISE À JOUR 3

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

Réglez xi sur les moyennes hebdomadaires d'origine.

noumenal
la source
Merci beaucoup pour cette réponse! Mon objectif de recherche est simple: j'ai des moyennes hebdomadaires et je veux obtenir des estimations quotidiennes et la moyenne des estimations quotidiennes (interpolées) sur une semaine doit être égale à la moyenne hebdomadaire (c'est-à-dire le point de données d'origine). Pensez-vous que c'est possible?. De plus, je ne comprends pas ce que "profondeur de bit" et "analyse de pic" signifient (je n'ai aucune expérience avec les transformations de Fourier que ce soit).
COOLSerdash
1
@COOLSerdash voir ma mise à jour. Oui c'est possible! Mais vous devrez trouver la meilleure façon de comparer les moyennes estimées de la série chronologique resynthétisée avec les moyennes originales de la série temporelle d'origine.
noumenal
(Btw: +1 demain, je ne peux plus voter aujourd'hui). Merci beaucoup pour la mise à jour, c'est plus clair maintenant. J'ai pensé à la procédure suivante: 1) ajuster une fonction de Fourier aux moyennes hebdomadaires par régression, 2) utiliser les prédictions de la régression pour "combler" les écarts entre les valeurs hebdomadaires (c'est-à-dire pour obtenir les valeurs quotidiennes) 3) pour chaque semaine, calculez la moyenne de toutes les valeurs quotidiennes et cette moyenne doit être égale à la valeur d'origine. Dans l'article, Epstein a utilisé une sorte de facteur de "correction" pour forcer la fonction à avoir les propriétés souhaitées, mais je ne sais toujours pas comment faire cela avec la régression.
COOLSerdash
@COOLSerdash Voir la mise à jour 2! Passez au dernier paragraphe.
noumenal
Absolument fantastique! Merci beaucoup pour la recherche. Notez que j'ai réussi à implémenter l'approche de Harzallah en utilisant des splines (linéaires et cubiques). Donc je suppose que j'aurais besoin interp. J'ai édité ma question. Merci encore beaucoup.
COOLSerdash