« Recension av Keith Devlin, Gary Lorden: The Numbers Behind Numb3rs: Solving Crime with Mathematics | Main | Bloggträff i Malmö torsdagen 4 oktober 18:30 på Kin Long »

augusti 21, 2007

Upptäcka kraftiga förändringar i tidsserier (changepoint detection)

I Recension av Keith Devlin, Gary Lorden: The Numbers Behind Numb3rs: Solving Crime with Mathematics nämndes att kapitel 4 When does the Writing First Appear on the Wall? - Changepoint Detection var det som gav mig mest ny information. Så naturligtvis har jag lekt lite med detta.

I Numb3rs-avsnittet Hardball (tredje säsongen) nämner Charlie Epps (matematikern i Numb3rs) Shiryayev-Roberts metod i samband med baseballresultat för att upptäcka förändringar i tidsserier, men den är "slightly more technical to describe" än den metod som E.S. Page skapade varpå man beskriver Pages. (Kapitlet är troligen skrivet av Gary Lorden eftersom han tidigare skrivit om changepoint detection. Och troligen är det E.S. Page paper "A test for a change in a parameter occurring at an unknown point" som avses. Men ingen av detta står uttryckligen i boken. Gnäll^2.)

Först presenteras kortfattat Pages metod, följt av min R-implementation. Hur väl den följer Page har jag faktiskt inte lyckats röna ut. Här används emellanåt "förändringspunkt" i stället för "changepoint"; "breakpoint" används också nedan.


Introduktion och beskrivning av Pages metod
En specifik händelse inträffar eller inte under en dag. Om den inträffar markeras detta i en tidsserie med 1, annars 0.

Bokens exempel är en händelse som inträffat c:a 1 gång per månad, dvs sannolikheten är 1/30. Sannolikheten för att händelsen inte inträffar är således 1-1/30 = 29/30. Detta är normalfallet. För att citera bokens lättsamma exempel på några sådana händelser som kan tänkas inträffa (i genomsnitt) en gång i månaden:


Examples abound - a New Yorker finds a parking place on the street in front of her apartment, a husband actually offers to take out the garbage, a local TV news show doesn't lead off with a natural disaster or violent crime, and so on.

Mer seriöst (och det beskrivs mer sådant i kapitlet) skulle man kunna tänka sig någon form av olycka eller sjukdomsfall som rapporteras varje dag. Notera att vi här endast arbetar med om händelsen inträffar eller inte; det är inte frågan om antalet olyckor etc. (Se i slutet av denna anteckning för hur man arbetar med numeriska värden.)

Vi vill nu upptäcka ("detektera") om/när det blir en drastisk förändring i hur ofta dessa händelser inträffar, t.ex. att de plötsligt inträffar i genomsnitt en gång i veckan i stället, dvs med sannolikheten 1/7. Låt oss kalla detta för extremfallet. Sannolikheten för att en extremhändelse inte inträffar en gång i veckan är således 1 - 1/7 = 6/7.

Falsklarm: Problemet med händelser som inträffar i genomsnitt en gång i månaden är att två händelser som i genomsnitt händer en gång i månaden faktiskt kan inträffa två intilliggande dagar eller klustra sig på andra sätt. Vi vill undvika falsklarm när slumpen gör det den ska (nämligen att slumpa).

Ett metod att minska antalet sådana falsklarm är alltså E.S. Pages metod som boken går igenom med ett exempel. Här har jag tagit beskrivningen från boken rakt av.

Det finns ett numeriskt index S som "spårar" aktiviteten dag för dag.
* Initialt är S = 1.
* För varje dag uppdateras värdet enligt följande:

* Till slut finns ett gränsvärde för när vi ska varna för att en s.k. changepoint har inträffat. Om S överstiger detta gränsvärde anses händelserna nu komma en gång i veckan (extremfallet). I boken har man satt 50 som detta värde. (Det förs en liten diskussion hur olika gränsvärden påverkar hur väl man kan upptäcka förändringarna, men det går jag inte in på här.)


Bokens exempel
Boken använder följande händelsesekvens som exempel, dvs de första två dagarna inträffade händelsen inte, på tredje dagen inträffade den osv.

0 0 1 0 0 0 0 0 1

Från denna serie får man motsvarande S-värden (tänk på att S återställs till 1 om S < 1 och att S är 1 från början). De första två dagarna blir 1 enligt:


