高知アグリ・データサイエンス・ラボでは、農家が大掛かりな設備投資をしなくても、データを活用した農業を実践できるにはどうすればいいのか、を研究しています。
その一環として、無償のプログラミング言語による収穫量予測に取り組んでいます。
本稿では、状態空間モデルによる収穫量予測について、簡単にレポートします。
なお以下の試みは、「手持ちのデータをローカル・レベル・モデルに適用してみた」ものであり、精度の高い予測ではありません。今後は、予測に必要なデータを追加していき、より精度の高い予測を実現したいと考えています。
1. 状態空間モデルとは
まず、状態空間モデルのおおまかな特徴を説明します。
- 状態空間モデルは、「互いに関連のある系列的なデータを確率的にとらえる」
- 状態空間モデルは、「直接知ることができない「真の状態」の変化と観測を分けて考える」
つまり状態空間モデルでは、直接観測されるデータに加え、直接的には観測されない潜在的な確率変数が導入されます。この潜在的な確率変数が状態と呼ばれます。この状態に関して、 状態は前時点のみと関連があるという関係性 (マルコフ性) を仮定します。さらに観測値に関して、ある時点の観測値はその時点の状態によってのみ決まると仮定します。
- 状態空間モデルでは「将来の予測も欠損値の補間も同じ枠組みで対応ができる」
観測値が欠損している場合は、欠損する直前の状態推定値から将来予測を行います。欠損している期間を予測して行うことで、欠損値を補間できます。
2. 状態空間モデルによる予測: R と Stan を用いて
状態空間モデルで実際に予測を行うために、プログラミング言語である R と Stan を使用しました。
3. グロリオーサの採花量予測
状態空間モデルにはさまざまな種類がありますが、ここでは、最も単純な「ローカルレベルモデル」を実装してみましょう。データは、高知県の特産花卉、グロリオーサの採花量 (2020 年 4 月 17 日 〜 2020 年 12 月 2 日) です。この期間の採花量を元に、ローカルレベルモデルを使って、20 日間の採花量予測を行いました。なおデータファイルは省略しますが、もし中身をみたい場合は、高知アグリデータサイエンスラボの Facebook ページから、「メッセージ機能」を通じてコンタクトください。
3.1 実装コードおよび結果
以下は、R と Stan による実装コードです。コードについてより詳しくは、参考文献を参照してください。
3.1.1 ローカルレベルモデル
# ローカルレベルモデル
# 分析の準備 -------------------------------------------------------------------
# パッケージの読み込み
library(rstan)
library(bayesplot)
library(ggfortify)
library(gridExtra)
# 計算の高速化
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# データの読み込みとPOSIXctへの変換 ----------------------------------------------------
# データの読み込み
sales_df <- read.csv("file_data_wide_misato.csv")
# 日付をPOSIXct形式にする
sales_df$date <- as.POSIXct(sales_df$date)
# データの先頭行を表示
head(sales_df, n = 3)
# POSIXctの補足
POSIXct_time <- as.POSIXct("1970-01-01 00:00:05", tz="UTC")
as.numeric(POSIXct_time)
# ローカルレベルモデルの推定 -------------------------------------------------------------
# データの準備
data_list <- list(
y = sales_df$misato,
T = nrow(sales_df)
)
# モデルの推定
local_level_stan <- stan(
file = "5-2-1-local-level.stan",
data = data_list,
seed = 1
)
# 収束の確認
mcmc_rhat(rhat(local_level_stan))
# 結果の表示
print(local_level_stan,
pars = c("s_w", "s_v","lp__"),
probs = c(0.025, 0.5, 0.975))
# 結果の図示 -------------------------------------------------------------------
# 生成された乱数を格納
mcmc_sample <- rstan::extract(local_level_stan)
# Stanにおける状態を表す変数名
state_name <- "mu"
# 1時点目の状態の95%ベイズ信用区間と中央値を得る
quantile(mcmc_sample[[state_name]][, 1],
probs=c(0.025, 0.5, 0.975))
# すべての時点の状態の、95%ベイズ信用区間と中央値
result_df <- data.frame(t(apply(
X = mcmc_sample[[state_name]],# 実行対象となるデータ
MARGIN = 2, # 列を対象としてループ
FUN = quantile, # 実行対象となる関数
probs=c(0.025, 0.5, 0.975) # 上記関数に入れる引数
)))
# 列名の変更
colnames(result_df) <- c("lwr", "fit", "upr")
# 時間軸の追加
result_df$time <- sales_df$date
# 観測値の追加
result_df$obs <- sales_df$misato
# 図示のためのデータの確認
head(result_df, n = 3)
# 図示
ggplot(data = result_df, aes(x = time, y = obs)) +
labs(title="result of local level model") +
ylab("misato") +
geom_point(alpha = 0.6, size = 0.9) +
geom_line(aes(y = fit), size = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3) +
scale_x_datetime(date_labels = "%Y-%m")
# 図示をする関数 -----------------------------------------------------------------
plotSSM <- function(mcmc_sample, time_vec, obs_vec = NULL,
state_name, graph_title, y_label,
date_labels = "%Y-%m"){
# 状態空間モデルを図示する関数
#
# Args:
# mcmc_sample : MCMCサンプル
# time_vec : 時間軸(POSIXct)のベクトル
# obs_vec : (必要なら)観測値のベクトル
# state_name : 図示する状態の変数名
# graph_title : グラフ
# y_label : y軸のラベル
# date_labels : 日付の書式
#
# Returns:
# ggplot2により生成されたグラフ
# すべての時点の状態の、95%区間と中央値
result_df <- data.frame(t(apply(
X = mcmc_sample[[state_name]],
MARGIN = 2, quantile, probs = c(0.025, 0.5, 0.975)
)))
# 列名の変更
colnames(result_df) <- c("lwr", "fit", "upr")
# 時間軸の追加
result_df$time <- time_vec
# 観測値の追加
if(!is.null(obs_vec)){
result_df$obs <- obs_vec
}
# 図示
p <- ggplot(data = result_df, aes(x = time)) +
labs(title = graph_title) +
ylab(y_label) +
geom_line(aes(y = fit), size = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3) +
scale_x_datetime(date_labels = date_labels)
# 観測値をグラフに追加
if(!is.null(obs_vec)){
p <- p + geom_point(alpha = 0.6, size = 0.9,
data = result_df, aes(x = time, y = obs))
}
# グラフを返す
return(p)
}
plotSSM(mcmc_sample = mcmc_sample, time_vec = sales_df$date,
obs_vec = sales_df$misato,
state_name = "mu", graph_title = "result of local level model",
y_label = "misato")
3.1.2 ローカルレベルモデルの結果
点 = 「・」が実際の観測値で、実線が「状態」= 予測値です。また、濃い灰色の部分は 95 %区間です。つまり、点と実線が近ければ、実際の観測値と予測された状態が近い、要するに「予測が当たっている」ということです。点が濃い灰色の範囲内であれば、予測が「惜しい」ということです。
3.1.3 ローカルレベルモデルによる予測
# 状態空間モデルによる予測
# 分析の準備 -------------------------------------------------------------------
# パッケージの読み込み
library(rstan)
library(bayesplot)
# 計算の高速化
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# 状態空間モデルの図示をする関数の読み込み
source("plotSSM.R", encoding="utf-8")
# データの読み込み
harvests_df_all <- read.csv("file_data_wide_misato.csv")
harvests_df_all$date <- as.POSIXct(harvests_df_all$date)
# ローカルレベルモデルによる予測の実行 ----------------------------------------------------------------------
# データの準備
data_list_pred <- list(
T = nrow(harvests_df_all),
y = harvests_df_all$misato,
pred_term = 20
)
# モデルの推定
local_level_pred <- stan(
file = "5-3-1-local-level-pred.stan",
data = data_list_pred,
seed = 1
)
# 参考：収束の確認
mcmc_rhat(rhat(local_level_pred))
# 参考：結果の表示
print(local_level_pred,
pars = c("s_w", "s_v","lp__"),
probs = c(0.025, 0.5, 0.975))
## 図示
# 予測対象期間も含めた日付を用意
date_plot <- seq(
from = as.POSIXct("2020-04-17"),
by = "days",
len = 250)
# 参考
seq(from = as.POSIXct("2020-04-17"),
by = 60*60*24,
len = 250)
# 生成された乱数を格納
mcmc_sample_pred <- rstan::extract(local_level_pred)
# 予測結果の図示
plotSSM(mcmc_sample = mcmc_sample_pred,
time_vec = date_plot,
state_name = "mu_pred",
graph_title = "Prediciton Result",
y_label = "misato",
date_labels = "%Y-%m")
3.1.3 ローカルレベルモデルによる予測結果
20 日間の予測を実行すると、状態値がほぼ変化せず、95 %区間が拡大するだけの結果になってしまいました。
3.2 考察と今後の展望
状態値がほぼ変化せず、「役に立つ」予測結果が得られなかった理由は、以下の 2 点だと考えられます。つまり、(1) データが少ないこと、(2) そもそも「採花量「のみ」」からしか予測していないことです。ローカル・レベル・モデルは、状態空間モデルのなかでも単純なモデルです。採花量と、採花量以外の時系列データを組み合わせ・相関関係を算出することで、より正確な採花量予測ができるかもしれません。
今後は、「「ふつうの農家」が収集しやすいデータ」= 「気温・湿度・日射量 (など)」を組み込んだ多変量時系列解析に取り組みたいと考えています。
