統計解析R演習11

演習の概要

状態空間モデル

注意事項

慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

分析の準備

パッケージのインストール

install.packages("datasets")
install.packages("ggplot2")
install.packages("dlm")
install.packages("rstan")
install.packages("forecast")
install.packages("tibble")

library(datasets)
library(ggplot2)
library(dlm)
library(rstan)
library(forecast)
library(tibble)

オプションの設定

options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
needs::prioritize(magrittr)

データの準備

datasetsパッケージに含まれるNileデータを用います。このデータは、1871-1970年の100年間のナイル川での年間河川流量(\(10^8 m^3\))に関する時系列データです。1898年のアスワンダム建設で、1989年と1899年の間に流量が変化しています。

data(Nile)

モデル推定

最尤推定

dlmパッケージを用いて、ローカルレベルモデル、フィルタリング、スムージングなどを行えます。
まず、分析用dataframeを作成します。

Nile.df <- data.frame(
  "flow" = c(Nile),
  "year" = as.numeric(c(1871:1970), units="years")
)

次に、作成したNile.dfからナイル川の流量を可視化します。

Nile.df %>% 
  ggplot(., aes(x=year)) +
  geom_line(aes(y=flow), col="blue") +
  labs(y="Annual Flow", x ="Year")

ローカルレベルモデル

まず、ローカルレベルモデルを推定します。そのために、以下のような関数を作成します。

buildModPoly1 <- function(theta){
  dlmModPoly(order=1,           # ローカルレベルモデル。2はローカルトレンドモデル
             dV=exp(theta[1]),  # 観測方程式のノイズ
             dW=exp(theta[2]),  # 状態の分散行列の対角成分
             m0=theta[3])       # 初期値mu0
}

dlmMLE()関数を使ってモデルを推定します。

Nile.local <- dlmMLE(Nile.df$flow, parm=c(1,1,Nile.df$flow[1]), buildModPoly1)
Nile.local

$par
[1]    9.622363    7.292343 1119.999910

$value
[1] 549.63

$counts
function gradient 
      29       29 

$convergence
[1] 0

$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

推定したパラメータを使ってモデルを組み直します。

mod.Nile <- buildModPoly1(Nile.local$par)
mod.Nile

$FF
     [,1]
[1,]    1

$V
         [,1]
[1,] 15098.68

$GG
     [,1]
[1,]    1

$W
         [,1]
[1,] 1469.009

$m0
[1] 1120

$C0
      [,1]
[1,] 1e+07

フィルタリング

次に、カルマンフィルターによるフィルタリングを行います。

Nile.filt <- dlmFilter(Nile.df$flow, mod.Nile)
Nile.filt$m

  [1] 1119.9999 1120.0000 1140.9141 1072.8136 1117.3106 1129.9719 1138.4573 1048.8600 1098.0306 1171.3000 1162.9020
 [12] 1117.9508 1069.0283 1079.9775 1057.0087 1047.1237 1023.8553 1065.5552  994.3711  984.6583 1026.1416 1045.8651
 [23] 1089.6964 1105.8001 1144.3077 1175.2026 1187.1655 1145.1955 1133.1263 1037.2240  984.5566  955.0332  885.3260
 [34]  899.9263  882.0541  833.7048  855.6812  811.9712  867.5239  916.2530  930.3387  903.8110  856.3277  749.4229
 [45]  769.3382  751.3560  849.7999  916.6142  894.0185  859.2979  849.0707  827.4213  832.1156  840.6301  846.3368
 [56]  806.7244  816.9457  797.4660  797.0745  861.9462  834.4551  820.1803  832.1491  835.5808  864.5334  896.4362
 [67]  896.5868  876.6689  912.2741  874.5478  821.5269  775.4547  794.2934  799.0218  783.7945  788.3891  855.5801
 [78]  856.7604  861.3641  857.7953  866.3954  833.7105  811.0891  818.2755  880.1560  890.2620  915.8282  884.0959
 [89]  894.4850  915.9860  889.0183  923.9961  919.1904  914.3328  982.6070  963.7519  905.6027  909.1803  858.1269
[100]  819.6388  798.3718

フィルタリングの結果を可視化します(オレンジ色)。

Nile.df %>% 
  ggplot(., aes(x=year)) +
  geom_line(aes(y=flow), col = "blue") +
  geom_line(aes(y=Nile.filt$m[2:101]), col="orange") +
  labs(y="Annual Flow", x="Year")

スムージング

最後にスムージングを適用します。

