Interpolasi Fourier / trigonometri

9

Latar Belakang

Dalam sebuah makalah dari Epstein (1991): Pada memperoleh nilai-nilai klimatologis harian dari rata-rata bulanan , formulasi dan algoritma untuk menghitung interpolasi Fourier untuk nilai periodik dan genap spasi diberikan.

Di koran, tujuannya adalah untuk mendapatkan nilai harian dari sarana bulanan dengan interpolasi.

Singkatnya, diasumsikan bahwa nilai harian yang tidak diketahui dapat diwakili oleh jumlah komponen harmonik: Dalam makalah (waktu) dinyatakan dalam bulan. t

y(t)=Sebuah0+j[Sebuahjcos(2πjt/12)+bjdosa(2πjt/12)]
t

Setelah beberapa penurunan, ditunjukkan bahwa persyaratan dapat dihitung dengan: Di mana menunjukkan sarana bulanan dan bulan.

a0=TYT/12aj=[(πj/12)/sin(πj/12)]×T[YTcos(2πjT/12)/6]       j=1,...,5bj=[(πj/12)/dosa(πj/12)]×T[YTdosa(2πjT/12)/6]       j=1,...,5Sebuah6=[(πj/12)/dosa(πj/12)]×T[YTcos(πT)/12]b6=0
YTT

Harzallah (1995) merangkum pendekatan ini sebagai berikut: "Interpolasi dilakukan dengan menambahkan nol pada koefisien spektral data dan dengan melakukan transformasi Fourier terbalik ke koefisien yang dihasilkan. Metode ini setara dengan menerapkan filter persegi panjang ke koefisien Fourier . "


Pertanyaan

Tujuan saya adalah menggunakan metodologi di atas untuk interpolasi sarana mingguan untuk mendapatkan data harian (lihat pertanyaan saya sebelumnya ). Singkatnya, saya memiliki rata-rata 835 data mingguan (lihat contoh dataset di bagian bawah pertanyaan). Ada beberapa hal yang tidak saya mengerti sebelum saya dapat menerapkan pendekatan yang diuraikan di atas:

  1. Bagaimana rumus harus diubah untuk situasi saya (nilai mingguan dan bukan bulanan)?
  2. Bagaimana waktu diungkapkan? Saya berasumsi (atau dengan titik data secara umum), apakah itu benar?tt/835t/nn
  3. Mengapa penulis menghitung 7 istilah (yaitu )? Berapa banyak syarat yang harus saya pertimbangkan?0j6
  4. Saya mengerti bahwa pertanyaannya mungkin dapat diselesaikan dengan menggunakan pendekatan regresi dan menggunakan prediksi untuk interpolasi (terima kasih kepada Nick). Namun, ada beberapa hal yang tidak jelas bagi saya: Berapa banyak istilah harmonik yang harus dimasukkan dalam regresi? Dan periode apa yang harus saya ambil? Bagaimana regresi dapat dilakukan untuk memastikan bahwa rata-rata mingguan dipertahankan (karena saya tidak ingin kecocokan harmonis yang tepat dengan data)?

Dengan menggunakan pendekatan regresi (yang juga dijelaskan dalam makalah ini ), saya berhasil mendapatkan kecocokan harmonis yang tepat dengan data ( dalam contoh saya akan berjalan melalui , jadi saya memasang 417 istilah). Bagaimana pendekatan ini dapat dimodifikasi - jika mungkin - untuk mencapai konservasi cara mingguan? Mungkin dengan menerapkan faktor koreksi untuk setiap istilah regresi?j1,...,417  

Plot harmonik yang pas adalah:

Fit harmonik yang tepat

EDIT

Menggunakan paket sinyal dan interp1fungsinya, inilah yang saya berhasil lakukan dengan menggunakan contoh data yang diatur dari bawah (terima kasih banyak kepada @noumenal). Saya menggunakan q=7karena kami memiliki data mingguan:

# 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

Ada dua masalah di sini:

  1. Kurva sepertinya pas "terlalu bagus": itu melewati setiap titik
  2. Cara mingguan tidak dilestarikan

