廿TT

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

ポアソン・ガンマモデルによる渋滞のシミュレーションと予測

今日の川柳

モチベーション

https://www.jcca.or.jp/kaishi/268/268_toku1.pdf に渋滞の伝播を表した図が乗っています。

f:id:abrahamcow:20190608152002p:plain
https://www.jcca.or.jp/kaishi/268/268_toku1.pdf より

LWRモデルを線形にしたやつ

では車が順々に移動していくような状況は再現できますが、渋滞の発生と伝播はシミュレーションできていませんでした。

そこで渋滞の発生と伝播を再現できるようななるべくシンプルなモデルを考えたいと思います。

上の図では速度が観測されていましたが、時速を逐一観測するのは難しいと思うので、単位時間ごとの単位区間ごと自動車の台数が観測されているような状況を仮定します。

モデル

時刻 t、区間 m の車の台数は平均 \lambda_{t,m}ポアソン分布に従うとします。

混雑具合 \lambda_{t,m} は、一時点前の一地点先の混雑と、一時点前の同じ地点の混雑に影響されて決まると考えられるので、以下の生成モデルを仮定します。

\lambda_{t,m} \sim \mathrm{Gamma}(\lambda_{t-1,m}+\lambda_{t-1,m+1},2)

このモデルで簡単なシミュレーションをやってみたのが下の図です。

f:id:abrahamcow:20190608153115p:plain

帯状の伝播が発生しているのがわかると思います。

推定

ギブスサンプリングでデータから \lambda_{t,m} を推定してみたいと思います。

一時点目の観測に関しては前の時点の情報がないので、パラメータ \alpha, \beta のガンマ分布を事前分布に入れます。

観測区間の最後のところ(M)には次の地点の情報がないので、一時点前の値と近い値をとるだろうという事前分布にしました。

事後分布は以下に比例します。

\left( \prod_{t=1}^{T}\prod_{m=1}^{M}( \lambda_{t,m})^{y_{t,m}}\exp (-\lambda_{t,m})\right) \\
\times \left(\prod_{t=2}^{T}\prod_{m=1}^{M-1}\lambda_{t,m}^{\lambda_{t-1,m}+\lambda_{t-1,m+1}-1}\exp(-2 \lambda_{t,m})\right) \\
\times \left(\prod_{m=1}^{M}\lambda^{\alpha-1}_{1,m}\exp(-\beta \lambda_{1,m})\right)\\
\times \left(\prod_{t=2}^{T}\lambda_{t,M}^{\lambda_{t-1,M}-1}\exp(-\lambda_{t,M})\right)

ギブスサンプリングに必要な条件付き分布は次のように求まります。

 \pi(\lambda_{t,m}) \propto \lambda_{t,m}^{y_{t,m}+\lambda_{t-1,m}+\lambda_{t-1,m+1}}\exp(-3\lambda_{t,m})
\pi(\lambda_{t,M}) \propto \lambda_{t-1,M}^{y_{t,M}-1}\exp(-2\lambda_{t-1,M})
\pi(\lambda_{1,m}) \propto \lambda_{1,m}^{y_{1,m}+\alpha-1}\exp(-(1+\beta)\lambda_{1,m})

予測

t=80 までを学習期間、80より後を予測期間として試してみました。

事後分布の平均をプロットしました。

f:id:abrahamcow:20190608155137p:plain

データに対してスムージングをしたような感じになっています。

地点ごとの予測はこんな感じです。

f:id:abrahamcow:20190608155352p:plain

右下のほうの新たな渋滞は予測できていません。

渋滞の伝播はそこそこわかるけど、渋滞の発生は突発的な事象で難しいということだと思います。

逆に言うと、予測が外れた部分はなんらかの突発的な事象(ものが落ちてるとか事故があるとか)が起きていそうなポイントであるとして、異常を発見するみたいな使い方ができるかもしれません。

コード

