שְׁאֵלָה:
איתור תקופות של סדרת זמן כללית
gianluca
2010-08-04 05:32:13 UTC
view on stackexchange narkive permalink

פוסט זה הוא המשך של פוסט נוסף הקשור ל שיטה כללית לזיהוי חריגים בסדרות זמן. בעיקרון, בשלב זה אני מעוניין בדרך חזקה לגלות את המחזוריות / עונתיות של סדרת זמן גנרית המושפעת מרעש רב. מנקודת מבט של מפתח, אני רוצה ממשק פשוט כגון:

int-inte_period לא חתום (vector<double> v);

כאשר v הוא המערך המכיל את הדגימות, וערך ההחזרה הוא תקופת האות. הנקודה העיקרית היא ששוב, אני לא יכול להניח שום הנחה לגבי ניתחתי את האות. כבר ניסיתי גישה המבוססת על התאמת האות האוטומטית (איתור פסגות המתאם), אך היא אינה חזקה כפי שהייתי רוצה.

ניסית xts :: מחזוריות?
שבע תשובות:
#1
+53
Rob Hyndman
2010-08-04 10:41:03 UTC
view on stackexchange narkive permalink

אם באמת אין לך מושג מה המחזוריות, כנראה שהגישה הטובה ביותר היא למצוא את התדר המתאים למקסימום של צפיפות הספקטרום. עם זאת, הספקטרום בתדרים נמוכים יושפע מהמגמה, לכן עליך לדחות את הסדרה תחילה. פונקציית R הבאה צריכה לעשות את העבודה ברוב הסדרות. זה רחוק מלהיות מושלם, אבל בדקתי את זה על כמה עשרות דוגמאות ונראה שזה עובד בסדר. הוא יחזיר 1 עבור נתונים שאין להם מחזוריות חזקה, ואורך תקופה אחרת.