Dataset contoh

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
sumber
2
Sastra yang lebih tua sering senang dengan formula kompleksitas rococo untuk hal semacam ini. Dalam praktiknya saya akan condong ke belakang untuk menganggap ini sebagai masalah regresi sehingga nilai-nilai yang diinterpolasi hanyalah nilai-nilai prediksi dari regresi. Lihat misalnya stats.stackexchange.com/questions/60500/... Prinsip utamanya adalah siklus utama berulang sekali per tahun.
Nick Cox
1
Dalam praktiknya Anda menginginkan kecocokan yang sangat dekat dengan data karena Anda ingin perataan lokal. Itu mungkin membutuhkan beberapa pasangan Fourier, tetapi pengembalian yang menurun ditetapkan dengan sangat cepat, sehingga setiap pasangan istilah baru segera menambahkan sangat sedikit. Anda hanya perlu menyedotnya dan melihatnya. Merencanakan segalanya membuatnya lebih jelas.
Nick Cox
1
Saya mencoba ini secara singkat dengan data Anda (menggunakan Stata, bukan R). Secara singkat, meskipun ada tanda musiman pada data Anda, pendekatan ini tidak cukup teratur untuk bekerja sama sekali: misalnya, tidak hanya waktu puncaknya sangat bervariasi, tetapi juga jumlah kasus pada puncaknya; dalam beberapa tahun tetapi tidak semua ada puncak sekunder di akhir tahun kalender. Selanjutnya, musim ini diperparah dengan tren jangka panjang yang ditandai. Dugaan saya adalah bahwa untuk mendapatkan case harian Anda harus menggunakan interpolasi lokal yang ketat atau memperlancar seri mingguan yang diperluas tujuh kali lipat.
Nick Cox
1
Dalam rekayasa sistem kontrol kriteria Nyquist digunakan sebagai dasar untuk laju pengambilan sampel. Dikatakan sampel lebih besar dari dua kali frekuensi tertinggi dalam data Anda. Dalam praktiknya, lebih konvensional untuk sampel di atas 5x frekuensi tertinggi yang ingin Anda selesaikan. Jika input Anda adalah data mingguan, maka Nyquist menyarankan frekuensi tertinggi yang dapat diselesaikan dalam urutan 2 minggu. Mungkin lebih baik jika Anda memiliki statistik mingguan lainnya untuk menginformasikan pengambilan sampel dan menopang mean. en.wikipedia.org/wiki/Nyquist-Shannon_sampling_theorem
EngrStudent
+1 Pertanyaan bagus! Apakah Anda tahu bagaimana mendeteksi sinyal dalam noise (lebih disukai di R), asalkan Anda memiliki sejumlah distribusi, di mana noise adalah Gaussian dan sinyal adalah Gaussian plus bagian lain? Saya melihat banyak paket dan fungsi (sinyal, fft (), dll) dan bahkan bermain dengan data, mencoba menerapkan transformasi Fourier dan bahkan langkah-langkah entropi, tetapi sejauh ini tidak berhasil. Saya mencoba menjawab pertanyaan (bukan milik saya) dan belajar sesuatu yang baru di sepanjang jalan, karena saya menemukan topik yang cukup menarik.
Aleksandr Blekh

Jawaban:

4

Saya bukan ahli transformasi Fourier, tapi ...

Total rentang sampel Epstein adalah 24 bulan dengan laju sampel bulanan: 1/12 tahun. Rentang sampel Anda adalah 835 minggu. Jika tujuan Anda adalah memperkirakan rata-rata selama satu tahun dengan data dari ~ 16 tahun berdasarkan data harian, Anda memerlukan tingkat sampel 1/365 tahun. Jadi gantilah 52 untuk 12, tetapi pertama-tama distandarisasi unit dan perluas 835 minggu Anda menjadi 835 * 7 = 5845 hari. Namun, jika Anda hanya memiliki titik data mingguan, saya menyarankan laju sampel 52 dengan kedalaman 16 atau 17 untuk analisis puncak, atau 32 atau 33 untuk perbandingan genap / ganjil. Jadi opsi input default meliputi: 1) untuk menggunakan cara mingguan (atau median deviasi absolut, MAD, atau sesuatu sejauh itu) atau 2) untuk menggunakan nilai harian, yang memberikan resolusi lebih tinggi.

Liebman et al. memilih titik cut-off jmax = 2. Oleh karena itu, Gambar 3. mengandung lebih sedikit parsial dan dengan demikian lebih simetris di bagian atas sinus dibandingkan dengan Gambar 2. (Parsial tunggal pada frekuensi dasar akan menghasilkan gelombang sinus murni. ) Jika Epstein akan memilih resolusi yang lebih tinggi (mis. Jmax = 12) transformasi mungkin hanya akan menghasilkan fluktuasi kecil dengan komponen tambahan, atau mungkin ia tidak memiliki kekuatan komputasi.

