需 要 予 測 のための 統 計 モテ ルの 研 究 異 常 値 検 知 のための 基 本 的 モテ ルの 考 察 平 成 26 年 3 月 東 京 大 学 大 学 院 情 報 理 工 学 系 研 究 科 教 授 博 士 課 程 修 士 課 程 竹 村 彰 通 小 川 光 紀 笹 井 健 行 特 定 非 営 利 活 動 法 人 ヒ ュー コミュニケーションス 副 理 事 長 小 松 秀 樹 主 任 研 究 員 池 谷 貞 彦 研 究 員 河 内 敦 雄

1 1 3 1.1......................... 3 1.2........................... 3 2 5 2.1......................... 5 2.1.1.......................... 5 2.1.2.......................... 5 2.1.3................. 6 2.2 (AR)...................... 6 3 8 3.1................. 8 3.2................. 8 3.3...................... 9 4 Change Finder 10 4.1............................... 10 4.2................. 11 4.3 Change Finder...................... 12 4.4 Change Finder....................... 13 4.5 SDAR....................... 15 5 18 5.1 Change Finder....... 18 5.1.1 1....................... 19 5.1.2 2....................... 20 5.1.3 3....................... 21 5.1.4 4....................... 22 5.1.5 5....................... 23 5.1.6 1........................ 24 5.1.7 2........................ 25 5.1.8 3........................ 26 5.1.9 4........................ 27

2 5.1.10 5........................ 28 5.1.11 6........................ 29 5.2 Change Finder............ 29 5.3 Change Finder....... 30 5.3.1 1....................... 30 5.3.2 5....................... 31 5.3.3 1........................ 31 5.4 Change Finder............ 32 6 33 A 35 A.1 -............. 35 A.1.1....................... 35 A.1.2.................. 35 A.2 -............. 36 A.3.......................... 37 B 38 C 43

3 1 1.1 ARIMA ( ) ARIMA ( ) ARIMA ( ) 1.2 (AR ) AR

1 4 Change Finder AR Change Finder