S * 0.8867 = 1 * 0.8867 = 0.8867 < 1, varpå S återställs till 1

Här är en liten tabell för att visa utvecklingen av S-värdet.


H S-värde
0 1
0 1
1 4.285714
0 3.800141
0 3.369583
0 2.987808
0 2.649287
0 2.349122
0 2.082965
1 8.926994

Inget av S-värdena är större än gränsvärdet (50), så ingen varning utfärdas för dessa första 10 dagar. Om händelsen inträffar även nästa dag (dag 11) skulle S-värdet vara 38.25854 (8.926994*4.285714) vilket inte heller ger någon varning. Om den även inträffar påföljande dag (dag 12) skulle det en varning utfärdas eftersom 8.926994*4.285714*4.285714 = 163.9652 >= 50.

Poängen med Pages S-serie är att den stiger kraftigt om händelsen inträffar, men avtar långsamt om händelsen inte inträffar. (Eller snarare: i detta exempel är det "kraftigt", det exakta beteendet beror på de specifika sannolikheterna.)


R-funktionen changepoint
Pages metod har implementerat i statistikpaketet R med följande funktion:


changepoint <- function(x,normal=1/30,non.normal=1/7,detect.value=50,plot=T){
occurs <- non.normal/normal # händelsen inträffar
not.occurs <- (1-non.normal)/(1-normal) # händelsen inträffar inte
S <- 1
vec = NULL # vektor med S-värdena
len <- length(x)
alerted <- FALSE;
cp <- NULL; # changepoint
S.cp <- 0; # värdet för S i changepoint
for (i in 1:len) {
obs <- x[i] # observationen
# justera S enligt metoden:
S <- ifelse(obs, S*occurs, S*not.occurs)
S <- ifelse(S<1,1,S) # reset om mindre än 1

# varna om S-värdet överstiger gränsvärdet
if (plot && !alerted && S >= detect.value) {
cat("Alert! observation #", i, ": ", S, " >= ", detect.value,"\n")
cp <- i; S.cp <- S; alerted <- TRUE
}
vec <- c(vec,S)
}
csum <- cumsum(x) # kummulativa summan av händelseserien
cumsum.val <- 0 # antal observationer
if (!is.null(cp) && cp > 0) { cumsum.val <- csum[cp] }
m<-max(csum) # för grafen
vec.len <- length(vec)
# plotta tre grafer
if(plot) {
op <- par(mfrow=c(3,1)) # plotta grafer
plot(x,type="l",
main=paste("Changepoint value: ", cp, " S-value:", round(S.cp,2) ))
if (!is.null(cp) && cp>0) { lines( rep( cp, 2), 0:1, col="red") }
plot(csum,type="l", main=paste("Cumulative sum. Num obs.:",cumsum.val) )
if (!is.null(cp) && cp>0) { lines( rep( cp, m+1), 0:m, col="red")}
plot(vec,type="l",main=paste("S values. Detect value: ", detect.value) )
abline( rep(detect.value, vec.len), 0:vec.len,col="red")
par(op)
}

# returvärde
list(vec=vec,cp=cp,cumsum.val=cumsum.val);
}


Funktionen har följande argument, där defaultvärdena är bokens exempel.

x: en vektor av respektive dagars värden om händelsen inträffade eller inte (1 eller 0)
normal:. Sannolikheten för att normalfallet inträffar, default 1/30 för att följa boken
non.normal: Sannolikheten för att extremfallet inträffar, default 1/7
detect.value: gränsvärde för changepoint-varning (default 50)
plot: huruvida man ska plotta respektive skriva ut en varning (när man gör flera simuleringar vill man inte att det plottas och håller på).

Följande värden returneras :

vec: vektor med S-värdena
cp: changepoint, dvs den observation där S-värdet överstiger varningsvärdet
cumsum.val: vilket är summan av antalet händelser som inträffat vid changepoint (egentligen endast för att jag ville kolla lite)


När man kör programmet visas tre grafer om man inte angivit plot=F i funktionen (visas nedan).

1) Först händelsegrafen, med 0 för ej inträffad och 1 för inträffad händelse. Den lodräta röda linjen visar changepointen om sådan finns. "S-value" är S-värdet för changepoint.
2) En graf för det kumulativa antalet inträffade händelser, samt en röd lodrät linje för changepoint. "Num. obs" är det antalet inträffade händelser för changepoint.
3) Den tredje grafen visar S-värdena. "Detect value" är gränsen för S-värdet (dvs då man ska slå larm) och visas som en röd vågrät linke. Notera att det kan vara svårt att se linjen om största S-värde är stort.