Melalui inspeksi visual terhadap data Anda, Anda tampaknya memiliki 16-17 puncak. Saya sarankan Anda mengatur jmax atau "bit depth" ke 6, 11, 16 atau 17 (lihat gambar) dan bandingkan hasilnya. Semakin tinggi puncak, semakin banyak kontribusinya pada bentuk gelombang kompleks asli. Jadi dengan asumsi resolusi 17-band atau kedalaman bit parsial ke-17 berkontribusi minimal terhadap pola gelombang asli dibandingkan dengan puncak ke-6. Namun, dengan resolusi 34 band Anda akan mendeteksi perbedaan antara genap dan puncak aneh seperti yang disarankan oleh lembah yang cukup konstan. Kedalaman bit tergantung pada pertanyaan penelitian Anda, apakah Anda tertarik pada puncak saja atau pada puncak dan lembah, tetapi juga seberapa tepatnya Anda ingin mendekati seri aslinya.

Analisis Fourier mengurangi poin data Anda. Jika Anda membalikkan fungsi pada kedalaman bit tertentu menggunakan transformasi Fourier, Anda mungkin bisa memeriksa silang apakah estimasi rata-rata baru sesuai dengan cara asli Anda. Jadi, untuk menjawab pertanyaan keempat Anda: parameter regresi yang Anda sebutkan tergantung pada sensitivitas dan resolusi yang Anda butuhkan. Jika Anda tidak menginginkan kecocokan yang tepat, maka dengan segala cara, cukup masukkan sarana mingguan dalam transformasi. Namun, berhati-hatilah karena kedalaman bit yang lebih rendah juga mengurangi data. Sebagai contoh, perhatikan bagaimana overlay harmonik Epstein pada analisis Lieberman dan rekan melewatkan titik tengah dari fungsi langkah, dengan kurva miring sedikit ke kanan (mis. Temp. Est. Terlalu tinggi), pada bulan Desember pada Gambar 3.


Parameter Liebman dan Kolega:

  • Kedalaman Bit: 2

Parameter Epstein:

  • Tingkat Sampel: 12 [setiap bulan]
  • Rentang Sampel: 24 bulan
  • Kedalaman Bit: 6

Parameter Anda:

  • Tingkat Sampel: 365 [setiap hari]

  • Rentang Sampel: 5845 hari

Pendekatan Kedalaman Bit Tepat

Cocok tepat berdasarkan inspeksi visual. (Jika Anda memiliki kekuatan, lihat saja apa yang terjadi dibandingkan dengan kedalaman bit yang lebih rendah.)

  • Spektrum Penuh (puncak): 17
  • Spektrum Penuh (datar / ganjil): 34

Pendekatan Kedalaman Bit Variabel

Ini mungkin yang ingin Anda lakukan:

  • Bandingkan Hanya Puncak: 6, 11, 16, 17
  • Bandingkan Even / Odd: 12, 22, 32, 34
  • Mensintesis ulang dan membandingkan cara

Pendekatan ini akan menghasilkan sesuatu yang mirip dengan perbandingan Angka di Epstein jika Anda membalikkan transformasi lagi, yaitu mensintesis parsial ke dalam perkiraan dari deret waktu asli. Anda juga bisa membandingkan titik-titik diskrit dari kurva yang disintesis ulang dengan nilai rata-rata, bahkan mungkin menguji perbedaan signifikan untuk menunjukkan sensitivitas pilihan kedalaman bit Anda.


PEMBARUAN 1:

Kedalaman Bit

Bit - kependekan dari digit biner - adalah 0 atau 1. Bit 010101 akan menggambarkan gelombang persegi. Kedalaman bit adalah 1 bit. Untuk menggambarkan gelombang gergaji, Anda akan membutuhkan lebih banyak bit: 0123210. Semakin kompleks gelombang, semakin banyak bit yang Anda butuhkan:

Ini adalah penjelasan yang agak disederhanakan, tetapi semakin kompleks deret waktunya, semakin banyak bit yang diperlukan untuk memodelkannya. Sebenarnya, "1" adalah komponen gelombang sinus dan bukan gelombang persegi (gelombang persegi lebih seperti 3 2 1 0 - lihat gambar terlampir). 0 bit akan menjadi garis datar. Informasi hilang dengan pengurangan kedalaman bit. Misalnya, audio berkualitas CD biasanya 16 bit, tetapi audio kualitas telepon darat sering sekitar 8 bit.

Harap baca gambar ini dari kiri ke kanan, dengan fokus pada grafik:

FFT

Anda sebenarnya baru saja menyelesaikan analisis spektrum daya (meskipun pada resolusi tinggi dalam gambar Anda). Tujuan Anda berikutnya adalah untuk mencari tahu: Berapa banyak komponen yang saya butuhkan dalam spektrum daya untuk secara akurat menangkap alat-alat deret waktu?

PEMBARUAN 2


Untuk Menyaring atau tidak Menyaring

