ts.length <- 100

ts
100 
a <- rep(0.5, ts.length)

print(length(a))
print(a)

[1] 100
  [1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
 [19] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
 [37] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
 [55] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
 [73] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
 [91] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

x  <- rep(0, ts.length)
v  <- rep(0, ts.length)
for (ts in 2:ts.length) { # 2 에서 100 까지 (ts - 1의 수행 때문)
  x[ts] <- v[ts - 1] * 2 + x[ts - 1] + 1/2 * a[ts-1] ^ 2 
  x[ts] <- x[ts] + rnorm(1, sd = 20)
  v[ts] <- v[ts - 1] + 2 * a[ts-1]
}

par(mfrow = c(3, 1))
plot(x,            main = "Position",     type = 'l')
plot(v,            main = "Velocity",     type = 'l')
plot(a,            main = "Acceleration", type = 'l')

z <-  x + rnorm(ts.length, sd = 300)
plot (x, ylim = range(c(x, z)))
lines(z)

kalman.motion <- function(z, Q, R, A, H) {
  dimState = dim(Q)[1]
  
  xhatminus <- array(rep(0, ts.length * dimState), 
                    c(ts.length, dimState)) 
  xhat      <- array(rep(0, ts.length * dimState), 
                    c(ts.length, dimState)) 
 
  Pminus <- array(rep(0, ts.length * dimState * dimState), 
                 c(ts.length, dimState, dimState))
  P      <- array(rep(0, ts.length * dimState * dimState), 
                 c(ts.length, dimState, dimState))  
    
  K <- array(rep(0, ts.length * dimState), 
            c(ts.length, dimState)) # Kalman gain # 칼만이득
  
  # intial guesses = starting at 0 for all metrics
  # 초기 추측 = 모든 지표는 0으로 시작합니다
  xhat[1, ] <- rep(0, dimState)
  P[1, , ]  <- diag(dimState)
  
  for (k in 2:ts.length) {
    # 시간의 갱신
    xhatminus[k, ] <- A %*% matrix(xhat[k-1, ])
    Pminus[k, , ] <- A %*% P[k-1, , ] %*% t(A) + Q
    
    K[k, ] <- Pminus[k, , ] %*% H %*% 
                            solve( t(H) %*% Pminus[k, , ] %*% H + R )
    xhat[k, ] <- xhatminus[k, ] + K[k, ] %*% 
                            (z[k]- t(H) %*% xhatminus[k, ])
    P[k, , ] <- (diag(dimState)-K[k,] %*% t(H)) %*% Pminus[k, , ]
  }    
  
  ## we return both the forecast and the smoothed value
  ## 예측과 평활화된 값 모두를 반환합니다
  return(list(xhat = xhat, xhatminus = xhatminus))
}

R <- 10^2 ## measurement variance - this value should be set 
          ## according to known physical limits of measuring tool 
          ## we set it consistent with the noise we added to x 
          ## to produce x in the data generation above
          ## 측정 분산 - 이 값은 측정 도구에 대하여 알려진,
          ## 물리적인 한계에 따라서 설정되어야만 합니다.
          ## x에 더해진 노이즈와 일관성있게 설정합니다.
Q <- 10   ## process variance - usually regarded as hyperparameter
          ## to be tuned to maximize performance
          ## 과정의 분산 - 일반적으로, 성능의 최대화를 위해서
          ## 조정되어야 하는 하이퍼파라미터로서 취급됩니다

## dynamical parameters
## 동적 파라미터
A <- matrix(1) ## x_t = A * x_t-1 (how prior x affects later x) (사전의 x가 나중의 x에 얼마나 영향을 미치는지)
H <- matrix(1) ## y_t = H * x_t   (translating state to measurement) (상태를 측정으로 변환)

## run the data through the Kalman filtering method
## 칼만필터 방법으로 데이터를 넣고 돌립니다
xhat <- kalman.motion(z, diag(1) * Q, R, A, H)[[1]]

result = kalman.motion(z, diag(1) * Q, R, A, H)
xhat = result[[1]]
xhatminus = result[[2]]

z <-  x + rnorm(ts.length, sd = 300)
plot (x, type="l", lty=2, ylim = range(c(x, z)))
lines(z, lty=1)
lines(xhat, lty=3)
lines(xhatminus, lty=4)