廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

dlm による LWR モデルのパラメータ推定(渋滞の予測)

シンプルな LWR モデルの近似解 - 廿TT のモデルでもう少し遊んでみる。

観測点が 6 個あって、それぞれが時間ごとにトラフィックの量を計測していると考える。

上記で離散化した、

\displaystyle \rho^{m+1}_{n} = (1+r)\rho^{m}_{n} - r \rho^{m}_{n-1}

状態方程式である。

\partial q(\rho)/\partial \rho が一定という仮定はちょっと不自然な気がするが、下手にいじると状態方程式が線形じゃなくなるので、このままいく。

ただし、道路の入り口からは常に一定数自動車の流入があるとして、

\displaystyle \rho^{m+1}_{0} = \rho^{m}_{0}

とする。

これは境界条件で左端の微分が 0 と置くことと対応する。

m=0,1,\ldots,50 までの観測値を使って、10期先まで予測する。

未知パラメータ r の真値は 0.1、システムノイズの分散は 1、観測ノイズの分散は25として、シミュレーションでデータを生成した。

観測値は整数になるように丸めた。

R の dlm パッケージを使ってパラメータを推定した。

各パラメータの推定値はこんな感じ。

> print(round(exp(fit1$par),2))
[1]  0.10 23.15  1.44

各観測地点の予測結果(スムージングの結果)のプロットはこんな感じ。点線の右が予測値です。

f:id:abrahamcow:20170924024226p:plain

以下に R のコードを貼ります。

library("tidyverse")
library("dlm")
diffeq <- function(q0, m, u0, sim=FALSE){ 
  n <- length(u0)
  r <- -q0
  A <- diag(1+r,n)
  A[1,1] <- 1
  A[cbind(2:n,1:(n-1))] <- -r
  b <- matrix(NA,m,n)
  b0 <- u0
  b[1,] <- u0
  for(i in 2:m){
    b1 <- A %*% b0
    if(sim){
      b0 <- b1+rnorm(n,0,1)
    }else{
      b0 <- b1
    }
    b[i,] <- b0
  }
  colnames(b) <- seq(0,by=1,length.out=n)
  rownames(b) <- seq(0,by=1,length.out=m)
  b
}
set.seed(10)
p <- gtools::rdirichlet(1,rep(1,6))
u0 <- drop(rmultinom(1,1000,p))
out <-diffeq(q0=0.1,61,u0=u0,sim = TRUE)
dat <- round(out + rnorm(length(out),0,5))
dat2 <- dat
dat2[51:61,] <- NA

ff <- diag(1,6)
build1 <- function(theta){
  r <- -exp(theta[1])
  gg <- diag(1+r,6)
  gg[1,1] <- 1
  gg[cbind(2:6,1:5)] <- -r
  dlm(FF=ff,
      V=diag(exp(theta[2]),6),
      GG=gg,
      W = diag(exp(theta[3]),6),
      m0 = numeric(6),
      C0 = diag(1e+7,6))
}

fit1 <- dlmMLE(dat2, parm=numeric(3), build1, method = "Nelder-Mead")

print(round(exp(fit1$par),2))

model1 <- build1(fit1$par)
Filt1 <- dlmFilter(dat2,model1)
Smooth1 <- dlmSmooth(Filt1)
smooth1 <-Smooth1$s
colnames(smooth1) <- 0:5

smooth_df <-as_data_frame(cbind(time=1:61,smooth1[-1,])) %>%
  gather(x,rho,-time)

dat_df <-as_data_frame(cbind(time=1:61,dat)) %>% 
  gather(x,rho,-time)

p1 <-ggplot(dat_df,aes(x=time,y=rho))+
  geom_point(alpha=0.5,size=2)+
  geom_line(data = smooth_df,colour="royalblue",size=1)+
  facet_wrap(~x,scales="free_y")+
  geom_vline(xintercept=50,linetype=2)
print(p1)

おすすめ書籍など

こちらは KFAS パッケージを使っている。

Rによるベイジアン動的線形モデル (統計ライブラリー)

Rによるベイジアン動的線形モデル (統計ライブラリー)

こちらは dlm パッケージの開発者による本。