Saya tidak sepenuhnya yakin bagaimana Anda akan memperkenalkan kendala dalam regresi karena saya hanya akrab dengan batasan interval, tetapi mungkin DSP adalah solusi Anda. Inilah yang saya pikir sejauh ini:

  • Langkah 1. Memecah seri menjadi komponen sinus melalui fungsi Fourier pada set data lengkap (dalam beberapa hari)

  • Langkah 2. Menciptakan deret waktu melalui transformasi Fourier terbalik, dengan batasan-rata tambahan yang digabungkan dengan data asli: penyimpangan interpolasi dari cara asli harus membatalkan satu sama lain (Harzallah, 1995).

Tebakan terbaik saya adalah Anda harus memperkenalkan autoregresi jika saya memahami Harzallah (1995, Gambar 2) dengan benar. Jadi itu mungkin sesuai dengan filter tanggapan tak terbatas (IIR)?

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

Singkatnya:

  1. Berasal berarti dari data mentah
  2. Fourier Transform data mentah
  3. Fourier Inverse Transform mentransformasikan data.
  4. Saring hasilnya menggunakan IIR

Mungkin Anda bisa menggunakan filter IIR tanpa melalui analisis Fourier? Satu-satunya keuntungan dari analisis Fourier seperti yang saya lihat adalah mengisolasi dan menentukan pola mana yang berpengaruh dan seberapa sering mereka terulang kembali (yaitu berosilasi). Anda kemudian dapat memutuskan untuk menyaring yang memberikan kontribusi lebih sedikit, misalnya menggunakan filter takik sempit di puncak kontribusi terendah (atau filter berdasarkan kriteria Anda sendiri). Sebagai permulaan, Anda bisa memfilter lembah aneh yang kurang berkontribusi yang lebih mirip noise di "sinyal". Kebisingan ditandai oleh sangat sedikit kasus dan tidak ada pola. Filter sisir pada komponen frekuensi aneh dapat mengurangi kebisingan - kecuali jika Anda menemukan pola di sana.

Inilah beberapa binning sewenang-wenang — hanya untuk tujuan penjelasan: Bisakah Anda melihat suara di lembah?

Ups - Ada Fungsi R untuk itu !?

Ketika mencari filter IIR saya kebetulan menemukan fungsi R interpolasi dalam paket Sinyal. Lupakan semua yang saya katakan sampai saat ini. Interpolasi harus bekerja seperti Harzallah: http://cran.r-project.org/web/packages/signal/signal.pdf

Bermain-main dengan fungsinya. Sebaiknya lakukan triknya.


PEMBARUAN 3

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

Set xi ke cara mingguan asli.

noumenal
sumber
Terima kasih banyak atas jawaban ini! Tujuan penelitian saya sederhana: Saya memiliki rata-rata mingguan dan ingin mendapatkan taksiran harian dan rata-rata taksiran harian (diinterpolasikan) selama satu minggu harus sama dengan rata-rata mingguan (yaitu titik data asli). Apakah Anda pikir itu mungkin? Lebih lanjut, saya tidak mengerti apa arti "kedalaman bit" dan "analisis puncak" (saya tidak punya pengalaman dengan transformasi Fourier apa pun).
COOLSerdash
1
@COOLSerdash lihat pembaruan saya. Ya itu mungkin! Tetapi Anda perlu mencari cara terbaik untuk membandingkan perkiraan rata-rata dari deret waktu yang disintesis ulang dengan sarana asli dari deret waktu asli.
noumenal
(Btw: +1 besok, saya tidak bisa lagi memilih hari ini). Terima kasih banyak untuk pembaruannya, sekarang lebih jelas. Saya memikirkan prosedur berikut ini: 1) cocokkan fungsi fourier dengan sarana mingguan dengan regresi, 2) gunakan prediksi regresi untuk "mengisi" kesenjangan antara nilai-nilai mingguan (yaitu untuk mendapatkan nilai harian) 3) untuk setiap minggu, hitung rata-rata semua nilai harian dan nilai rata-rata ini harus sama dengan nilai aslinya. Di koran, Epstein menggunakan semacam faktor "koreksi" untuk memaksa fungsi memiliki sifat yang diinginkan, tapi saya masih tidak yakin bagaimana melakukannya dengan regresi.
COOLSerdash
@COOLSerdash Lihat pembaruan 2! Lewati ke paragraf terakhir.
noumenal
Benar-benar fantastis! Terima kasih banyak atas penelitiannya. Perhatikan bahwa saya telah berhasil menerapkan pendekatan Harzallah menggunakan splines (linier dan kubik). Jadi saya kira saya perlu interp. Saya telah mengedit pertanyaan saya. Sekali lagi terima kasih banyak.
COOLSerdash