Körning av R-koden
Här är bokens exempel:

> x = c(0,0,1,0,0,0,0,0,0,1)
> changepoint(x, 1/30, 1/7, 50, T)
$vec
[1] 1.000000 1.000000 4.285714 3.800141 3.369583 2.987808 2.649287 2.349122
[9] 2.082965 8.926994

$cp
NULL

$cumsum.val
[1] 0


Som vi såg tidigare blev det alltså ingen changepoint eftersom S-värdet inte översteg 50. cp är därför NULL och antalet händelser vid changepoint är därför 0.

Graferna ser ut på följande sätt:


Ett mer intressant exempel (slumpat)
Här är ett mer intressant exempel. Det är en längre sekvens som jag helt enkelt slumpar fram med utgångspunkt från exemplets förutsättningar (normalfall resp. extremfall man vill detektera). I ärlighetens namn körde jag funktionen tills det blev en tidsserie med en changepoint.

För de första 30 dagarna inträffar händelsen (1) med sannolikheten 1/30, sedan kommer 30 dagar med sannolikheten 1/7 för 1. Här är slumpfunktionen:


c(sample(0:1,30,prob=c(29/30,1/30),rep=T),sample(0:1,30,prob=c(6/7,1/7),rep=T));

Så kör vi!


> x <-c(0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,
0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0)
> changepoint(x)
Alert! observation # 47 : 168.0480 >= 50
$vec
[1] 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
[7] 1.000000 1.000000 1.000000 1.000000 1.000000 4.285714
[13] 3.800141 3.369583 2.987808 2.649287 2.349122 2.082965
[19] 1.846964 1.637702 1.452150 1.287621 1.141732 4.893139
[25] 4.338744 3.847162 3.411277 3.024778 2.682069 2.378189
[31] 10.192239 43.681023 38.731942 34.343594 30.452447 27.002170
[37] 23.942811 21.230079 18.824700 16.691853 14.800657 13.123736
[43] 11.636810 10.318354 44.221515 39.211196 168.047983 149.008064
[49] 132.125377 566.251614 502.095027 445.207413 394.765194 1691.850832
[55] 1500.163299 1330.194059 1179.482416 1045.846478 927.351557 822.282168

$cp
[1] 47

$cumsum.val
[1] 6

Kommentar: Här ser vi alltså att det blev en varning vid observation #47 (dag 47) och då hade händelsen inträffat totalt 6 gånger. Det var på den 16 dagen (31 + 16 = 47) efter att "extremfallet" började gälla (dvs med 1/7 som sannolikhet). Man kan f.ö. notera att just i denna sekvens blev det 2 stycken händelser efter varandra precis i början på "extremfallet".

Grafen för detta exempel:


Not: Om man kör programmet på ovanstående sätt, dvs batchar hela serien på en gång, så förfelar det sitt syfte. Man ska naturligtvis köra programmet varje dag efter man fyllt i dagens "händelsevärde" och kontrollera om det kommer en varning. Programmet fortsätter alltså att köra efter varningen. I stället borde kanske lampor i lustiga färger börja blinka, sirener börja tjuta och den stora valvdörren börjar stängas så att man måste skynda sig in i valvet innan dörren stängs helt med ett tungt och djupt "swoonk". Minsann kommer där en jeep i högsta fart genom stängslet och fyra personer rusar fram till dörren. Kommer de att klara sig? Ojoj, vad spännade det är...

Andra changepoint-metoder
Pages metod (enligt bokens förklaring och exempel i alla fall) hanterar endast binära tidserier (0/1). Det finns även metoder som detekterar changepoints i tidsserier med numeriska värden. Här nämns några R-funktioner.

Ett standardexempel är datamängden Nile: (via ?Nile i R) Measurements of the annual flow of the river Nile at Ashwan 1871-1970. Det intressanta med denna tidsserie är att år 1898 minskade flödet drastiskt eftersom man byggde första Ashwan-dammen. Kan denna förändring uppptäckas med changepoint detection-metoderna (även kallad "change point analysis", "change detection" etc). Svaret är naturligtvis jakande.


