練習参加者数の予測
UQ_baritone (71期)
2024-08-04
2. 実はよくわかってないので素人質問されるときつい
必要なpackageの追加
# 使用するpackageの追加
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(tseries)
library(TSSS)## Registered S3 method overwritten by 'TSSS':
## method from
## plot.arma tserieslibrary(KFAS)## Please cite KFAS in publications by using:
##
## Jouni Helske (2017). KFAS: Exponential Family State Space Models in R. Journal of Statistical Software, 78(10), 1-39. doi:10.18637/jss.v078.i10.library(ggplot2)2.1 データの概要と予測の方針
当然だが，私は緑会合唱団に所属している．基本的に週二(月曜日と木曜日)で練習を行っていて，毎回練習に参加する者を記録している．記録を遡ることができた2023年6月27日から現在に至るまでの練習参加人数を時系列データとして扱い，将来の練習参加人数の予測をこのレポートで行う． また，夏と春に合宿を行っている．合宿期間中における月曜日と木曜日の参加人数は合宿参加人数で代用している．
予測の方針としては，
ARモデルによる予測
状態空間モデルによる予測
説明変数を複数加えた予測(実際の予測はAR modelをベースに行う)
の3つを試みる．
2024年6月までのデータを用いて直近一ヶ月(7月)の練習参加者数を予測することでモデルの予測性能を比較することは候補として挙がったが，7月は
発表会が終わった直後であること
テスト期間であること
主に2つの理由により，練習参加者が少ないのである．その特殊性ゆえに，直近一ヶ月は予測性能を測ることは適切であるとは言えないため，既存のデータを用いて予測性能を測ることはしなかった．
今回取り扱う時系列について図示すると以下のようになる．
# データの読み込み
data_mcc <- read.csv("/Users/suzukikouhei/Library/Mobile Documents/com~apple~CloudDocs/10_授業/3S/3_5_数理手法7/data分析/data.csv",
header = FALSE)
#xtsデータとして変換
xts_sample <- as.xts(
matrix(
data_mcc[,2],
dimnames = list(
X = data_mcc[,1],
Y = "participants")
)
)
#図示
plot(xts_sample[,'participants'],
ylab = 'participants',
main = 'MCC participants from:2023-06-29 to:2024-07-29')
2.2 データを観察する
弱定常であるか確認するために，ADF検定を実行すると
p-value は
0.09542となり，定常であることは棄却される．
#単位根検定
adf.test(xts_sample)##
## Augmented Dickey-Fuller Test
##
## data: xts_sample
## Dickey-Fuller = -3.1755, Lag order = 4, p-value = 0.09542
## alternative hypothesis: stationaryforecast package の
ggtsdisplay()関数を用いると，
時系列
自己相関係数(ACF)
偏自己相関係数(P-ACF)
をみることができる．
ggtsdisplay(xts_sample[,'participants'], main = 'MCC participants from:2023-06-29 to:2024-07-29' )
ここで，lagが5とその整数倍付近(11, 15, 20,
21)のとき，自己相関が大きくなっていることがわかる．これは周期を5とした周期性を示唆するがそれの原因となる事象は思い浮かばない．
(練習参加のモチベーションが5回の周期性で変化しているのであろうか？)
また，lagが9, 10のときも若干ピークが観察される，これは月に練習がおよそ9,10回開催されるため，一ヶ月を周期とした周期性が存在することが示唆される．
このデータには月曜日と木曜日の練習以外にも他の曜日に臨時で開催される練習日も含んでいる．臨時練習は特別な練習メニューであることが多いため，月木の練習日だけを抽出して自己相関を観察してみる．
# 読み込み
mon_thu <- read.csv("/Users/suzukikouhei/Library/Mobile Documents/com~apple~CloudDocs/10_授業/3S/3_5_数理手法7/data分析/mon_thu_only.csv")
# xtsへの変換
mon_thu.xts <- as.xts(
matrix(
mon_thu[,2],
dimnames = list(
X = mon_thu[,1],
Y = "participants")
)
)
ggtsdisplay(mon_thu.xts[,'participants'], main = 'MCC participants from:2023-06-29 to:2024-07-29' )
先ほど見られた，lagが5のときのピークは弱まっているが存在する．どのような練習形態に起因するのかは私には不明である．
2.3 ARモデルによる予測
これからは時系列のデータのみからAR modelを用いて今後20日間の練習参加人数を予測する．
3つの手法を試している．
何も手を加えずにARmodelをつかった予測
月曜日と木曜日の練習だけを抽出してARmodelを当てはめた予測
前週同曜日の練習からの差分でAR modelを当てはめた予測
2.3.1 AR modelをつかった予測
弱定常であることは棄却されたが，平均と分散についての変化は見て取れないので，AR modelへの当てはめを試みる．
AICによると，AR
modelの係数が5のときが最も適当である．
arfit(xts_sample, lag=10, method=1 ) # method = 1のときユールウォーカー法を使用する
20日間分の予測は以下のようになる．
# 20日分予測したときにそれに対応する日付
data_20 <-
c("2024-08-01",
"2024-08-05",
"2024-08-08",
"2024-08-12",
"2024-08-15",
"2024-08-19",
"2024-08-22",
"2024-08-26",
"2024-08-29",
"2024-09-02",
"2024-09-05",
"2024-09-09",
"2024-09-12",
"2024-09-16",
"2024-09-19",
"2024-09-23",
"2024-09-26",
"2024-09-30",
"2024-10-03",
"2024-10-07"
)
# AR modelをもとに予測してみる
MCC.ar <- ar(xts_sample)
MCC.pr <- predict(MCC.ar, n.ahead=20)
SE1 <- MCC.pr$pred+ 1*MCC.pr$se
SE2 <- MCC.pr$pred- 1*MCC.pr$se
#今後20日間の予測
pred_future20 <- as.xts(matrix(
MCC.pr$pred[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
se1_future20 <- as.xts(
matrix(
SE1[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
se2_future20 <- as.xts(
matrix(
SE2[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
combined_xts <- merge(xts_sample, pred_future20 , se1_future20, se2_future20, all = TRUE)
plot.zoo(combined_xts, screens = 1,col = c(1, 2, 4 ,4), lty = c(1, 1, 2,2),xlab = "Date", ylab = "participants", main = "Time Series(AR model)")
# 月単位の軸を設定
axis(1, at = seq(start(combined_xts), end(combined_xts), by = "months"), labels = format(seq(start(combined_xts), end(combined_xts), by = "months"), "%Y-%m"), las = 2, cex.axis = 0.7)
legend("topleft", col = c(1, 2, 3,3),lty = c(1, 1, 2,2),legend = c("real", "pred", "pred + se1","pred - se1"))
2.3.2 月曜日と木曜日だけのデータで予測
月木の練習だけのデータ方が
週ごとの周期性が見やすい
特異的な練習を除くことができる
という2つのメリットがあると考えられ，予測精度が高まる可能性がある．
先ほどの月曜日と木曜日だけのデータを用いる．
定常性を確認するため，ADF検定を実行する．p-value = 0.1357となるため定常性は棄却された．
adf.test(mon_thu.xts)##
## Augmented Dickey-Fuller Test
##
## data: mon_thu.xts
## Dickey-Fuller = -3.0632, Lag order = 4, p-value = 0.1357
## alternative hypothesis: stationaryこれも弱定常であることは棄却されたが，平均と分散についての変化は目では見て取れないので，AR modelへの当てはめを試みる． AICを最小化する係数数を調べると，AR modelに当てはめたときの最適な係数の数は3つである．
arfit(mon_thu.xts, lag=6, method=1 )
これで今後20日間分の予測をしてみると，
# 予測図をつくる
mon_thu.ar <- ar(mon_thu.xts)
mon_thu.pr <- predict(mon_thu.ar, n.ahead=20)
mon_thu.SE1 <- mon_thu.pr$pred+ 1*mon_thu.pr$se
mon_thu.SE2 <- mon_thu.pr$pred- 1*mon_thu.pr$se
mon_thu.pred_20 <- as.xts(matrix(
mon_thu.pr$pred[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
mon_thu.se1_20 <- as.xts(
matrix(
mon_thu.SE1[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
mon_thu.se2_20 <- as.xts(
matrix(
mon_thu.SE2[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
combined_xts <- merge(mon_thu.xts,mon_thu.pred_20 , mon_thu.se1_20, mon_thu.se2_20, all = TRUE)
plot.zoo(combined_xts, screens = 1,col = c(1, 2, 4 ,4), lty = c(1, 1, 2,2),xlab = "Date", ylab = "Values", main = "participants on Monday and Thursday")
# 月単位の軸を設定
axis(1, at = seq(start(combined_xts), end(combined_xts), by = "months"), labels = format(seq(start(combined_xts), end(combined_xts), by = "months"), "%Y-%m"), las = 2, cex.axis = 0.7)
legend("topright", col = c(1, 2, 3,3),lty = c(1, 1, 2,2),legend = c("real", "pred", "pred + se1","pred - se1"))
2.3.3 前週同曜日との差分で予測
lag = 2のときの自己相関が大きいこと，実際に授業等により月曜日だけ(または木曜日だけ)練習に来る人も一定数存在すること,を踏まえて，前週同曜日からの差分を時系列データとして解析する． プロットしたところ分散の時間的な変化は見受けられるが，概ね定常であると考えられる．
#前週同曜日からの増減
week_ratio <- (mon_thu.xts - lag(mon_thu.xts,k = 2))[3:length(mon_thu.xts)]
#プロット
plot(week_ratio)
ADF検定を実行すると，p-value = 0.01となり定常性は採択される．
adf.test(week_ratio)## Warning in adf.test(week_ratio): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: week_ratio
## Dickey-Fuller = -5.192, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryarfit( week_ratio, lag=15, method=1 )
これで今後20日間分の予測をしてみると，
#予測図を作成
week_ratio.ar <- ar(week_ratio)
week_ratio.pr <- predict(week_ratio.ar, n.ahead=20)
week_ratio.SE1 <- week_ratio.pr$pred+ 1*week_ratio.pr$se
week_ratio.SE2 <- week_ratio.pr$pred- 1*week_ratio.pr$se
week_ratio.pred_20 <- as.xts(matrix(
week_ratio.pr$pred[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
week_ratio.se1_20 <- as.xts(
matrix(
week_ratio.SE1[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
week_ratio.se2_20 <- as.xts(
matrix(
week_ratio.SE2[1:20],
dimnames = list(
X = data_20,
Y = "participants")
)
)
combined_xts <- merge(week_ratio,week_ratio.pred_20 , week_ratio.se1_20, week_ratio.se2_20, all = TRUE)
plot.zoo(combined_xts, screens = 1,col = c(1, 2, 4 ,4), lty = c(1, 1, 2,2),xlab = "Date", ylab = "Values", main = "The difference in the number of participants from the same day last week")
# 月単位の軸を設定
axis(1, at = seq(start(combined_xts), end(combined_xts), by = "months"), labels = format(seq(start(combined_xts), end(combined_xts), by = "months"), "%Y-%m"), las = 2, cex.axis = 0.7)
legend("topright", col = c(1, 2, 3,3),lty = c(1, 1, 2,2),legend = c("real", "pred", "pred + se1","pred - se1"))
人数ベースの予測に直す．
week_number.pr <- 1:20
week_number.SE1 <- 1:20
week_number.SE2 <- 1:20
week_number.pr[1] <- mon_thu.xts["2024-07-25"] + week_ratio.pr$pred[1]
week_number.pr[2] <- mon_thu.xts["2024-07-29"] + week_ratio.pr$pred[2]
week_number.SE1[1] <- mon_thu.xts["2024-07-25"] + week_ratio.SE1[1]
week_number.SE1[2] <- mon_thu.xts["2024-07-29"] + week_ratio.SE1[2]
week_number.SE2[1] <- mon_thu.xts["2024-07-25"] + week_ratio.SE2[1]
week_number.SE2[2] <- mon_thu.xts["2024-07-29"] + week_ratio.SE2[2]
for (i in 3:20){
week_number.pr[i] <- week_number.pr[i-2] + week_ratio.pr$pred[i]
week_number.SE1[i] <- week_number.pr[i-2] + week_ratio.SE1[i]
week_number.SE2[i] <- week_number.pr[i-2] + week_ratio.SE2[i]
}
week_number.pred_20 <- as.xts(matrix(
week_number.pr,
dimnames = list(
X = data_20,
Y = "participants")
)
)
week_number.se1_20 <- as.xts(
matrix(
week_number.SE1,
dimnames = list(
X = data_20,
Y = "participants")
)
)
week_number.se2_20 <- as.xts(
matrix(
week_number.SE2,
dimnames = list(
X = data_20,
Y = "participants")
)
)
combined_xts <- merge(mon_thu.xts,week_number.pred_20 , week_number.se1_20, week_number.se2_20, all = TRUE)
plot.zoo(combined_xts, screens = 1,col = c(1, 2, 4 ,4), lty = c(1, 1, 2,2),xlab = "Date", ylab = "participants", main = "Time Series predicted by week trend")
# 月単位の軸を設定
axis(1, at = seq(start(combined_xts), end(combined_xts), by = "months"), labels = format(seq(start(combined_xts), end(combined_xts), by = "months"), "%Y-%m"), las = 2, cex.axis = 0.7)
legend("topleft", col = c(1, 2, 3,3),lty = c(1, 1, 2, 2),legend = c("real", "pred", "pred + se0.5","pred - se0.5"))
2.3.4 3つのモデルを同一時間軸上で表示
ARモデルによる3つの予測を同日時間軸上で表示する． 3つ目の予測が一練習日ごとに変化が激しいのは，前週からの差分を予測しているため，直近２練習日の結果に大きく依存してしまうからである．
combined_3model <- merge(xts_sample, pred_future20 ,mon_thu.pred_20 , week_number.pred_20, all = TRUE)
plot.zoo(combined_3model, screens = 1,col = c(1, 2, 3 ,4), lty = c(1, 1, 1,1),xlab = "Date", ylab = "participants", main = "Three Time Series predicted by week trend")
# 月単位の軸を設定
axis(1, at = seq(start(combined_xts), end(combined_xts), by = "months"), labels = format(seq(start(combined_xts), end(combined_xts), by = "months"), "%Y-%m"), las = 2, cex.axis = 0.7)
legend("topleft", col = c(1, 2, 3,4),lty = c(1, 1, 1, 1),legend = c("real", "AR(5)", "mon/thu","week"))
2.4 状態空間モデルによる予測
2.4.1 ARmodelを状態を記述する方程式として考えている
過去の練習参加者数についてAR
modelをarfitで推定し，関数からの出力，すなわち，AICを最小とするAR次数，AR係数，イノベーション分散を用いて状態空間モデルを定義する．その状態空間モデルを用いて長期予測を行う．これはTSSSpackageのtsmooth()関数という長期予測を行う関数によって計算できる．
なお，このコードは『Rによる時系列モデリング入門』の9.5節「時系列の予測」の食品産業従事者数データの長期予測を参考に記述した．
z1 <- arfit(xts_sample,
plot = FALSE)
z1$maice.order## [1] 5よって，AR modelへの当てはめを考えたとき，AICを最小にする次数は5であり，AR(5) modelを状態方程式として状態空間モデルを組み立てる．
m1 <- z1$maice.order
tau2 <- z1$sigma2
arcoef <- z1$arcoef[[m1]]
f <- matrix(0.0e0 , m1 ,m1)
f[1,] <- arcoef
if (m1 != 1){
for (i in 2:m1){
f[i,i-1] <- 1
}
g <- c(1,rep(0.0e0,m1 -1))
h <- c(1,rep(0.0e0,m1 -1))
q <- tau2[m1 + 1]
r <- 0.0e0
x0 <- rep(0.0e0,m1)
v0 <- NULL
}
s1 <- tsmooth(mon_thu.xts,
f,g,h,q,r,x0,v0,
filter.end = length(mon_thu.xts) ,
predict.end = length(mon_thu.xts)+20,
)
2.5 説明変数として複数パラメータを追加して予測する
2.5.1 加えた説明変数について
練習参加人数はテスト期間か，先生指導の練習か，本番までの残り日数，練習場所などに依存する．そのため，それらの変数も考慮して予測したい． 今回加えた変数は，
- lesson: 先生による練習であれば，1．そうでなければ0となるダミー変数
- VT: 三ヶ月に一回ほど各パートごとにVoice Trainingを行っている．Voice Trainingがある練習日は1, そうでない練習日は0となるダミー変数
- daysLeft: 本番までの残り日数．なおこのデータにおける「本番」としては2024年1月18日に行われた「定期演奏会」および2024年6月30日に行われた「ジョイントコンサート」を採用した．
- not_EX: たまに執行に責任をもつ代(執行代と呼称している)以外が練習を担うことがある．その練習日について1，そうでないときを0とするダミー変数
- new_comer: 5月から新入生が正式に団員となったため，練習参加者として記載されるようになった．それを反映するため2024年4月までの練習日(新入生が計上されていない)を0，5月以降の練習日(新入生が計上されている)を1とするダミー変数
- rainfall: その当日の降水量(なお降水量に関する天気予報が見当たらなかったため予測には用いていない)
- tokyo: 練習場所が東京大学であったときに1, そうでなかったときは0となるダミー変数
- ocha: 練習場所がお茶の水女子大学であったときに1, そうでなかったときは0となるダミー変数
- test: テスト期間中(1, 5, 7, 11月の下旬)の練習であったときに1, そうでなかったときは0となるダミー変数
- camp: 合宿期間中の練習であったときに1, そうでなかったときは0となるダミー変数
- GP: リハーサルの練習であったときに1, そうでなかったときは0となるダミー変数
以上である．降水量のデータは気象庁のサイトから取得した．
また，パラメータdaysLeftについては，本番に近づくほど危機感が募り練習に参加するモチベーションが高まるという経験則から対数変換した．この変換によって，本番が近いほどdaysLeftの変化が残り日数に対して鋭敏になるため，先に述べた「経験則」をより反映すると考えられる．
実際のデータとして2023-06-29から2024-07-29までのデータを，予測用のパラメータを2024年8月からの24日分用意した.(この分析を行った時点で24日先までは各種パラメータの値が確定しているからである．)
metadata_all <- read.csv("/Users/suzukikouhei/Library/Mobile Documents/com~apple~CloudDocs/10_授業/3S/3_5_数理手法7/data分析/metadata_all.csv")
#daysLeftの増減について実感にあうように対数変換
metadata_all$daysLeft <- log(-metadata_all$daysLeft + 1)
#訓練用データ
traindata <- metadata_all[1:117,]
#予測用の説明変数として24日分用意した
newdata1 <- metadata_all[118:158,]2.5.2 分析の方針
まず練習参加者数と加えた説明変数とで重回帰分析を実行する．そして，重回帰分析における誤差項について抽出して時系列として分析することを試みる．練習の参加者数はリハや先生の参加によって大きく変動するものであり，その影響を取り除いた時系列データ，すなわち誤差項を時系列解析の手法で予測することによって，重回帰モデルで説明しきれない時間的な傾向性を捉えることができると考えた．
まず重回帰分析を実行してパラメータに対するp値を出力する．
# 重回帰分析
regression_model <- lm(participants ~ lesson + VT + daysLeft + not_Ex + newcomer + tokyo + ocha + test + camp + GP, data = traindata)
summary(regression_model)##
## Call:
## lm(formula = participants ~ lesson + VT + daysLeft + not_Ex +
## newcomer + tokyo + ocha + test + camp + GP, data = traindata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.0717 -3.1334 -0.2467 2.7471 13.1837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.5532 3.0184 12.773 < 2e-16 ***
## lesson -1.9757 1.3098 -1.508 0.134417
## VT 1.7427 1.4040 1.241 0.217261
## daysLeft -3.2679 0.6655 -4.911 3.31e-06 ***
## not_Ex 4.3618 2.3143 1.885 0.062201 .
## newcomer 8.0876 1.2900 6.270 7.97e-09 ***
## tokyo 1.0602 1.6853 0.629 0.530640
## ocha 2.7888 1.0851 2.570 0.011556 *
## test -6.4271 1.7452 -3.683 0.000365 ***
## camp 19.4845 2.7886 6.987 2.58e-10 ***
## GP 20.7133 3.9316 5.268 7.28e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.208 on 106 degrees of freedom
## Multiple R-squared: 0.6732, Adjusted R-squared: 0.6424
## F-statistic: 21.83 on 10 and 106 DF, p-value: < 2.2e-16パラメータのうち，p-valueが小さいものを除いて，すべてのパラメータのp-valueが0.05を下回るようにする．
# 重回帰分析
regression_model <- lm(participants ~ daysLeft + not_Ex + newcomer + ocha + test + camp + GP, data = traindata)
summary(regression_model)##
## Call:
## lm(formula = participants ~ daysLeft + not_Ex + newcomer + ocha +
## test + camp + GP, data = traindata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1226 -3.1825 -0.6266 3.2657 14.6614
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 38.0726 3.0024 12.681 < 2e-16 ***
## daysLeft -3.1723 0.6631 -4.784 5.42e-06 ***
## not_Ex 4.6805 2.2763 2.056 0.042155 *
## newcomer 8.4321 1.2350 6.828 5.10e-10 ***
## ocha 2.8734 1.0339 2.779 0.006420 **
## test -6.3274 1.7204 -3.678 0.000367 ***
## camp 18.5076 2.7804 6.656 1.17e-09 ***
## GP 19.8049 3.9552 5.007 2.14e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.268 on 109 degrees of freedom
## Multiple R-squared: 0.6561, Adjusted R-squared: 0.634
## F-statistic: 29.71 on 7 and 109 DF, p-value: < 2.2e-16残差をプロットすると，下記のようになり弱定常性は言えると考えられる．
# 残差の抽出
residuals <- residuals(regression_model)
residuals.ts <- as.ts(residuals)
plot(residuals.ts)
出力された残差について弱定常であると目視で確認できる．AR modelによる当てはめを試みると，次数が3のときAICが谷となる．ゆえにAR(3) modelで残差を当てはめる．
arfit(residuals.ts)
AR(3)モデルのフィッティングにはforecastpackageのarima()関数を用いている．
# AR(3)モデルのフィッティング
ar_model <- arima(residuals.ts, order = c(3, 0, 0))
# 残差の予測
n_ahead <- 41
ar_forecast <- predict(ar_model, n.ahead = n_ahead)
forecast_values <- ar_forecast$pred
forecast_se <- ar_forecast$se
regression_forecast <- predict(regression_model, newdata = newdata1)
# 総合予測値の計算
total_forecast <- regression_forecast + forecast_values
total_forecast_lwr <- total_forecast - forecast_se
total_forecast_upr <- total_forecast + forecast_se
# 予測結果をデータフレームに変換
forecast_df <- data.frame(
Date = as.Date(newdata1$date),
Fit = as.numeric(total_forecast),
Lwr = as.numeric(total_forecast_lwr),
Upr = as.numeric(total_forecast_upr)
)
# 観測データの準備
observed_df <- data.frame(
Date = as.Date(traindata$date),
Participants = traindata$participants
)
# プロットの作成
ggplot() +
geom_line(data = observed_df, aes(x = Date, y = Participants), color = "blue") +
geom_line(data = forecast_df, aes(x = Date, y = Fit), color = "red") +
geom_ribbon(data = forecast_df, aes(x = Date, ymin = Lwr, ymax = Upr), alpha = 0.2, fill = "red") +
scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m-%d") +
labs(title = "Participants Forecast",
x = "Date",
y = "Participants") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1))
2.6 中の仕組みはよくわからないけど一番予測性能高い方法(おまけ)
自分はどういう仕組で動いているのかは検討がつかないが，世の中にはKFAS
packageがあり，これを使うと説明変数が複数あり，かつ時系列分析を行うことができる．(おそらく2.5
の方法でAR
modelではなく状態空間モデルをいい感じに使ったものであろう)
#モデルの構造を決める
model <- SSModel(
traindata$participants ~
traindata$lesson +
traindata$VT +
traindata$daysLeft +
traindata$not_Ex +
traindata$newcomer +
traindata$tokyo +
traindata$ocha +
traindata$test +
traindata$camp +
traindata$GP
+
SSMtrend(degree = 1, Q = NA),
H = NA,
)
# パラメータ推定
fit <- fitSSM(model, inits = rep(1,12))
# 推定結果の表示
fit$model## Call:
## SSModel(formula = traindata$participants ~ traindata$lesson +
## traindata$VT + traindata$daysLeft + traindata$not_Ex + traindata$newcomer +
## traindata$tokyo + traindata$ocha + traindata$test + traindata$camp +
## traindata$GP + SSMtrend(degree = 1, Q = NA), H = NA)
##
## State space model object of class SSModel
##
## Dimensions:
## [1] Number of time points: 117
## [1] Number of time series: 1
## [1] Number of disturbances: 1
## [1] Number of states: 11
## Names of the states:
## [1] traindata$lesson traindata$VT traindata$daysLeft
## [4] traindata$not_Ex traindata$newcomer traindata$tokyo
## [7] traindata$ocha traindata$test traindata$camp
## [10] traindata$GP level
## Distributions of the time series:
## [1] gaussian
##
## Object is a valid object of class SSModel.#fit$optim.out
# フィルタリング
#filtering_result <- KFS(fit$model, filtering = "state")
# 平滑化
smoothing_result <- KFS(fit$model, smoothing = "state")
# 結果の表示
print(smoothing_result)## Smoothed values of states and standard errors at time n = 117:
## Estimate Std. Error
## traindata$lesson -0.723 1.287
## traindata$VT 2.093 1.376
## traindata$daysLeft -3.665 0.803
## traindata$not_Ex 3.748 2.203
## traindata$newcomer 7.645 2.998
## traindata$tokyo 2.896 1.808
## traindata$ocha 4.908 1.239
## traindata$test -6.579 1.896
## traindata$camp 19.907 2.715
## traindata$GP 20.159 3.716
## level 37.404 5.286 future_steps <- 20
# 観測データと将来の説明変数を結合
#all_explanatory_vars <- rbind(explanatory_vars, future_explanatory_vars)
new_model <- SSModel(
rep(NA, 158) ~ metadata_all$lesson + metadata_all$VT + metadata_all$daysLeft + metadata_all$not_Ex + metadata_all$newcomer + metadata_all$tokyo + metadata_all$ocha + metadata_all$test + metadata_all$camp + metadata_all$GP +SSMtrend(1, Q = list(fit$model$Q)),
H = fit$model$H
)
new_model1 <- SSModel(
rep(NA, 41) ~ newdata1$lesson + newdata1$VT + newdata1$daysLeft + newdata1$not_Ex + newdata1$newcomer + newdata1$tokyo + newdata1$ocha + newdata1$test + newdata1$camp +newdata1$GP + SSMtrend(1, Q = list(fit$model$Q)),
H = fit$model$H
)
future_steps <- 24
forecast <- predict(fit$model,
n.ahead = future_steps,
newdata = new_model1,
interval = "confidence")結果を図示すると以下のようになる．
# 予測結果をデータフレームに変換
forecast_df <- data.frame(
Date = as.Date(newdata1$date),
Fit = as.numeric(forecast[, "fit"]),
Lwr = as.numeric(forecast[, "lwr"]),
Upr = as.numeric(forecast[, "upr"])
)
# 観測データの準備
observed_df <- data.frame(
Date = as.Date(traindata$date),
Participants = traindata$participants
)
# 観測データと予測データを結合
plot_data <- rbind(
data.frame(Date = as.Date(traindata$date), Participants = traindata$participants, Type = "Observed"),
data.frame(Date = as.Date(newdata1$date), Participants = forecast_df$Fit, Type = "Forecast")
)
# 予測区間のデータフレーム作成
forecast_interval <- data.frame(
Date = as.Date(newdata1$date),
Lwr = forecast_df$Lwr,
Upr = forecast_df$Upr
)
# プロットの作成
ggplot() +
geom_line(data = observed_df, aes(x = Date, y = Participants), color = "blue") +
geom_line(data = forecast_df, aes(x = Date, y = Fit), color = "red") +
geom_ribbon(data = forecast_interval, aes(x = Date, ymin = Lwr, ymax = Upr), alpha = 0.2) +
scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m-%d") +
labs(title = "Participants Forecast",
x = "Date",
y = "Participants") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 45, hjust = 1))