対数尤度 を利用して Box-Cox変換 における 最適ラムダ(λ) の算出を R で試みます。
始めにサンプルデータを作成します。なお、サンプルデータは意図的に「分散が非定常」かつ「季節性(周期性)の無い」時系列データとします。
library(ggplot2)
seed <- 20250517
set.seed(seed)
n <- 200 # データ点数
t <- 1:n # 時間インデックス
# 線形トレンドを作成 (時間とともに増加)
trend <- 50 + 0.5 * t
# トレンドの値に比例する標準偏差を持つ正規ノイズを生成
noise_sd_ratio <- 0.15 # トレンドに対するノイズ標準偏差の比率
# 標準偏差が trend * noise_sd_ratio となるノイズ
noise <- rnorm(n, mean = 0, sd = trend * noise_sd_ratio)
# 時系列データを作成 (トレンド + ノイズ)
x <- trend + noise
# Box-Cox変換の前提である「データがすべて正の値である」が満たされない場合は
# x <- pmax(x, 1e-6) によって非常に小さい正の値に置き換える
if (any(x <= 0)) {
x <- pmax(x, 1e-6)
}
# tsオブジェクト (時系列オブジェクト) に変換
x_ts <- ts(x, frequency = 1) # 今回のサンプルデータには周期性がないため frequency = 1
# 生成したデータの確認
p_original <- forecast::autoplot(x_ts) +
ggtitle("Original Time Series (Non-constant Variance)") +
ylab("Value") +
xlab("Time") +
theme_minimal()
print(p_original)
最適なラムダを選択する関数 boxcox_loglik_lambda_calculator を作成します。
なお、同関数は forecast::BoxCox.lambda を引用参照しています。
関数の処理の流れは以下の通りです。
- データ前処理
- 計画行列の作成
- ラムダ候補の評価
- 最適なラムダの選択
与えられたデータ x に対して、指定された範囲 lower から upper まで0.05刻みでラムダ候補を評価し、対数尤度に比例する値を最大化するラムダを探索します。
boxcox_loglik_lambda_calculator <- function(x, lower = -1, upper = 2) {
# 元のxの属性を保持
x_orig_is_ts <- inherits(x, "ts")
original_freq <- if (x_orig_is_ts) stats::frequency(x) else 1
# --- 1. データ前処理 ---
x_clean <- stats::na.omit(c(x)) # ベクトルとして扱いNA除去
n <- length(x_clean)
logx_clean <- log(x_clean)
xdot <- exp(mean(logx_clean)) # x_clean の幾何平均
# --- 2. 計画行列の作成 ---
# 計画行列 X_matrix は元の x の特性（時系列かどうか、周期）に依存し、
# 行数は n (x_clean の長さ) となります。
X_matrix <- NULL
if (x_orig_is_ts) {
if (original_freq > 1) {
# トレンドと季節性を持つ時系列データ
trend_vec <- 1:n
# x_clean の長さに合わせて季節性ファクターを作成
temp_ts_for_season <- stats::ts(seq_len(n), frequency = original_freq)
season_fac <- factor(stats::cycle(temp_ts_for_season))
# 季節性ファクターのレベルが複数あり、かつデータが十分ある場合のみ季節性項を考慮
if (length(levels(season_fac)) > 1 && n >= original_freq) {
df_for_matrix <- data.frame(trend_vec = trend_vec, season_fac = season_fac)
X_matrix <- stats::model.matrix(~ trend_vec + season_fac, data = df_for_matrix)
} else {
# 季節性を考慮できない場合はトレンドのみ
df_for_matrix <- data.frame(trend_vec = trend_vec)
X_matrix <- stats::model.matrix(~trend_vec, data = df_for_matrix)
}
} else {
# トレンドのみの時系列データ (周期 <= 1)
trend_vec <- 1:n
df_for_matrix <- data.frame(trend_vec = trend_vec)
X_matrix <- stats::model.matrix(~trend_vec, data = df_for_matrix)
}
} else {
# 非時系列データ: 切片のみのモデル
df_for_matrix <- data.frame(dummy_for_intercept = rep(1, n))
X_matrix <- stats::model.matrix(~1, data = df_for_matrix)
}
# --- 3. ラムダ候補の評価 ---
lambdas <- seq(lower, upper, by = 0.05)
loglik_values <- numeric(length(lambdas))
for (i in seq_along(lambdas)) {
la <- lambdas[i]
xt_transformed <- numeric(n)
# Box-Cox変換の適用
if (abs(la) > 0.02) { # 標準的なBox-Cox変換
xt_transformed <- (x_clean^la - 1) / la
} else { # laが0に近い場合 (対数変換)、テイラー展開で近似
xt_transformed <- logx_clean * (1 + (la * logx_clean) / 2 *
(1 + (la * logx_clean) / 3 *
(1 + (la * logx_clean) / 4)))
}
# 変換後のデータをスケーリング
scale_factor <- xdot^(la - 1)
z_lambda <- xt_transformed / scale_factor
# z_lambda にモデルをフィットさせ、残差平方和 (RSS) を計算
current_fit <- stats::lm.fit(X_matrix, z_lambda)
rss <- sum(current_fit$residuals^2)
loglik_values[i] <- -n / 2 * log(rss)
}
# --- 4. 最適なラムダの選択 ---
best_idx <- which.max(loglik_values)
return(lambdas[best_idx])
}それでは、サンプルデータに対して最適なラムダを求めます。
boxcox_loglik_lambda_calculator(x = x_ts, lower = -1, upper = 2)[1] 0.45関数 BoxCox.lambda {forecast} による結果と比較します。
forecast::BoxCox.lambda(x = x_ts, method = "loglik", lower = -1, upper = 2)[1] 0.45一致しています。
以上です。