R-paketet struccchange
Med functionen breakpoints() i paketet strucchange visas flera olika changepoints (benämns breakpoints och breakdates för respektive observationsnummer och årtal). Brytpunkterna sorteras på hur stora (viktiga) de är. Nedan ser man att observation 28 (år 1898) är den viktigaste förändringspunkten; det hände också något vid observation 83 (år 1953) etc.


> library(strucchange)
> data(Nile)
> summary(breakpoints(Nile~1))

Optimal (m+1)-segment partition:

Call:
breakpoints.formula(formula = Nile ~ 1)

Breakpoints at observation number:

m = 1 28
m = 2 28 83
m = 3 28 68 83
m = 4 28 45 68 83
m = 5 15 30 45 68 83

Corresponding to breakdates:

m = 1 1898
m = 2 1898 1953
m = 3 1898 1938 1953
m = 4 1898 1915 1938 1953
m = 5 1885 1900 1915 1938 1953

Fit:

m 0 1 2 3 4 5
RSS 2835156.750 1597457.194 1552923.616 1538096.513 1507888.476 1659993.500
BIC 1318.242 1270.084 1276.467 1284.718 1291.944 1310.765

R paketet bcp
Ennan metod är funktionen bcp() från paketet bcp (Bayesian Change Point) som visar sannolikheterna för att de olika observationerna är en förändringspunkter. Den mest sannolika punkten har markerats med ett gäng utropstecken.


> library(bcp)
> Nile.bcp <- bcp(Nile)
> summary(Nile.bcp)
> Nile.bcp

Probability of a change in mean, Posterior Mean,
and SD for each position:

Probability Mean SD
1 0.016 1085.8 14.15
2 0.012 1085.6 14.06
3 0.014 1085.1 13.61
4 0.018 1085.5 13.80
5 0.018 1085.2 14.49
6 0.026 1084.8 15.67
7 0.046 1079.2 37.90
8 0.038 1089.4 29.73
9 0.062 1093.0 40.96
10 0.058 1082.8 26.11
11 0.028 1075.1 25.60
12 0.016 1072.0 27.79
13 0.016 1071.6 27.32
14 0.010 1071.0 28.02
15 0.022 1070.6 28.80
16 0.026 1069.5 30.92
17 0.022 1070.0 30.39
18 0.044 1067.7 38.25
19 0.088 1071.5 35.46
20 0.062 1084.1 29.90
21 0.054 1091.7 29.46
22 0.038 1098.5 31.98
23 0.026 1102.0 32.26
24 0.024 1104.9 34.27
25 0.028 1105.3 36.10
26 0.136 1102.9 36.36
27 0.184 1071.1 77.60
28 0.618 1029.3 108.61 !!!!
29 0.078 869.8 66.26
30 0.034 856.8 40.52
31 0.042 852.5 28.16
32 0.020 849.5 14.42
33 0.022 850.4 14.29
34 0.008 849.9 14.03
35 0.018 849.6 14.19
36 0.018 850.2 14.20
...


Den 28:e observationen har alltså högst sannolikhet. Den sista observationen har sannolikhet 1 och som jag ignorerar. En graf plottas också där man tydligt ser förändringen.

> which.max(Nile.bcp$posterior.prob[-c(100)])
[1] 28

> plot(Nile.bcp)

Grafen:


(För övrigt är detta blogganteckning #1200. Yeay!)

Uppdatering
Håkan Karlsson (aka hakke) uppmärksammade uppmärksamt att "normalfallet" och "extremfallet" hade förväxlats i pratbeskrivningen av hur S uppdateras. Tack hakke.

Posted by hakank at augusti 21, 2007 07:39 EM Posted to Statistik/data-analys

Comments

Detta påminner om en metod jag jobbade med för många år sen, AFMM, Adaptive Forgetting through Multiple Models, när man gjorde tidsserieanalys rekursivt, med användandet av linjära och stationära modeller, men då data faktiskt var icke-stationära. Funkade rätt bra för att snabbt hänga med i snabba förändringar av processen!

Hittade ingen riktigt bra länk, men här kanske man får en uppfattning om det: http://www.mathworks.de/products/sysid/demos.html?file=/products/demos/shipping/ident/iddemo5.html

Posted by: thebe at oktober 6, 2007 10:09 EM

Thebe: Tack för tipset.

Posted by: Håkan Kjellerstrand at oktober 7, 2007 09:43 FM