Rの関数:trendfilter {genlasso}

# パッケージの読み込み
library(ggplot2)
library(dplyr)
library(tidyr)

# 乱数シードの固定
seed <- 20260109
set.seed(seed)

# データ数
n <- 100
x <- 1:n

# 真のトレンド関数の作成（区分的に定義）
# 1-30: 平坦 (値: 10)
# 31-60: 線形上昇
# 61-100: 二次関数的な曲線
y_true <- numeric(n)
y_true[1:30] <- 10
y_true[31:60] <- seq(10, 20, length.out = 30)
# 61-100は滑らかにつなげるための二次関数
y_true[61:100] <- 20 - 0.05 * (1:40)^2 + 1.5 * (1:40)

# 観測データ y_obs (真のトレンド + ノイズ)
# ノイズは正規分布に従う
y_obs <- y_true + rnorm(n, mean = 0, sd = 1.5)

# データフレーム化
df_sim <- data.frame(
  x = x,
  y_obs = y_obs,
  y_true = y_true
)

# データの可視化
p_raw <- ggplot(df_sim, aes(x = x, y = y_obs)) +
  geom_point(color = "grey50", alpha = 0.7) +
  geom_line(aes(y = y_true), color = "black", linewidth = 1, linetype = "dashed") +
  labs(
    title = "シミュレーションデータ（観測値と真のトレンド）",
    subtitle = "黒破線が真のトレンド、灰色点が観測値",
    x = "時点 (Index)",
    y = "値 (Value)"
  ) +
  theme_minimal()

print(p_raw)

# 異なる次数(ord)でモデルを適合
# λ（正則化パラメータ）はクロスバリデーションで最適値を選択します

# 1. 0次トレンドフィルタリング (Fused Lasso: 階段状)
fit_ord0 <- trendfilter(y = y_obs, ord = 0)
cv_ord0 <- cv.trendfilter(fit_ord0) # 最適なλを選択
# 最適なλでの予測値を抽出
pred_ord0 <- predict(fit_ord0, lambda = cv_ord0$lambda.min)$fit

# 2. 1次トレンドフィルタリング (Linear Trend Filtering: 折れ線)
fit_ord1 <- trendfilter(y = y_obs, ord = 1)
cv_ord1 <- cv.trendfilter(fit_ord1)
pred_ord1 <- predict(fit_ord1, lambda = cv_ord1$lambda.min)$fit

# 3. 2次トレンドフィルタリング (Quadratic Trend Filtering: 曲線)
fit_ord2 <- trendfilter(y = y_obs, ord = 2)
cv_ord2 <- cv.trendfilter(fit_ord2)
pred_ord2 <- predict(fit_ord2, lambda = cv_ord2$lambda.min)$fit

# 結果を結合
df_results <- df_sim %>%
  mutate(
    Ord0_Step = as.numeric(pred_ord0),
    Ord1_Linear = as.numeric(pred_ord1),
    Ord2_Quad = as.numeric(pred_ord2)
  ) %>%
  pivot_longer(
    cols = c("Ord0_Step", "Ord1_Linear", "Ord2_Quad"),
    names_to = "Model",
    values_to = "Fitted"
  )

# モデル名のラベルを見やすく変更
df_results$Model <- factor(df_results$Model,
  levels = c("Ord0_Step", "Ord1_Linear", "Ord2_Quad"),
  labels = c("次数0 (階段状)", "次数1 (折れ線)", "次数2 (曲線)")
)

# 結果の可視化
p_res <- ggplot(df_results, aes(x = x)) +
  # 背景に観測データ
  geom_point(aes(y = y_obs), color = "grey70", alpha = 0.5) +
  # 推定されたトレンド
  geom_line(aes(y = Fitted, color = Model), linewidth = 1) +
  # 真のトレンド（比較用）
  geom_line(aes(y = y_true), color = "black", linetype = "dashed", linewidth = 1) +
  scale_color_brewer(palette = "Set1") +
  facet_wrap(~Model, ncol = 1) +
  labs(
    title = "各次数(ord)によるトレンド推定結果の比較",
    x = "時点 (Index)",
    y = "推定値 (Fitted Value)",
    color = "モデル"
  ) +
  theme_minimal()

print(p_res)