ts1011

1 paper 範例

• Time Update (“Predict”)

\[\hat{x}_k^- ~=~ \hat{x}_{k-1}\]

\[P_k^- ~=~ P_{k-1}+Q\]

• Measurement Update (“Correct”)

\[K_k ~=~ \frac{P_k^-}{P_k^-+R}\]

\[\hat{x}_k ~=~ \hat{x}_k^-+K_k(Z_k-\hat{x}_k^-)\]

\[P_k ~=~ (1-K_k)P_k^-\]

---

KF <- function(r=1){
  A=1
  B=0
  samples=100
  timeidx=1:100
  timeidx=t(timeidx)
  H=1
  xbar=c()
  #set.seed(23)
  y=-0.37727+matrix(rnorm(100, 0,0.025))
  q=1
  xhat=0
  P0=1
  r=r
  u=matrix(0,samples+1,1)

  nu=dim(u)[2]
  ny=dim(y)[2]
  nx=1
  nv=1
  nw=1

  Pxbar=P0
  xhat_data=matrix(0,samples+1,1)
  Pmat=matrix(0,samples+1,nx*nx)
  yidx=1
  
  Px=Pxbar
  uk1=0
  yidx=1
  
  xhat_data[1]=t(xhat)
  
  
  for(k in 1:length(y)){
    ## Time Update ("Predict")
    # (1) Project the state ahead
    if(nu>0) uk1=u[k+1]
    xbar=A*xhat+B*uk1
    # (2) Project the error covariance ahead
    Pxbar=A*Px*t(A)+q
    
    ###
    ## Measurement Update ("Correct")
    if(k<=timeidx[100]&timeidx[yidx]==k){
      # (1) Compute the Kalman gain
      K=Pxbar*t(H)/(H*Pxbar*t(H)+r)
      
      # (2) Update estimate with measurement zk
      xhat = xbar + K*(y[yidx] - H*xbar)
      # (3) Update the error covariance
      Px = Pxbar-K*H*Pxbar
      yidx=yidx+1
    }
    xhat_data[k+1,]=xhat
    Pmat[k+1,]=Px
  }
  plot(y,ylim = c(-0.5,-0.28), main=paste0("R=",r))
  matplot(xhat_data[2:101], add = T, type = "l", col=2)
  abline(h = -0.37727)
}

par(mfcol=c(1,2))
set.seed(100)
KF(10)
set.seed(100)
KF(1)

set.seed(100)
KF(0.1)
set.seed(100)
KF(0.01)

---

2 Kalman Filter

KF2 <- function(r){
  y <- read.csv("/home/yuan/Desktop/meeting/處理後來店人數資料.csv")$x
  A=1
  B=0
  samples=2548
  timeidx=1:2548
  timeidx=t(timeidx)
  H=1
  xbar=c()
  y=t(y)
  # set.seed(23)
  # y=-0.37727+matrix(rnorm(100, 0,0.025))
  q=1
  xhat=0
  P0=1
  r=r
  u=matrix(0,samples+1,1)

  nu=dim(u)[2]
  ny=dim(y)[2]
  nx=1
  nv=1
  nw=1

  Pxbar=P0
  xhat_data=matrix(0,samples+1,1)
  Pmat=matrix(0,samples+1,nx*nx)
  yidx=1

  Px=Pxbar
  uk1=0
  yidx=1

  xhat_data[1]=t(xhat)


  for(k in 1:length(y)){
    # (1) Project the state ahead
    if(nu>0) uk1=u[k+1]
    xbar=A*xhat+B*uk1
    
    # (2) Project the error covariance ahead
    Pxbar=A*Px*t(A)+q
    
    ###
    ## Measurement Update ("Correct")
    if(k<=timeidx[2548]&timeidx[yidx]==k){
      # (1) Compute the Kalman gain
      K=Pxbar*t(H)/(H*Pxbar*t(H)+r)
      # (2) Update estimate with measurement zk
      xhat = xbar + K*(y[yidx] - H*xbar)
      # (3) Update the error covariance
      Px = Pxbar-K*H*Pxbar
      yidx=yidx+1
    }
    xhat_data[k+1,]=xhat
    Pmat[k+1,]=Px
  }
  plot(t(y), type = "l", main=paste0("R=",r))
  matplot(xhat_data[2:2549], add = T, type = "l", col=2) 
}

par(mfcol=c(1,2))
KF2(r = 10)
KF2(r = 1)

KF2(r = 0.1)
KF2(r = 0.01)

---

3 Holt-Winter

\(\hat{y} ~=~ l_t+hb_t+s_{t-m+h_m^+}\)

\(l_t ~=~ \alpha(y_t-s_{t-m})+(1-\alpha)(l_t+b_{t-1})\)

\(b_t ~=~ \beta^*(l_t-l_{t-1})+(1-\beta^*)b_{t-1}\)

\(s_t ~=~ \gamma(y_t-l_{t-1}-b_{t-1})+(1-\gamma)s_{t-m}\)

where \(h_m^+ ~=~ \lfloor (h-1) ~~~ mod ~~~ m \rfloor +1\)

---

4 exponential smoothing with regressor

\(y_t ~=~ \begin{bmatrix} 1 & Z_t \end{bmatrix} \begin{bmatrix} l_{t-1}\\P_1 \end{bmatrix} +\varepsilon_t\)

\(\begin{bmatrix} l_t \\P_t \end{bmatrix} ~=~ \begin{bmatrix} 1&0\\0&1 \end{bmatrix} \begin{bmatrix} l_{t-1}\\P_t \end{bmatrix} + \begin{bmatrix} \alpha\\0 \end{bmatrix} \varepsilon_t\)

---

\(K_k ~=~ \frac{P_k^-H^T}{HP_k^-H^T+R}\)

\(~=~ P_k^- \begin{bmatrix} 1\\ Z_t \end{bmatrix} \begin{bmatrix} \begin{bmatrix}1 & Z_t \end{bmatrix} P_k^- \begin{bmatrix} 1 \\ Z_t \end{bmatrix} +\sigma^2 \end{bmatrix}^{-1}\)