Rcpp です。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
List poisson_gamma_gibbs(NumericMatrix ymat, int iter, int ahead = 0) {
  int Tmax = ymat.nrow();
  int M = ymat.ncol();
  List lambdasamples(iter-1);
  NumericMatrix lambdamat(Tmax+ahead,M);
  lambdamat.fill(1);
  for(int i=1; i<iter; i++){
    lambdamat(0,0) = R::rgamma(ymat(0,0)+0.001,1/(1+0.001));
    for(int m=1; m<M; m++){
      lambdamat(0,m) = R::rgamma(ymat(0,m)+0.001,1/(1+0.001));
    }
    for (int t=1; t<Tmax; t++) {
      for(int m=0; m<M-1; m++){
        lambdamat(t,m) = R::rgamma(ymat(t,m)+lambdamat(t-1,m)+lambdamat(t-1,m+1),1.0/(2.0+1.0)); 
      }
      lambdamat(t,M-1) = R::rgamma(ymat(t,M-1)+lambdamat(t-1,M-1),0.5);
    }
    for (int t=Tmax; t<Tmax+ahead; t++) {
      for(int m=0; m<M-1; m++){
        lambdamat(t,m) = R::rgamma(lambdamat(t-1,m)+lambdamat(t-1,m+1),0.5); 
      }
      lambdamat(t,M-1) = R::rgamma(lambdamat(t-1,M-1),1.0);
    }
    NumericMatrix Lambda(Tmax+ahead,M,lambdamat.begin());
    lambdasamples[i-1] = Lambda;
  }
  return lambdasamples;
}

Rです。

library(tidyverse)
library(Rcpp)
sourceCpp("poisson_gamma_gibbs.cpp")
M<-20
Tim<-100
set.seed(1)
lam <- obs <-matrix(NA,Tim,M)
lam[1,] <-rgamma(M,1,0.01)
obs[1,] <-rpois(M,lam[1,])
for(t in 2:Tim){
  lam[t,1] <- rgamma(1,lam[t-1,1]+lam[t-1,2],2)
  obs[t,1] <- rpois(1,lam[t,1])
  for(m in 2:(M-1)){
    lam[t,m] <- rgamma(1,lam[t-1,m]+lam[t-1,m+1],2)
    obs[t,m] <- rpois(1,lam[t,m])
  } 
  lam[t,M] <- rgamma(1,lam[t-1,M],1)
  obs[t,M] <- rpois(1,lam[t,M])
}

out_gibbs <- poisson_gamma_gibbs(obs[1:80,],iter = 2001, ahead = 20)
burnin <- 1:1000
plot(simplify2array(out_gibbs)[3,3,],type="l")
lambdahat <- apply(simplify2array(out_gibbs)[,,-burnin],1:2,mean)
lambda_lower <- apply(simplify2array(out_gibbs)[,,-burnin],1:2,quantile,0.025)
lambda_upper <- apply(simplify2array(out_gibbs)[,,-burnin],1:2,quantile,0.975)

fitdf<-as.data.frame(lambdahat) %>%
  set_names(1:M) %>% 
  mutate(time=1:Tim) %>% 
  gather(pos,fit,-time) %>% 
  mutate(pos=as.integer(pos))

lowerdf<-as.data.frame(lambda_lower) %>%
  set_names(1:M) %>% 
  mutate(time=1:Tim) %>% 
  gather(pos,lower,-time) %>% 
  mutate(pos=as.integer(pos))

upperdf<-as.data.frame(lambda_upper) %>%
  set_names(1:M) %>% 
  mutate(time=1:Tim) %>% 
  gather(pos,upper,-time) %>% 
  mutate(pos=as.integer(pos))

CIdf <- left_join(lowerdf,upperdf)

obsdf<-as.data.frame(obs) %>%
  set_names(1:M) %>% 
  mutate(time=1:Tim) %>% 
  gather(pos,count,-time) %>% 
  mutate(pos=as.integer(pos),type="observed")

theme_set(theme_bw(15))

ggplot(obsdf,aes(x=pos,y=time,fill=count))+
  geom_tile()+
  scale_fill_continuous(low="white",high="red")+
  scale_y_reverse()

ggplot(fitdf,aes(x=pos,y=time,fill=fit))+
  geom_tile()+
  scale_fill_continuous(low="white",high="red")+
  scale_y_reverse()

ggplot(obsdf,aes(x=time))+
  geom_point(aes(y=count))+
  geom_line(data=fitdf,aes(y=fit),colour="royalblue")+
  geom_ribbon(data=CIdf,aes(ymin=lower,ymax=upper),fill="royalblue",linetype=2,alpha=0.5)+
  geom_vline(xintercept = 80,linetype=2)+
  facet_wrap(~pos)