Nile.smooth <- dlmSmooth(Nile.filt) # スムージングする
Nile.smooth$s

  [1] 1111.6727 1111.6715 1110.8600 1105.2674 1113.5166 1112.3785 1106.6072 1095.6411 1112.1742 1117.2436 1097.7213
 [12] 1074.0856 1058.1445 1054.1845 1044.7940 1040.3455 1037.8765 1042.9844 1034.7615 1049.4768 1073.0921 1090.1978
 [23] 1106.3497 1112.4172 1114.8280 1104.0874 1078.1776 1038.4692  999.5849  950.9308  919.4909  895.7851  874.1988
 [34]  870.1447  859.2942  851.0019  857.3038  857.8950  874.6267  877.2145  862.9914  838.4543  814.6424  799.4549
 [45]  817.6834  835.2973  865.8803  871.7390  855.3893  841.3152  834.7634  829.5507  830.3265  829.6747  825.6833
 [56]  818.1584  822.3242  824.2838  834.0545  847.5277  842.2745  845.1235  854.2112  862.2493  871.9656  874.6734
 [67]  866.7444  855.8717  848.2946  824.9844  806.9267  801.6073  811.1356  817.2719  823.9210  838.5406  856.8127
 [78]  857.2618  857.4444  856.0163  855.3682  851.3505  857.7775  874.7877  895.3773  900.9229  904.8071  900.7917
 [89]  906.8746  911.3887  909.7137  917.2538  914.7973  913.1967  912.7828  887.3432  859.5046  842.7094  818.4916
[100]  804.0510  798.3718

スムージングの結果を可視化します。(赤色)

Nile.df %>% 
  ggplot(., aes(x=year)) +
  geom_line(aes(y=flow), col="blue") +
  geom_line(aes(y=Nile.filt$m[2:101]), col="orange") +
  geom_line(aes(y=Nile.smooth$s[2:101]), col="red") +
  labs(y="Annual Flow", x="Year")

時系列モデルの推定

以下の方法により、AICを基準に最適な時系列モデルを選択することができます。ここでは、\(ARIMA(1,1,1)\)が選択されました。

model.Nile <- auto.arima(
  Nile.df$flow,
  ic="aic",
  trace=T,
  stepwise=F,
  approximation=F
)


 ARIMA(0,1,0)                    : 1296.697
 ARIMA(0,1,0) with drift         : 1298.645
 ARIMA(0,1,1)                    : 1269.091
 ARIMA(0,1,1) with drift         : 1270.309
 ARIMA(0,1,2)                    : 1267.957
 ARIMA(0,1,2) with drift         : 1268.544
 ARIMA(0,1,3)                    : 1269.208
 ARIMA(0,1,3) with drift         : 1269.604
 ARIMA(0,1,4)                    : 1270.472
 ARIMA(0,1,4) with drift         : 1270.69
 ARIMA(0,1,5)                    : 1272.462
 ARIMA(0,1,5) with drift         : 1272.689
 ARIMA(1,1,0)                    : 1281.48
 ARIMA(1,1,0) with drift         : 1283.346
 ARIMA(1,1,1)                    : 1267.255
 ARIMA(1,1,1) with drift         : 1267.637
 ARIMA(1,1,2)                    : 1268.923
 ARIMA(1,1,2) with drift         : 1269.134
 ARIMA(1,1,3)                    : 1270.858
 ARIMA(1,1,3) with drift         : 1271.085
 ARIMA(1,1,4)                    : Inf
 ARIMA(1,1,4) with drift         : 1272.689
 ARIMA(2,1,0)                    : 1277.484
 ARIMA(2,1,0) with drift         : 1279.282
 ARIMA(2,1,1)                    : 1268.896
 ARIMA(2,1,1) with drift         : 1269.132
 ARIMA(2,1,2)                    : Inf
 ARIMA(2,1,2) with drift         : Inf
 ARIMA(2,1,3)                    : Inf
 ARIMA(2,1,3) with drift         : Inf
 ARIMA(3,1,0)                    : 1278.196
 ARIMA(3,1,0) with drift         : 1279.952
 ARIMA(3,1,1)                    : 1270.873
 ARIMA(3,1,1) with drift         : 1271.057
 ARIMA(3,1,2)                    : 1272.746
 ARIMA(3,1,2) with drift         : 1272.914
 ARIMA(4,1,0)                    : 1277.333
 ARIMA(4,1,0) with drift         : 1278.989
 ARIMA(4,1,1)                    : 1272.236
 ARIMA(4,1,1) with drift         : 1272.579
 ARIMA(5,1,0)                    : 1277.125
 ARIMA(5,1,0) with drift         : 1278.674



 Best model: ARIMA(1,1,1)