5 2 AR [1] 2.1 2.1.1 T t T x t p (x t ) {p (x t ) t T t T = {1, 2,... 2.1.2 {x t {x t E (x t ) = µ Cov (x t, x t+1 ) = γ ( h ) {x t V (y t ) = γ (0) ρ (h) = γ (h) γ (0) = γ (t) V (y t ) = ρ ( h) h ρ (h), h = 0, 1, 2,... ρ (h) 1 1

2 6 2.1.3 S = Cov (x t, x t+h ) = γ (h) h= h= S 2.2 (AR) (AR) {x t t T x t = ϕ 1 x t 1 + + ϕ p x t p + ε t, {ε t i.i.d.n ( 0, σ 2) (2.1) {x t p AR(p) AR(p) AR(p) ϕ (x) = 1 ϕ 1 x ϕ 2 x 2 ϕ p x p = 0 1 AR(p) (2.1) y t j (j > 0) γ (j) = ϕ 1 γ (j 1) + ϕ 2 γ (j 2) + + ϕ p γ (j p) (2.2) γ (0) p ρ (j) = ϕ 1 ρ (j 1) + ϕ 2 ρ (j 2) + + ϕ p ρ (j p), (j > 0) (2.3) λ 1, λ 2,, λ p c 1, c 2,, c p ρ (j) = p c i i=1 λ j i c 1, c 2,..., c p (2.3) j = 1, 2,..., p 1 ρ (1), ρ (2),..., ρ (p 1)

2 7 ρ (1) = ϕ 1 + ϕ 2 ρ (1) + + ϕ p ρ (p 1) ρ (2) = ϕ 1 ρ (1) + ϕ 2 + + ϕ p ρ (p 2). ρ (p 1) = ϕ 1 ρ (p 2) + ϕ 2 ρ (p 3) + + ϕ p ρ (1). c 1, c 2,.., c p 1 = c 1 + c 2 + + c p ρ (1) = c 1 + c 2 + + c p λ 1 λ 2 λ p. ρ (p 1) = c 1 λ p 1 1 + c 2 λ p 1 2 + + c p λ p 1 p AR(p) ϕ 1, ϕ 2,..., ϕ p (2.3) j = 1, 2,..., p p. ρ (1) = ϕ 1 + ϕ 2 ρ (1) + + ϕ p ρ (p 1) ρ (2) = ϕ 1 ρ (1) + ϕ 2 + + ϕ p ρ (p 2). ρ (p) = ϕ 1 ρ (p 1) + ϕ 2 ρ (p 2) + + ϕ p. 1 ρ (1) ρ (p 1) ϕ 1 ρ (1) ρ (1) 1 ρ (p 2)....... =.... ρ (p 1) ρ (p 2) 1 ρ (p) ϕ p

8 3 [2] 3.1 x t, (t = 1, 2,...) µ x τ τ = 7 γ xx (τ) = E [(x t µ x ) (x t+τ µ x )] 3.2 A B A B x t, y t (t = 1, 2,...) µ x, µ y τ x t = 1 y t = 2 τ = 7 1

3 9 2 γ xy (τ) = E [(x t µ x ) (y t+τ µ y )] 3.3

10 4 Change Finder [3, 4] Change Finder AR AR AR [3, 4] 4.1 x 1 x 2... p x t 1 = x 1 x 2...x t 1 p ( x t x t 1) p p (1), p (2) p ( x t x t 1) = p (1) ( x t x t 1), t < a, p ( x t x t 1) = p (2) ( x t x t 1), t a. t = a -

4 Change Finder 11 ( D p (2) p (1)) [ 1 := lim n n E p log p(2) (x n ] ), (2) p (1) (x n ) p (1) = p (2) D = 0 - - 1. jumping mean p (1), p (2) (µ 1, σ), (µ 2, σ) D ( p (1) p (2)) = (µ 1 µ 2 ) 2 / ( 2σ 2) 2. jumping variance p (1), p (2) (µ, σ 1 ), (µ, σ 2 ) D ( p (1) p (2)) = 1 2 ( σ 2 2 σ 2 1 1 log σ2 2 σ 2 1 ) 2 4.2 x n 1 = x 1,..., x n t x i, (i = 1,, n) d x t = x 1,..., x t x n t = x t+1,..., x n t x t 1 t k = x t k,..., x t 1 x t p ( [ ] x t x t 1 ) 1 t k = exp (x t w t ) T Σ 1 (x t w t ) (2π) d/2 1/2 Σ 2 Σ AR k α 1, α 2,..., α k, µ w t = k α i (x t i µ) + µ i=1

4 Change Finder 12 AR w t ŵ t I (x n 1 ) = n x t ŵ t 2. t=1 t I (x t 1) Error1 t I ( x n t+1) Error2 t = t δ > 0 I (x n 1 ) ( I ( ) ( x t 1 + I x n t +1)) > δ n t Change Finder 4.3 Change Finder Change Finder Change Finder Step 1 x 1, x 2,... {p t (x) : t = 1, 2,... p t 1 (x) x t 1 = x 1,..., x t 1 t x t Score (x t ), t = 1, 2,... Step 2 W > 0 Step1 {y t : t = 1, 2,... Step {Score (x i ) : i = t W + 1,..., t W - y t

4 Change Finder 13 y t = 1 W t t=t W +1 Score (x i ). Step 3 Step1 Step2 y t, t = 1, 2,... {q t (y t ) : t = 1, 2,... t y t Score (y t ), t = 1, 2,... W > 0 t W - Score (t) = 1 t W Score (y i ) i=t T W +1 Change Finder Step2 Step3 Step2 T Change Finder W W W 4.4 Change Finder Change Finder AR Step 1 x t, t = 1, 2,... SDAR {p t (x) : t = 1, 2,... SDAR

4 Change Finder 14 t x t Score (x t ) = log p t 1 (x t ) Score (x t ) = d (p t 1, p t ) = ( p (xt 1 ) ) 2 p (x t ) d AR z t = k A i z t i + ε i=1 A i ε 0 Σ {x t : t = 1, 2,... x t = z t + µ x t 1 t k = x t k...x t 1 x t 1 t k x t p ( ( x t x t 1 t k : θ) 1 = exp 1 ) (2π) k/2 Σ 1/2 2 (x t w) T Σ 1 (x t w), k w = A i (x t i µ) + µ i=1 θ = (A 1,..., A k, µ, Σ) x t θ θ t p t = p (, θ t ) θ (SDAR ) Step 2 k (1 r) t i log p ( x i x i 1, θ ). i=1 Step1 {y t : t = 1, 2,... AR SDAR {q t : t = 1, 2,... W > 0 t W -

4 Change Finder 15 Score (t) = 1 t W ( log q i 1 (y i )) i=t W +1 Step1 Score (t) = 1 t W d (q i 1, q i ) i=t W +1 Score (t) t 4.5 SDAR AR SDAR (Sequencially Discounting AR model learning) n x 1,..., x n x t = z t + µ z 1,..., z n k p (z t θ) t=1 n t=k+1 p ( z t z t 1 t k : θ). k log p (z t θ) + t=1 n t=k+1 log p ( z t z t 1 t k : θ). n k AR ( (n k) log (2π) 1/2 Σ 1/2) 1 2 n t=k+1 ( z t ) T k w i z t i Σ (z 1 t w w i (i = 1,..., k) i=1 k w i C j i = C j (j = 1,..., k). i=1 C j ) k w i z t i. i=1

4 Change Finder 16 C j = 1 n k = 1 n k n t=k+1 n t=k+1 z t z T t j (x t µ) (x t j µ) T. s C s = C T s. µ Σ ˆµ = 1 n k n t=k+1 C j µ ˆµ ŵ 1,..., ŵ k x t, ˆΣ = 1 n k = 1 n k n t=k+1 n t=k+1 ( ( z t ) k ŵ i z t i (z t i=1 x t ˆµ ) T k ŵ i z t i i=1 ) ( k ŵ i (x t i µ) x t ˆµ i=1 k ŵ i (x t i µ) SDAR (1 r) : r 0 < r < 1 r SDAR AR SDAR AR SDAR i=1 ) T

4 Change Finder 17 SDAR :0 < r < 1 ˆµ, C j, ŵ j (j = 1, 2,...), ˆΣ t (= 1, 2,...) x t ˆµ := (1 r) ˆµ + rx t, C j := (1 r) C j + r (x t ˆµ) (x t j ˆµ) T. k w i C i i = C j (j = 1,..., k). i=1 ŵ 1,..., ŵ k k ˆx t := ŵ i (x t i ˆµ) + ˆµ, i=1 ˆΣ := (1 r) ˆΣ + r (x t ˆx t ) (x t ˆx t ) T

18 5 Change Finder 11 Change Finder 5.1 Change Finder SDAR AR 2 0.02 Score [x (t)] E [Score [x (t)]] σ [Score [x (t)]] E [Score [x (t)]] + 4σ [Score [x (t)]] Score Score [x (t)] x (t) E [x (t)] 2σ [x (t)] E [x (t)] + 2σ [x (t)] SDAR

5 19 5.1.1 1 yasai1 sales quantity 0 20 60 5.1.1: 1 yasai1 score 0 2 4 6 8 5.1.2: 1 1 46

5 20 5.1.2 2 yasai2 sales quantity 0 100 200 5.1.3: 2 yasai2 score 0 10 20 30 5.1.4: 2 2 59 13, 17, 25, 52, 29 13, 59

5 21 5.1.3 3 yasai3 sales quantity 0 100 200 300 5.1.5: 3 yasai3 score 0 2 4 6 8 5.1.6: 3 3 45, 46, 59, 60 45, 46, 59, 60

5 22 5.1.4 4 yasai4 sales quantity 0 100 200 5.1.7: 4 yasai4 score 0 4 8 5.1.8: 4 4 11

5 23 5.1.5 5 yasai5 sales quantity 0 100 200 300 5.1.9: 5 yasai5 score 0 2 4 6 8 5.1.10: 5 5

5 24 5.1.6 1 niku1 sales quantity 0 5 15 5.1.11: 1 T = 2 AR 3 niku1 score 0 5 10 20 5.1.12: 1 1 62 63

5 25 5.1.7 2 niku2 sales quantity 0 10 30 5.1.13: 2 niku2 score 0 2 4 6 8 5.1.14: 2 2 7, 67

5 26 5.1.8 3 niku3 sales quantity 0 400 800 5.1.15: 3 niku3 score 0 2 4 6 8 5.1.16: 3 AR 3 3

5 27 5.1.9 4 niku4 sales quantity 0 20 40 5.1.17: 4 niku4 score 0 5 10 15 5.1.18: 4 4 60 60

5 28 5.1.10 5 niku5 sales quantity 0 40 80 120 5.1.19: 5 niku5 score 0 2 4 6 5.1.20: 5 5

5 29 5.1.11 6 niku6 sales quantity 0 200 400 5.1.21: 6 niku6 score 0 4 8 12 5.1.22: 6 6 5.2 Change Finder AR

5 30 7 2 3 1 5.3 Change Finder 11 1 2 1 Change Finder SDAR 5.3.1 1 yasai1 yasai1 sales quantity 0 20 60 score 0 20 40 60 5.3.1: AR 3 0.02 4 AR 2 0.02 4 1 50

5 31 5.3.2 5 yasai5 yasai5 sales quantity 0 100 200 300 score 0 40 80 5.3.2: AR 2 0.02 3 AR 2 0.002 3 5 50 5.3.3 1 niku1 niku1 sales quantity 0 5 15 score 0 2 4 6 8 5.3.3: AR 3 0.02 4 AR 3 0.002 4 1 60

5 32 5.4 Change Finder

33 6 Change Finder Change Finder

34 [1] 2006 [2] J. P. Peckmann, S. Oliffson Kamphorst, and D. Ruelle, Recurrence Plots of Dynamical Systems, Europhysics Letters, Vol. 4, No. 9, pp. 973-977, 1987. [3] 2009 [4] J. Takeuchi and K. Yamanishi, A Unifying Framework for Detecting Outliers and Change Points from Time Series, IEEE Transaction on Knowledge and Data Engineering, Vol. 18, No. 4, pp482-492, 2006. [5] C. M. 2008.

35 A A.1 - - [5] A.1.1 x h (x) x x p (x) 1 h (x) p (x) 2 x, y h (x, y) = h (x) + h (y) p (x, y) = p (x) p (y) h (x) p (x) h (x) = log p (x) x h (x) 2 A.1.2 H [x] = x p (x) log p (x). N i n i n i n i! N

A 36 W = N! i n i! H = 1 N ln W = 1 N ln N! 1 N ln n i! n i /N N i i n i = N H = lim N ln N! N ln N N i ( ni ) ( ni ) ln = N N i p i ln p i p i = lim N (n i /N) i n i /N W H [x] = p (x) ln p (x) dx A.2 - - p (x) q (x) p (x) q (x) q (x) x x KL (p q) = = p (x) ln q (x) dx { q (x) p (x) ln p (x) dx ( ) p (x) ln p (x) dx p (x) q (x) -

A 37 A.3 Change Finder - p (x) θ q (x θ) θ p (x) q (x θ) - p (x) p (x) x n (n = 1, 2,..., N) - KL (p q) 1 N N { ln q (x n θ) + ln p (x n ) n=1 - Change Finder

38 B Change Finder R ###.txt ###scorext ###score ##### ###########step1 x <- read.table("yasai1.txt") ### n <- nrow(x) scorext <-matrix(0, nrow=n, ncol=1) pt <-matrix(0, nrow=n, ncol=1) ### Tst1 <- 3 Tst3 <- 5 ###AR kst1 <- 3 kst3 <- 2 ### rst1 <- 0.02 rst3 <- 0.002 xmean <- mean(x[,1]) Cst1 <- matrix(0, nrow=kst1+1, ncol=n) Cst1e <- matrix(0, nrow=kst1, ncol=1) Cmatst1 <- matrix(0, nrow=kst1, ncol=kst1) nmatst1 <- matrix(0, nrow=kst1, ncol=kst1) wst1 <- matrix(0, nrow=kst1, ncol=n) xt <- matrix(0, nrow=n, ncol=1) xtemp <- 0 #must1 <- xmean must1 <- 0 sigmast1 <- matrix(0, nrow=n, ncol=1)

B 39 kstartst1 <- kst1+1 for (t in kstartst1 : n ) { must1 <- (1-rst1)*must1 + rst1*x[t,1] for (j in 1 : kstartst1){ Cst1[j,t] <- (1-rst1)*Cst1[j,t-1]+rst1*(x[t,1]-must1)*(x[t-j+1,1]-must1) Cmatst1 <- matrix(0, nrow=kst1, ncol=kst1) for (j in 1 : kst1){ nmatst1 <- matrix(0, nrow=kst1, ncol=kst1) for (i in 1 : kst1){ for (l in 1 : kst1){ if (abs(l-i) == j-1){ nmatst1 <- matrix(0, nrow=kst1, ncol=kst1) nmatst1[l,i] <- 1 Cmatst1 <- Cmatst1 + Cst1[j,t]*nmatst1 for (i in 1:kst1){ Cst1e[i]<-Cst1[i+1,t] wst1[,t] <- solve(cmatst1,cst1e) for (i in 1 : kst1){ xtemp <- xtemp + wst1[i,t]*(x[t-i,1]-must1) xt[t] <- xtemp + must1 xtemp <- 0 if (t == kstartst1){ sigmast1[t] <- (x[t,1]-xt[t])*(x[t,1]-xt[t]) if (t!= kstartst1){ sigmast1[t] <- (1-rst1)*sigmast1[t-1] + rst1*(x[t,1]-xt[t])*(x[t,1]-xt[t]) pt[t] <- exp((-(xt[t]-x[t,1])^2)/(2*sigmast1[t]))*(1/(2*pi*sqrt(sigmast1[t]))) scorext[t] <- -log(pt[t])

B 40 #######step2 y <- matrix(0, nrow=n, ncol=1) ytemp <- 0 step2kst1 <- kst1+tst1 for (t in step2kst1 : n ) { kstartst2 <- t-tst1+1 for (i in kstartst2 : t){ ytemp <- ytemp + scorext[i] y[t] = ytemp/tst1 ytemp <- 0 #########step3 scoreyt <-matrix(0, nrow=n, ncol=1) qt <-matrix(0, nrow=n, ncol=1) ytempst3 <- 0 Cst3 <- matrix(0, nrow=kst3+1, ncol=n) Cmatst3 <- matrix(0, nrow=kst3, ncol=kst3) Cst3e <- matrix(0, nrow=kst3, ncol=1) nmatst3 <- matrix(0, nrow=kst3, ncol=kst3) wst3 <- matrix(0, nrow=kst3, ncol=n) yt <- matrix(0, nrow=n, ncol=1) ysum <- sum(y) #must3 <- ysum/(n-tst1-kst1+1) must3 <- 0 sigmast3 <- matrix(0, nrow=n, ncol=1) kstartst3 <- kst1+tst1+1+kst3 kstartst32 <- kst3+1 for (t in kstartst3 : n ) { must3 <- (1-rst3)*must3 + rst3*y[t]

B 41 for (j in 1 : kstartst32){ Cst3[j,t] <- (1-rst3)*Cst3[j,t-1] + rst3*(y[t,1]-must3)*(y[t-j+1,1]-must3) Cmatst3 <- matrix(0, nrow=kst3, ncol=kst3) for (j in 1 : kst3){ nmatst3 <- matrix(0, nrow=kst3, ncol=kst3) for (i in 1 : kst3){ for (l in 1 : kst3){ if (abs(l-i) == j-1){ nmatst3 <- matrix(0, nrow=kst3, ncol=kst3) nmatst3[l,i] <- 1 Cmatst3 <- Cmatst3 + Cst3[j,t]*nmatst3 for (i in 1:kst3){ Cst3e[i]<-Cst3[i+1,t] wst3[,t] <- solve(cmatst3,cst3e) ytempst3 <- 0 for (i in 1 : kst3){ ytempst3 <- ytempst3 + wst3[i,t]*(y[t-i]-must3) yt[t] <- ytempst3 + must3 if (t==kstartst3){ sigmast3[t] <- (y[t]-yt[t])*(y[t]-yt[t]) if (t!=kstartst3){ sigmast3[t] <- (1-rst3)*sigmast3[t-1] + rst3*(y[t]-yt[t])*(y[t]-yt[t]) qt[t] <- exp(((yt[t]-y[t])^2)/(2*sigmast3[t])*(-1))*(1/(2*pi*sqrt(sigmast3[t]))) scoreyt[t] <- -log(qt[t]) score <- matrix(0, nrow=n, ncol=1) scoretemp <- 0 winsst3 <- kst1+tst1+tst3 for (t in winsst3 : n ) {

B 42 kstartscore <- t - Tst3 +1 for (i in kstartscore : t){ scoretemp <- scoretemp + scoreyt[i] score[t] = scoretemp/tst3 scoretemp <- 0 #plot(x[,1],xlab="",ylab="sales quantity",type="o",main="niku6") #par(new=t) plot(score,xlab="",ylab="score",col="red",main="niku1",type="o")

43 C 2013 10 17 12 25 70