Rの関数:dlmSmooth {dlm}

Rの関数から dlmSmooth {dlm} を確認します。

本ポストはこちらの続きです。

const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...

`dlmSmooth`関数の引数

dlmSmooth関数は、dlmパッケージにおいて平滑化（Smoothing）を実行するための中心的な関数です。

フィルタリングが各時点tまでのデータを用いてその時点の状態（レベル、トレンドなど）を推定するのに対し、平滑化は全ての観測期間（t=1からt=Tまで）のデータを用いて、過去の各時点tの状態を再推定する処理です。

それにより、フィルタリングよりも多くの情報を利用できるため、ノイズの影響がさらに除去され、より滑らかで確からしい状態の推定値を得ることができます。

以下に、関数の引数を詳細に説明します。

library(dlm)
args(dlmSmooth)

function (y, ...) 
NULL

`y`

分析対象の時系列データを指定します。後述する「パターン1」のように、dlmモデルオブジェクトを直接dlmSmoothに渡す際に使用します。
データ形式: tsオブジェクト、数値ベクトル (vector)、または行列 (matrix) を指定できます。
省略可能な場合: 後述する「パターン2」のように、dlmFilter関数の結果（dlmFilteredオブジェクト）をdlmSmoothに渡す場合、そのオブジェクト内に既にデータyの情報が含まれているため、この引数は不要です。

`...`

dlmSmooth関数の挙動を決定する部分です。ここには、dlmモデルオブジェクト、またはdlmFilter関数の実行結果のいずれかを渡します。

パターン1：`dlm`オブジェクトを渡す場合 (`dlmSmooth(y, mod, ...)` )

この使い方では、引数yにデータ、そして...の部分にdlmオブジェクトをmodという名前（または名前なし）で渡します。

mod: dlmクラスのオブジェクトを指定します。
- このmodオブジェクトは、観測ノイズの分散Vやシステムノイズの共分散行列Wといった、全てのパラメータが具体的な数値として設定済みである必要があります。通常、dlmMLEで推定したパラメータを使って構築したモデルをここに渡します。
- この形式で呼び出された場合、dlmSmoothはまず内部でdlmFilter(y, mod)を実行してフィルタリングを行い、その結果を用いて平滑化計算を行います。ですので、一行でフィルタリングから平滑化までを実行していることになります。

パターン2：`dlmFiltered`オブジェクトを渡す場合 (`dlmSmooth(object, ...)` )

この使い方では、dlmFilter関数の実行結果であるdlmFilteredオブジェクトを...の部分に渡します。

object: dlmFilter関数の返り値（dlmFilteredクラスのオブジェクト）を指定します。
- この方法は、既にフィルタリングを実行済みの場合、フィルタリングの結果と平滑化の結果の両方を個別に評価したい場合に有用です。
- この場合、dlmFilteredオブジェクト内にデータyやモデルの情報が含まれているため、引数yを明示的に指定する必要はありません。

戻り値について

dlmSmooth関数は、dlmSmoothedクラスのリストを返します。このリストには多くの情報が含まれていますが、特に重要なのは以下の2つです。

s: 各時点における平滑化された状態ベクトルの時系列（行列）です。
- 行列の各行が時点、各列が状態（例：1列目がレベル、2列目がトレンド）に対応します。
- これが分析の主目的となる結果であり、時系列の成分分解などに使用されます。
S: 各時点における、平滑化された状態ベクトルの共分散行列のリストです。
- これは、平滑化された状態推定値の不確実性（分散）を表します。状態の信頼区間を計算する際に使用できます。

コードによる具体例

サンプルデータとモデルを用いて、上記2つのパターンの使い方を示します。

# サンプルデータとモデルの準備
seed <- 20250930
set.seed(seed)
y_sample <- ts(cumsum(rnorm(50, sd = 2)) + rnorm(50, sd = 5), start = 2020)

# モデル構築関数 (状態が1つのモデル)
build_func <- function(par) {
  dlmModPoly(order = 1, dV = exp(par[1]), dW = exp(par[2]))
}

# パラメータの最尤推定
mle_fit <- dlmMLE(y_sample, parm = c(0, 0), build = build_func)

# 推定されたパラメータで最終モデルを構築
final_model <- build_func(mle_fit$par)

# --- パターン1：データ(y)とモデル(mod)を直接渡す ---
smoothed_result1 <- dlmSmooth(y = y_sample, mod = final_model)

# as.matrix()でオブジェクトを行列に変換してから要素を抽出します。
# また、smoothed_result1$s には t=0 の初期状態も含まれるため、
# [-1, ] でその最初の行を除外しています。
s_matrix <- as.matrix(smoothed_result1$s)
smooth_level1 <- ts(s_matrix[-1, 1], start = start(y_sample), frequency = frequency(y_sample))

cat("--- パターン1で得られた平滑化系列の一部を確認 ---\n")
print(head(smooth_level1))

# --- パターン2：dlmFilterの結果を渡す ---
filtered_result <- dlmFilter(y = y_sample, mod = final_model)
smoothed_result2 <- dlmSmooth(filtered_result)

s_matrix2 <- as.matrix(smoothed_result2$s)
smooth_level2 <- ts(s_matrix2[-1, 1], start = start(y_sample), frequency = frequency(y_sample))

cat("\n--- パターン2で得られた平滑化系列の一部を確認 ---\n")
print(head(smooth_level2))


# --- 結果の可視化 ---
plot(y_sample,
  type = "p", pch = 16, col = "grey",
  main = "観測データと平滑化系列",
  xlab = "年", ylab = "値"
)
lines(smooth_level1, col = "firebrick", lwd = 2, lty = 1)
lines(smooth_level2, col = "blue", lwd = 2, lty = "dashed")
legend("topleft",
  legend = c("観測値", "平滑化系列 (パターン1)", "平滑化系列 (パターン2)"),
  col = c("grey", "firebrick", "blue"),
  lty = c(NA, 1, 2),
  pch = c(16, NA, NA),
  lwd = 2,
  bty = "n"
)
grid()

cat("\n--- 両パターンの結果は一致しているか ---\n")
print(all.equal(target = smooth_level1, current = smooth_level2))

--- パターン1で得られた平滑化系列の一部を確認 ---
Time Series:
Start = 2020 
End = 2025 
Frequency = 1 
[1] 2.470889 3.133549 2.464554 3.767059 5.041644 8.113619

--- パターン2で得られた平滑化系列の一部を確認 ---
Time Series:
Start = 2020 
End = 2025 
Frequency = 1 
[1] 2.470889 3.133549 2.464554 3.767059 5.041644 8.113619

--- 両パターンの結果は一致しているか ---
[1] TRUE