summary(model.Nile)

Series: Nile.df$flow 
ARIMA(1,1,1) 

Coefficients:
         ar1      ma1
      0.2544  -0.8741
s.e.  0.1194   0.0605

sigma^2 estimated as 20177:  log likelihood=-630.63
AIC=1267.25   AICc=1267.51   BIC=1275.04

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
Training set -16.06603 139.8986 109.9998 -4.005967 12.78745 0.825499 -0.03228482

推定された\(ARIMA(1,1,1)\)を用いて50年後の将来推計を行うと、以下のようになります。

plot(forecast(model.Nile, level = c(95), h = 50))

ベイズ推定

ここでは、rstanパッケージを用いて状態空間モデルをベイズ推定します。
まず、以下のようにデータを定義します。

Nile.data.stan <- list(
  year = nrow(Nile.df),
  flow = Nile.df$flow
)

次に以下のようなstanコードを作成し、Nile.stanとして保存した後、rstanコード用の適当なフォルダに置きます。

// stan program for stated space model (local level model) for Nile data
data {
  int<lower=1> year;          // length year
  vector<lower=0>[year] flow; // Nile Annual Flow
}
parameters {
  real mu0;                   // average of state variable at t0
  vector[year] mu;            // Local level average of parameter
  real<lower=0> sigmaV;       // observed noise
  real<lower=0> sigmaW;       // Local level noise
}
model{
  // State equation
  mu[1] ~ normal(mu0, sigmaW);  
  // State transition
  for(i in 2:year) {
    mu[i] ~ normal(mu[i-1], sigmaW);
  }
  // Observed equation
  for(i in 1:year){
    flow[i] ~ normal(mu[i], sigmaV);
  }
  sigmaV ~ cauchy(0,5);
  sigmaW ~ cauchy(0,5);
}

setwd("~/(フォルダの指定)/rstan")などとして、stanコードを置いたフォルダを作業用フォルダとして指定します。
optionを指定します。

rstan_options(auto_write = TRUE)
options(mc.cores = 3)

rstan()関数を用いてstanプログラムを実行します。

Nile.stan.model1 <- rstan::stan(
  file = 'Nile.stan',
  data = Nile.data.stan,
  iter = 4500,
  warmup = 500,
  chains = 3,
  thin = 10,
  seed = 1234
)

\(Rhat\)の結果を可視化します。

stan_rhat(Nile.stan.model1)

以下の方法でtrace plotを可視化できます。

traceplot(Nile.stan.model1, pars = "mu0", inc_warmup = FALSE)
traceplot(Nile.stan.model1, pars = "sigmaV", inc_warmup = FALSE)
traceplot(Nile.stan.model1, pars = "sigmaW", inc_warmup = FALSE)

ベイズ推定の結果を格納する方法を記述します。

Nile.stan.pars1 <- rstan::extract(Nile.stan.model1)
Nile.df2 <- data.frame(
  "Year" = Nile.df$year,
  "Flow" = Nile.df$flow,
  "muM"  = apply(Nile.stan.pars1$mu, 2, mean),
  "muL"  = apply(Nile.stan.pars1$mu, 2, quantile, 0.025),
  "muU"  = apply(Nile.stan.pars1$mu, 2, quantile, 0.925)
)

推定結果を可視化します。

Nile.df2 %>%
  ggplot(., aes(x=Year, y=Flow)) +
  geom_point(col = "black") +
  labs(y = "Nile Flow", x = "Year") + 
  geom_smooth(aes(ymin=muL, ymax=muU, y=muM),
              lwd=1.2,color="blue",stat="identity"
  )

参考資料

多くの資料を参考にさせていただきました。ありがとうございました。
資料1
資料2
資料3
資料4
資料5