עדכון: גרסה 2 של הפונקציה. זה הרבה יותר מהיר ונראה חזק יותר.

  find.freq <- function (x) {n <- length (x) spec <- spec.ar (c (x), plot = FALSE) if (max (spec $ spec) >10) # סף שרירותי שנבחר על ידי ניסוי וטעייה. {period <- round (1 / spec $ freq [which.max (spec $ spec)]) if (period == Inf) # מצא את המקסימום המקומי הבא {j <- אשר (diff (spec $ spec) >0) אם ( אורך (j) >0) {nextmax <- j [1] + which.max (spec $ spec [j [1]: 500]) תקופה <- עגול (1 / spec $ freq [nextmax])} תקופה אחרת <- 1}} תקופה אחרת <- 1 החזר (נקודה)}  
תודה. שוב אנסה גישה זו בהקדם האפשרי ואכתוב כאן את התוצאות הסופיות.
הרעיון שלך די טוב, אבל במקרה שלי, הוא לא מצליח לזהות את המחזוריות של סדרת זמן ממש (ולא כל כך רועשת) כמו http://dl.dropbox.com/u/540394/chart.png. לפי הגישה ה"אמפירית "שלי (המבוססת על מתאם אוטומטי), האלגוריתם הפשוט שכתבתי מחזיר תקופה מדויקת של 1008 (בעל מדגם כל עשר דקות, זה אומר 1008/24/6 = 7, אז מחזוריות שבועית). הבעיות העיקריות שלי הן: 1) זה איטי מדי להתכנס (זה דורש הרבה נתונים היסטוריים) ואני זקוק לגישה תגובתית ומקוונת; 2) זה לא יעיל כמו לעזאזל מנקודת מבט של שימוש בזיכרון; 3) זה לא חזק את כל;
תודה. למרבה הצער, זה עדיין לא עובד כמו שהייתי מצפה. באותה סדרת זמן של ההערה הקודמת היא מחזירה את 166, וזה צודק רק באופן חלקי (מנקודת מבטי, התקופה השבועית הניכרת מעניינת יותר). ובאמצעות סדרת זמן רועשת מאוד, כמו זו http://dl.dropbox.com/u/540394/chart2.png (ניתוח חלון מקלט TCP), הפונקציה מחזירה 10, בעוד שהייתי מצפה ל -1 (אני יכול ' לא רואה כל מחזוריות ברורה). BTW אני יודע שיהיה ממש קשה למצוא את מה שאני מחפש, מכיוון שאני מתמודד עם אותות שונים מדי.
166 אינו אומדן גרוע של 168. אם ידוע לך שהנתונים נצפים מדי שעה עם תבנית שבועית, אז מדוע בכלל להעריך את התדירות?
כי אני צריך לנתח הרבה סדרות זמן (נניח 100 מדדי רשת), ורק לחלקם יש מחזוריות שבועית. בכל מקרה, אני מניח שביישומי אשתמש באלגוריתם הדומה לפונקציה שלך, ואבחין ידנית בין המחזוריות השבועית. באמת תודה על התמיכה שלך, אני מאוד מעריך (וממשיך בעבודה הטובה עם ספריית התחזיות :-))
בדקתי פונקציה זו בדוגמה פשוטה: x = c (58.89446, 37.31097, 53.99865, 26.13904, 34.74298) ו- y = ts (rep_len (x, 15 * אורך (x)). באמצעות ההגדרה לעיל ציפיתי לקבל15 כמו find.freq (y) (או משהו קרוב), אבל אני מקבל 3. מה חסר לי כאן?
מדוע לא לכלול זאת באיזו חבילה?יש הרבה משימות כאשר המחזוריות אינה ידועה ..
גרסה משופרת נמצאת בחבילת התחזית כ- 'findfrequency'
#2
+10
Rich
2010-08-10 23:41:11 UTC
view on stackexchange narkive permalink

אם אתה מצפה שהתהליך יהיה נייח - המחזוריות / עונתיות לא תשתנה עם הזמן - אז משהו כמו מחזור צ 'ריבועי (ראה למשל סוקולובה ובושל, 1978) עשוי להיות בחירה טובה. הוא משמש בדרך כלל לניתוח נתונים ביממה שיכולים להכיל כמויות גדולות מאוד של רעש, אך הוא צפוי להיות בעל מחזוריות יציבה מאוד.

גישה זו אינה מניחה את צורת צורת הגל (מלבד זאת. הוא עקבי ממחזור למחזור), אך דורש שכל רעש יהיה ממוצע קבוע ולא מתואם לאות.

  chisq.pd <- function (x, min. period, max. נקודה, אלפא) {N <- אורך (x) שונות = NULL תקופות = seq (מינימום תקופה, מקסימום תקופה) רשימת שורה = NULLfor (lc בתקופות) {ncol = lc nrow = רצפה (N / ncol) רשימת שורה = c ( רשימת שורות, מס ') x.trunc = x [1: (ncol * nrow)] x.reshape = t (array (x.trunc, c (ncol, nrow))) שונות = c (שונות, var (colMeans (x. צורה מחדש))}} Qp = (רשימת שורה * תקופות * שונות) / var (x) df = תקופות - 1pvals = 1-pchisq (Qp, df) לעבור. תקופות = תקופות [pvals<alpha] pass.pvals = pvals [pvals<alpha] # return (cbind (pass.periods, pass.pvals)) return (cbind (נקודות [pvals == min (pvals)], pvals [pvals == min (pvals)]))} x = cos ((2 * pi / 37) * (1: 1000)) + rnorm (1000) chisq.pd (x, 2, 72, .05)  

שתי השורות האחרונות הן רק דוגמה המראה שהוא יכול לזהות את התקופה של פונקציה טריגונומטרית טהורה, אפילו עם הרבה רעש תוסף.

כפי שנכתב, הטיעון האחרון ( אלפא ) בשיחה מיותר, הפונקציה פשוט מחזירה את התקופה 'הטובה ביותר' שהיא יכולה למצוא; בטל את ההצהרה על הצהרת ההחזרה והגב את השנייה כדי שתחזיר רשימה של כל התקופות המשמעותיות ברמה alpha .

פונקציה זו אינה מבצעת שום סוג של בדיקת שפיות כדי לוודא שהכנסת תקופות שניתן לזהות, וגם אינה (יכולה) לעבוד עם תקופות חלקיות, ואין שום סוג של בקרת השוואה מרובה מובנית אם אתה מחליט להסתכל על מספר תקופות. אך מלבד זאת הוא אמור להיות חזק יחסית.

נראה מעניין אבל אני לא מבין את התפוקה, זה לא אומר לי איפה מתחילה התקופה ורוב הערכים של 1.
#3
+4
Wesley Burr
2010-08-06 07:48:10 UTC
view on stackexchange narkive permalink

כדאי להגדיר את מה שאתה רוצה בצורה ברורה יותר (לעצמך, אם לא כאן). אם מה שאתה מחפש הוא התקופה הנייחת הכי משמעותית סטטיסטית הכלולה בנתונים הרועשים שלך, יש למעשה שני מסלולים שצריך לקחת:

1) חישב הערכת אוטוקורלציה חזקה, וקח את המקדם המרבי
2) חישב אומדן צפיפות ספקטרלי כוח חזק, וקח את המקסימום של הספקטרום

הבעיה עם מספר 2 היא שבכל סדרת זמן רועשת, תקבל כמות גדולה של כוח בתדרים נמוכים, מה שהופך קשה להבחין. ישנן כמה טכניקות לפתרון בעיה זו (כלומר להלבין מראש, ואז להעריך את ה- PSD), אך אם התקופה האמיתית מהנתונים שלך ארוכה מספיק, הגילוי האוטומטי יהיה בטוח.

ההימור הטוב ביותר שלך הוא כנראה ליישום שגרת אוטוקורלציה חזקה כמו ניתן למצוא בפרקים 8.6, 8.7 ב סטטיסטיקה חזקה - תיאוריה ושיטות מאת מרונה, מרטין ויוחאי. חיפוש בגוגל "דורבין-לוינסון חזק" יניב גם כמה תוצאות.

אם אתה רק מחפש תשובה פשוטה, אני לא בטוח שקיימת תשובה. גילוי תקופות בסדרות זמן יכול להיות מסובך, ולבקש שגרה אוטומטית שיכולה לבצע קסמים עשויה להיות יותר מדי.

תודה על המידע היקר שלך, אני אסתכל על הספר הזה בוודאות.
#4
+4
babelproofreader
2010-08-10 22:29:28 UTC
view on stackexchange narkive permalink

אתה יכול להשתמש בהילברט טרנספורמציה מתיאוריית DSP כדי למדוד את התדירות המיידית של הנתונים שלך. באתר http://ta-lib.org/ יש קוד קוד פתוח למדידת תקופת המחזור הדומיננטית של נתונים פיננסיים; הפונקציה הרלוונטית נקראת HT_DCPERIOD; ייתכן שתוכל להשתמש בזה או להתאים את הקוד למטרותיך.

#5
+3
Fabrizio Maccallini
2016-12-29 19:15:47 UTC
view on stackexchange narkive permalink

גישה אחרת יכולה להיות פירוק מצב אמפירי.חבילת R נקראת EMD שפותחה על ידי ממציא השיטה:

  require (EMD) ndata <- 3000 tt2 <- seq (0, 9, length =ndata) xt2 <- sin (pi * tt2) + sin (2 * pi * tt2) + sin (6 * pi * tt2) + 0.5 * tt2 נסה <- emd (xt2, tt2, גבול = "גל") ### התוויית ה- par של ה- IMF (mfrow = c (נסה $ nimf + 1, 1), mar = c (2,1,2,1)) rangeimf <- range (נסה $ imf) עבור (i in 1: נסה $ nimf) {עלילה (tt2, נסה $ imf [, i], type = "l", xlab = "", ylab = "", ylim = rangeimf, main = הדבק (i, "-th IMF", sep = ""));abline (h = 0)} מגרש (tt2, נסה $ residue, xlab = "", ylab = "", main = "residue", type = "l", axes = FALSE);box ()  

השיטה סומנה כ'אמפירית 'מסיבה טובה ויש סיכון שהפונקציות Intrinsic Mode (רכיבי התוסף הבודדים) יתערבבו.מצד שני השיטה מאוד אינטואיטיבית ועשויה להועיל לבדיקה חזותית מהירה של מחזוריות.

#6
  0
Chris
2015-05-02 20:24:14 UTC
view on stackexchange narkive permalink

בהתייחס לפוסט של רוב הינדמן לעיל https://stats.stackexchange.com/a/1214/70282

פונקציית find.freq עובדת בצורה מעולה. במערכת הנתונים היומית שבה אני משתמש, זה הסתדר כראוי שהתדירות תהיה 7.

כשניסיתי את זה רק בימי השבוע, הוא הזכיר שהתדר הוא 23, שהוא קרוב להפליא ל 21.42857 = 29.6 * 5/7 שהוא המספר הממוצע של ימי עבודה בחודש. (או להיפך 23 * 7/5 הוא 32.)

במבט לאחור על הנתונים היומיים שלי, התנסיתי בתחושה של נטילת התקופה הראשונה, ממוצע לפי זה ואז מציאת התקופה הבאה וכו '. להלן:

 find.freq.all = function (x) {f = find.freq (x); freqs = c (f); בעוד (f> 1) {start = 1; # נסה גם להתחיל = f; x = period.apply (x, seq (התחלה, אורך (x), f), ממוצע); f = find.freq (x); freqs = c (freqs, f); } אם (אורך (freqs) == 1) {return (freqs); } עבור (i ב- 2: אורך (freqs)) {freqs [i] = freqs [i] * freqs [i-1]; } freqs [1: (length (freqs) -1)];} find.freq.all (dailyts) #using data daily 

האמור לעיל נותן (7,28) או (7,35) תלוי מופעל אם ה- seq מתחיל ב- 1 או f. (ראה הערה למעלה.)

מה שמשמעותו שהתקופות העונתיות עבור msts (...) צריכות להיות (7,28) או (7,35).

ההיגיון נראה רגיש לתנאים ראשוניים בהתחשב ברגישות הפרמטרים של האלגוריתם. הממוצע של 28 ו -35 הוא 31.5 שזה קרוב לאורך הממוצע של חודש.

אני חושד שהמצאתי מחדש את הגלגל, מה שמו של האלגוריתם הזה? האם יש יישום טוב יותר ב- R איפשהו?

מאוחר יותר, הפעלתי את הקוד לעיל בניסיון כל ההתחלות של 1 עד 7 וקיבלתי 35,35,28,28,28,28,28 לשנייה פרק זמן. הממוצע מסתדר ל -30 שזה מספר הימים הממוצע בחודש. מעניין ...

יש מחשבות או הערות?

#7
  0
ali
2016-09-27 17:10:59 UTC
view on stackexchange narkive permalink

אפשר גם להשתמש במבחן Ljung-Box כדי להבין איזה הבדל עונתי מגיע למצב הכי טוב. עבדתי בנושא אחר והשתמשתי בזה למעשה לאותן מטרות. נסה תקופות שונות כגון 3 עד 24 לקבלת נתונים חודשיים. ובדוק כל אחד מהם על ידי ליונג-בוקס ואחסן תוצאות צ'י-סקוור. ובחר את התקופה עם הערך הנמוך ביותר של כיכר הצ'י.

הנה קוד פשוט לעשות זאת.

  minval0 <- 5000 # הקצה מספר גדול כדי להיות בטוח שערכי Chi הם קטנים minindex0 <- 0periyot <- 0for (i ב 3:24) {# מצא תקופה אופטימלית על ידי Qtests על פני נתונים מקוריים d0D1 <- diff (a, lag = i) #store results Qtest_d0D1 [[i]] <- Box.test (d0D1, lag = 20, type = "Ljung-Box") סטטיסטיקה צ'י-ריבוע # חנות sira0 [i] <- Qtest_d0D1 [[i]] [1]} # רשימת סיבובים למסגרת נתונים, ואז matrixdatam0 <- data.frame (מטריצה (unlist (sira0), nrow = length (Qtest_d0D1) -2, byrow = T)) datamtrx0 <- as.matrix (datam0 []) # get min value's indexminindex0 <- which (datamtrx0 == min (datamtrx0), arr. ind = F) periyot <- minindex0 + 2  


שאלה ותשובה זו תורגמה אוטומטית מהשפה האנגלית.התוכן המקורי זמין ב- stackexchange, ואנו מודים לו על רישיון cc by-sa 2.0 עליו הוא מופץ.
Loading...