LS0tCnRpdGxlOiAi57Wx6KiI6Kej5p6QUua8lOe/kjExIgojYXV0aG9yOiAiVG9tb3l1a2kgRnVydXRhbmkiCiNkYXRlOiAiMjAxOS84LzI2IgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIGh0bWxfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQprbml0cjo6b3B0c19jaHVuayRzZXQoY2FjaGUgPSBUUlVFKQprbml0cjo6b3B0c19jaHVuayRzZXQoZXJyb3IgPSBUUlVFKQpgYGAKCiMjIOa8lOe/kuOBruamguimgQotIOeKtuaFi+epuumWk+ODouODh+ODqwoKIyMjIOazqOaEj+S6i+mghQotIOaFtuaHiee+qeWhvuWkp+WtplNGQ+OBp+mWi+ism+OBl+OBpuOBhOOCi+OAjOe1seioiOino+aekOOAjeOBruaOiOalreWxpeS/ruiAheWQkeOBkeOBrua8lOe/kueUqOODmuODvOOCuOOBp+OBmeOAggotIOW/heOBmuOBl+OCguWFqOOBpuOBruODkOODvOOCuOODp+ODs+OBrlLjgoRPU+OBp+WLleS9nOeiuuiqjeOCkuihjOOBo+OBpuOBhOOBvuOBm+OCk+OAguOBk+OBrua8lOe/kueUqOODmuODvOOCuOOCkuS9nOaIkOOBl+OBpuOBhOOCi+autemajuOBp+OBr+OAgVIzLjYuMOOCkuS9v+OBo+OBpuOBhOOBvuOBmeOAggotIFLjga7mm7TmlrDjgarjganjgavjgojjgorjgIFS44Kz44O844OJ44Gq44Gp44KS5LqI5ZGK54Sh44GX44Gr5aSJ5pu044GZ44KL5aC05ZCI44GM44GC44KK44G+44GZ44CCIAoKIyMg5YiG5p6Q44Gu5rqW5YKZCiMjIyDjg5Hjg4PjgrHjg7zjgrjjga7jgqTjg7Pjgrnjg4jjg7zjg6sKYGBge3IgaW5zdGFsbCBwYWNrYWdlcywgZXZhbD1GQUxTRX0KaW5zdGFsbC5wYWNrYWdlcygiZGF0YXNldHMiKQppbnN0YWxsLnBhY2thZ2VzKCJnZ3Bsb3QyIikKaW5zdGFsbC5wYWNrYWdlcygiZGxtIikKaW5zdGFsbC5wYWNrYWdlcygicnN0YW4iKQppbnN0YWxsLnBhY2thZ2VzKCJmb3JlY2FzdCIpCmluc3RhbGwucGFja2FnZXMoInRpYmJsZSIpCmBgYAoKYGBge3IgbGlicmFyeSwgZXZhbD1GQUxTRX0KbGlicmFyeShkYXRhc2V0cykKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGRsbSkKbGlicmFyeShyc3RhbikKbGlicmFyeShmb3JlY2FzdCkKbGlicmFyeSh0aWJibGUpCmBgYAoKIyMjIOOCquODl+OCt+ODp+ODs+OBruioreWumgpgYGB7ciBvcHRpb25zfQpvcHRpb25zKG1jLmNvcmVzID0gcGFyYWxsZWw6OmRldGVjdENvcmVzKCkpCnJzdGFuX29wdGlvbnMoYXV0b193cml0ZSA9IFRSVUUpCm5lZWRzOjpwcmlvcml0aXplKG1hZ3JpdHRyKSAKYGBgCgojIyMg44OH44O844K/44Gu5rqW5YKZCi0gYGRhdGFzZXRzYOODkeODg+OCseODvOOCuOOBq+WQq+OBvuOCjOOCi2BOaWxlYOODh+ODvOOCv+OCkueUqOOBhOOBvuOBmeOAguOBk+OBruODh+ODvOOCv+OBr+OAgTE4NzEtMTk3MOW5tOOBrjEwMOW5tOmWk+OBruODiuOCpOODq+W3neOBp+OBruW5tOmWk+ays+W3nea1gemHjygkMTBeOCBtXjMkKeOBq+mWouOBmeOCi+aZguezu+WIl+ODh+ODvOOCv+OBp+OBmeOAgjE4OTjlubTjga7jgqLjgrnjg6/jg7Pjg4Djg6Dlu7roqK3jgafjgIExOTg55bm044GoMTg5OeW5tOOBrumWk+OBq+a1gemHj+OBjOWkieWMluOBl+OBpuOBhOOBvuOBmeOAggoKYGBge3IgZGF0YU5pbGV9CmRhdGEoTmlsZSkKYGBgCgojIyDjg6Ljg4fjg6vmjqjlrpoKIyMjIOacgOWwpOaOqOWumgotIGBkbG1g44OR44OD44Kx44O844K444KS55So44GE44Gm44CB44Ot44O844Kr44Or44Os44OZ44Or44Oi44OH44Or44CB44OV44Kj44Or44K/44Oq44Oz44Kw44CB44K544Og44O844K444Oz44Kw44Gq44Gp44KS6KGM44GI44G+44GZ44CCCi0g44G+44Ga44CB5YiG5p6Q55SoZGF0YWZyYW1l44KS5L2c5oiQ44GX44G+44GZ44CCCmBgYHtyIE5pbGUuZGZ9Ck5pbGUuZGYgPC0gZGF0YS5mcmFtZSgKICAiZmxvdyIgPSBjKE5pbGUpLAogICJ5ZWFyIiA9IGFzLm51bWVyaWMoYygxODcxOjE5NzApLCB1bml0cz0ieWVhcnMiKQopCmBgYAotIOasoeOBq+OAgeS9nOaIkOOBl+OBn2BOaWxlLmRmYOOBi+OCieODiuOCpOODq+W3neOBrua1gemHj+OCkuWPr+imluWMluOBl+OBvuOBmeOAggpgYGB7ciBOaWxlLnBsb3QxfQpOaWxlLmRmICU+JSAKICBnZ3Bsb3QoLiwgYWVzKHg9eWVhcikpICsKICBnZW9tX2xpbmUoYWVzKHk9ZmxvdyksIGNvbD0iYmx1ZSIpICsKICBsYWJzKHk9IkFubnVhbCBGbG93IiwgeCA9IlllYXIiKSAKYGBgCiMjIyMg44Ot44O844Kr44Or44Os44OZ44Or44Oi44OH44OrCi0g44G+44Ga44CB44Ot44O844Kr44Or44Os44OZ44Or44Oi44OH44Or44KS5o6o5a6a44GX44G+44GZ44CC44Gd44Gu44Gf44KB44Gr44CB5Lul5LiL44Gu44KI44GG44Gq6Zai5pWw44KS5L2c5oiQ44GX44G+44GZ44CCCmBgYHtyIGJ1aWxkTW9kUG9seTF9CmJ1aWxkTW9kUG9seTEgPC0gZnVuY3Rpb24odGhldGEpewogIGRsbU1vZFBvbHkob3JkZXI9MSwgICAgICAgICAgICMg44Ot44O844Kr44Or44Os44OZ44Or44Oi44OH44Or44CCMuOBr+ODreODvOOCq+ODq+ODiOODrOODs+ODieODouODh+ODqwogICAgICAgICAgICAgZFY9ZXhwKHRoZXRhWzFdKSwgICMg6Kaz5ris5pa556iL5byP44Gu44OO44Kk44K6CiAgICAgICAgICAgICBkVz1leHAodGhldGFbMl0pLCAgIyDnirbmhYvjga7liIbmlaPooYzliJfjga7lr77op5LmiJDliIYKICAgICAgICAgICAgIG0wPXRoZXRhWzNdKSAgICAgICAjIOWIneacn+WApG11MAp9CmBgYAotIGBkbG1NTEUoKWDplqLmlbDjgpLkvb/jgaPjgabjg6Ljg4fjg6vjgpLmjqjlrprjgZfjgb7jgZnjgIIKYGBge3IgTmlsZS5sb2NhbH0KTmlsZS5sb2NhbCA8LSBkbG1NTEUoTmlsZS5kZiRmbG93LCBwYXJtPWMoMSwxLE5pbGUuZGYkZmxvd1sxXSksIGJ1aWxkTW9kUG9seTEpCk5pbGUubG9jYWwKYGBgCi0g5o6o5a6a44GX44Gf44OR44Op44Oh44O844K/44KS5L2/44Gj44Gm44Oi44OH44Or44KS57WE44G/55u044GX44G+44GZ44CCCmBgYHtyIG1vZC5OaWxlfQptb2QuTmlsZSA8LSBidWlsZE1vZFBvbHkxKE5pbGUubG9jYWwkcGFyKQptb2QuTmlsZQpgYGAKIyMjIyDjg5XjgqPjg6vjgr/jg6rjg7PjgrAKLSDmrKHjgavjgIHjgqvjg6vjg57jg7Pjg5XjgqPjg6vjgr/jg7zjgavjgojjgovjg5XjgqPjg6vjgr/jg6rjg7PjgrDjgpLooYzjgYTjgb7jgZnjgIIKYGBge3IgTmlsZS5maWx0fQpOaWxlLmZpbHQgPC0gZGxtRmlsdGVyKE5pbGUuZGYkZmxvdywgbW9kLk5pbGUpCk5pbGUuZmlsdCRtCmBgYAotIOODleOCo+ODq+OCv+ODquODs+OCsOOBrue1kOaenOOCkuWPr+imluWMluOBl+OBvuOBmSjjgqrjg6zjg7PjgrjoibIp44CCCmBgYHtyIE5pbGUucGxvdDJ9Ck5pbGUuZGYgJT4lIAogIGdncGxvdCguLCBhZXMoeD15ZWFyKSkgKwogIGdlb21fbGluZShhZXMoeT1mbG93KSwgY29sPSJibHVlIikgKwogIGdlb21fbGluZShhZXMoeT1OaWxlLmZpbHQkbVsyOjEwMV0pLCBjb2w9Im9yYW5nZSIpICsKICBsYWJzKHk9IkFubnVhbCBGbG93IiwgeD0iWWVhciIpIAoKYGBgCiMjIyMg44K544Og44O844K444Oz44KwCi0g5pyA5b6M44Gr44K544Og44O844K444Oz44Kw44KS6YGp55So44GX44G+44GZ44CCCmBgYHtyIE5pbGUuc21vb3RofQpOaWxlLnNtb290aCA8LSBkbG1TbW9vdGgoTmlsZS5maWx0KSAjIOOCueODoOODvOOCuOODs+OCsOOBmeOCiwpOaWxlLnNtb290aCRzCmBgYAotIOOCueODoOODvOOCuOODs+OCsOOBrue1kOaenOOCkuWPr+imluWMluOBl+OBvuOBmeOAgijotaToibIpCmBgYHtyIE5pbGUucGxvdDN9Ck5pbGUuZGYgJT4lIAogIGdncGxvdCguLCBhZXMoeD15ZWFyKSkgKwogIGdlb21fbGluZShhZXMoeT1mbG93KSwgY29sPSJibHVlIikgKwogIGdlb21fbGluZShhZXMoeT1OaWxlLmZpbHQkbVsyOjEwMV0pLCBjb2w9Im9yYW5nZSIpICsKICBnZW9tX2xpbmUoYWVzKHk9TmlsZS5zbW9vdGgkc1syOjEwMV0pLCBjb2w9InJlZCIpICsKICBsYWJzKHk9IkFubnVhbCBGbG93IiwgeD0iWWVhciIpIAoKYGBgCiMjIyMg5pmC57O75YiX44Oi44OH44Or44Gu5o6o5a6aCi0g5Lul5LiL44Gu5pa55rOV44Gr44KI44KK44CBQUlD44KS5Z+65rqW44Gr5pyA6YGp44Gq5pmC57O75YiX44Oi44OH44Or44KS6YG45oqe44GZ44KL44GT44Go44GM44Gn44GN44G+44GZ44CC44GT44GT44Gn44Gv44CBJEFSSU1BKDEsMSwxKSTjgYzpgbjmip7jgZXjgozjgb7jgZfjgZ/jgIIKYGBge3J9Cm1vZGVsLk5pbGUgPC0gYXV0by5hcmltYSgKICBOaWxlLmRmJGZsb3csCiAgaWM9ImFpYyIsCiAgdHJhY2U9VCwKICBzdGVwd2lzZT1GLAogIGFwcHJveGltYXRpb249RgopCmBgYApgYGB7ciBzdW1tYXJ5Lm1vZGVsLk5pbGV9CnN1bW1hcnkobW9kZWwuTmlsZSkKYGBgCi0g5o6o5a6a44GV44KM44GfJEFSSU1BKDEsMSwxKSTjgpLnlKjjgYTjgaY1MOW5tOW+jOOBruWwhuadpeaOqOioiOOCkuihjOOBhuOBqOOAgeS7peS4i+OBruOCiOOBhuOBq+OBquOCiuOBvuOBmeOAggpgYGB7ciBmb3JlY2FzdC5tb2RlbC5OaWxlfQpwbG90KGZvcmVjYXN0KG1vZGVsLk5pbGUsIGxldmVsID0gYyg5NSksIGggPSA1MCkpCmBgYAoKIyMjIOODmeOCpOOCuuaOqOWumgotIOOBk+OBk+OBp+OBr+OAgWByc3RhbmDjg5Hjg4PjgrHjg7zjgrjjgpLnlKjjgYTjgabnirbmhYvnqbrplpPjg6Ljg4fjg6vjgpLjg5njgqTjgrrmjqjlrprjgZfjgb7jgZnjgIIKLSDjgb7jgZrjgIHku6XkuIvjga7jgojjgYbjgavjg4fjg7zjgr/jgpLlrprnvqnjgZfjgb7jgZnjgIIKYGBge3IgTmlsZS5kYXRhLnN0YW59Ck5pbGUuZGF0YS5zdGFuIDwtIGxpc3QoCiAgeWVhciA9IG5yb3coTmlsZS5kZiksCiAgZmxvdyA9IE5pbGUuZGYkZmxvdwopCmBgYAotIOasoeOBq+S7peS4i+OBruOCiOOBhuOBqmBzdGFuYOOCs+ODvOODieOCkuS9nOaIkOOBl+OAgWBOaWxlLnN0YW5g44Go44GX44Gm5L+d5a2Y44GX44Gf5b6M44CBcnN0YW7jgrPjg7zjg4nnlKjjga7pganlvZPjgarjg5Xjgqnjg6vjg4Djgavnva7jgY3jgb7jgZnjgIIKYGBge3N0YW59Ci8vIHN0YW4gcHJvZ3JhbSBmb3Igc3RhdGVkIHNwYWNlIG1vZGVsIChsb2NhbCBsZXZlbCBtb2RlbCkgZm9yIE5pbGUgZGF0YQpkYXRhIHsKICBpbnQ8bG93ZXI9MT4geWVhcjsgICAgICAgICAgLy8gbGVuZ3RoIHllYXIKICB2ZWN0b3I8bG93ZXI9MD5beWVhcl0gZmxvdzsgLy8gTmlsZSBBbm51YWwgRmxvdwp9CnBhcmFtZXRlcnMgewogIHJlYWwgbXUwOyAgICAgICAgICAgICAgICAgICAvLyBhdmVyYWdlIG9mIHN0YXRlIHZhcmlhYmxlIGF0IHQwCiAgdmVjdG9yW3llYXJdIG11OyAgICAgICAgICAgIC8vIExvY2FsIGxldmVsIGF2ZXJhZ2Ugb2YgcGFyYW1ldGVyCiAgcmVhbDxsb3dlcj0wPiBzaWdtYVY7ICAgICAgIC8vIG9ic2VydmVkIG5vaXNlCiAgcmVhbDxsb3dlcj0wPiBzaWdtYVc7ICAgICAgIC8vIExvY2FsIGxldmVsIG5vaXNlCn0KbW9kZWx7CiAgLy8gU3RhdGUgZXF1YXRpb24KICBtdVsxXSB+IG5vcm1hbChtdTAsIHNpZ21hVyk7ICAKICAvLyBTdGF0ZSB0cmFuc2l0aW9uCiAgZm9yKGkgaW4gMjp5ZWFyKSB7CiAgICBtdVtpXSB+IG5vcm1hbChtdVtpLTFdLCBzaWdtYVcpOwogIH0KICAvLyBPYnNlcnZlZCBlcXVhdGlvbgogIGZvcihpIGluIDE6eWVhcil7CiAgICBmbG93W2ldIH4gbm9ybWFsKG11W2ldLCBzaWdtYVYpOwogIH0KICBzaWdtYVYgfiBjYXVjaHkoMCw1KTsKICBzaWdtYVcgfiBjYXVjaHkoMCw1KTsKfQpgYGAKLSBgc2V0d2QoIn4vKOODleOCqeODq+ODgOOBruaMh+WumikvcnN0YW4iKWDjgarjganjgajjgZfjgabjgIFgc3RhbmDjgrPjg7zjg4njgpLnva7jgYTjgZ/jg5Xjgqnjg6vjg4DjgpLkvZzmpa3nlKjjg5Xjgqnjg6vjg4DjgajjgZfjgabmjIflrprjgZfjgb7jgZnjgIIKLSBgb3B0aW9uYOOCkuaMh+WumuOBl+OBvuOBmeOAggpgYGB7ciByc3Rhbl9vcHRpb25zfQpyc3Rhbl9vcHRpb25zKGF1dG9fd3JpdGUgPSBUUlVFKQpvcHRpb25zKG1jLmNvcmVzID0gMykKYGBgCi0gYHJzdGFuKClg6Zai5pWw44KS55So44GE44GmYHN0YW5g44OX44Ot44Kw44Op44Og44KS5a6f6KGM44GX44G+44GZ44CCCmBgYHty44CATmlsZS5zdGFuLm1vZGVsMX0KTmlsZS5zdGFuLm1vZGVsMSA8LSByc3Rhbjo6c3RhbigKICBmaWxlID0gJ05pbGUuc3RhbicsCiAgZGF0YSA9IE5pbGUuZGF0YS5zdGFuLAogIGl0ZXIgPSA0NTAwLAogIHdhcm11cCA9IDUwMCwKICBjaGFpbnMgPSAzLAogIHRoaW4gPSAxMCwKICBzZWVkID0gMTIzNAopCmBgYAotICRSaGF0JOOBrue1kOaenOOCkuWPr+imluWMluOBl+OBvuOBmeOAggpgYGB7cn0Kc3Rhbl9yaGF0KE5pbGUuc3Rhbi5tb2RlbDEpCmBgYAotIOS7peS4i+OBruaWueazleOBp3RyYWNlIHBsb3TjgpLlj6/oppbljJbjgafjgY3jgb7jgZnjgIIKYGBge3IgdHJhY2VwbG90Lk5pbGUuc3Rhbi5tb2RlbDF9CnRyYWNlcGxvdChOaWxlLnN0YW4ubW9kZWwxLCBwYXJzID0gIm11MCIsIGluY193YXJtdXAgPSBGQUxTRSkKdHJhY2VwbG90KE5pbGUuc3Rhbi5tb2RlbDEsIHBhcnMgPSAic2lnbWFWIiwgaW5jX3dhcm11cCA9IEZBTFNFKQp0cmFjZXBsb3QoTmlsZS5zdGFuLm1vZGVsMSwgcGFycyA9ICJzaWdtYVciLCBpbmNfd2FybXVwID0gRkFMU0UpCmBgYAoKLSDjg5njgqTjgrrmjqjlrprjga7ntZDmnpzjgpLmoLzntI3jgZnjgovmlrnms5XjgpLoqJjov7DjgZfjgb7jgZnjgIIKYGBge3LjgIBOaWxlLnN0YW4ucGFyczF9Ck5pbGUuc3Rhbi5wYXJzMSA8LSByc3Rhbjo6ZXh0cmFjdChOaWxlLnN0YW4ubW9kZWwxKQpOaWxlLmRmMiA8LSBkYXRhLmZyYW1lKAogICJZZWFyIiA9IE5pbGUuZGYkeWVhciwKICAiRmxvdyIgPSBOaWxlLmRmJGZsb3csCiAgIm11TSIgID0gYXBwbHkoTmlsZS5zdGFuLnBhcnMxJG11LCAyLCBtZWFuKSwKICAibXVMIiAgPSBhcHBseShOaWxlLnN0YW4ucGFyczEkbXUsIDIsIHF1YW50aWxlLCAwLjAyNSksCiAgIm11VSIgID0gYXBwbHkoTmlsZS5zdGFuLnBhcnMxJG11LCAyLCBxdWFudGlsZSwgMC45MjUpCikKYGBgCgotIOaOqOWumue1kOaenOOCkuWPr+imluWMluOBl+OBvuOBmeOAggpgYGB7ciBOaWxlLnBsb3Q0fQpOaWxlLmRmMiAlPiUKICBnZ3Bsb3QoLiwgYWVzKHg9WWVhciwgeT1GbG93KSkgKwogIGdlb21fcG9pbnQoY29sID0gImJsYWNrIikgKwogIGxhYnMoeSA9ICJOaWxlIEZsb3ciLCB4ID0gIlllYXIiKSArIAogIGdlb21fc21vb3RoKGFlcyh5bWluPW11TCwgeW1heD1tdVUsIHk9bXVNKSwKICAgICAgICAgICAgICBsd2Q9MS4yLGNvbG9yPSJibHVlIixzdGF0PSJpZGVudGl0eSIKICApCmBgYAoKCiMjIOWPguiAg+izh+aWmQotIOWkmuOBj+OBruizh+aWmeOCkuWPguiAg+OBq+OBleOBm+OBpuOBhOOBn+OBoOOBjeOBvuOBl+OBn+OAguOBguOCiuOBjOOBqOOBhuOBlOOBluOBhOOBvuOBl+OBn+OAggotIFvos4fmlpkxXShodHRwOi8vcXVhbnRkZXZlbC5jb20vcHVibGljL1N0YXRlU3BhY2VNb2RlbHMvaW5kZXgucGRmKQotIFvos4fmlpkyXShodHRwczovL3d3dy5qc3RhdHNvZnQub3JnL2FydGljbGUvdmlldy92MDQxaTA0L3Y0MWkwNC5wZGYpCi0gW+izh+aWmTNdKGh0dHBzOi8vbG9naWNzLW9mLWJsdWUuY29tL2RsbSVlMyU4MSVhZSVlNCViZCViZiVlMyU4MSU4NCVlNiU5NiViOS8pCi0gW+izh+aWmTRdKGh0dHBzOi8vc2h1bWFnaXQuZ2l0aHViLmlvL215YmxvZy8yMDE4LzA4LzA5L3N0YW4tbG9jYWwtbGV2ZWwtd2l0aC1hci8pCi0gW+izh+aWmTVdKGh0dHA6Ly93d3cua29zdWdpdHRpLm5ldC9hcmNoaXZlcy81Nzg2KQ==