Rの関数:supsmu {stats}

Rの関数から supsmu {stats} を確認します。

関数 supsmu とは

supsmu は、Friedman (1984) によって考案された「SuperSmoother（スーパースムーザー）」と呼ばれる散布図平滑化（ノンパラメトリック回帰）手法を実行する関数です。

本手法の特徴は、クロスバリデーション（交差検証）を利用して、データの局所的な曲率に合わせて平滑化のウィンドウ幅（スパン）を自動的かつ動的に変化させる点にあります。

一般的な平滑化手法（固定スパン）では、データ全体に単一のウィンドウ幅を適用します。

しかし、波長が長い緩やかな変化の領域と、波長が短い激しい変化の領域が混在するデータに対しては、「スパンを広げると激しい変化が潰れてしまう」「スパンを狭めると緩やかな領域のノイズを拾ってしまう」というジレンマに陥ります。

supsmu は内部で異なる3つの固定スパンで平滑化を行い、各データ点の周辺において最も予測誤差が小さくなるスパンの計算結果をシームレスに結合します。

それゆえ、変化の乏しい部分ではノイズ除去が優勢し、変化の激しい部分では元の波形に追随するフィッティングを実現します。

関数 supsmu の引数

args(supsmu)
function (x, y, wt = rep(1, n), span = "cv", periodic = FALSE, 
    bass = 0, trace = FALSE) 
NULL
  1. x

    • 説明変数を表す数値ベクトル。
    • 観測データの横軸となる値を指定します。欠損値や無限大の数値は内部で除外されます。
  2. y

    • 応答変数を表す数値ベクトル。
    • 観測データの縦軸となる値を指定します。x と長さが一致している必要があります。
  3. wt

    • 各観測値に対する重みを指定する数値ベクトル。
    • デフォルト: すべての観測値に対して等しい重み rep(1, n)。
    • 測定精度が異なるデータ等において、特定の観測値の影響力を強めたり弱めたりする際に使用します。
  4. span

    • 平滑化のウィンドウ幅（データ全体の範囲に対する割合）を指定するパラメータ。
    • "cv" (デフォルト): クロスバリデーションにより、局所的に最適なスパンを自動選択します。
    • 0 から 1 までの数値: データ全体に対して固定のスパンを適用します。
  5. periodic

    • データが周期的な境界条件を持つかどうかの論理値。
    • デフォルト: FALSE。
    • TRUE を指定する場合、x の値は 0 から 1 の範囲に収まっている必要があります。両端のデータが滑らかに繋がるように平滑化が行われます。
  6. bass

    • 平滑化の度合い（低音強調）を制御するパラメータ。
    • デフォルト: 0 （最大 10 まで）。
    • 数値を大きくするほど、交差検証においてより大きなスパンが選ばれやすくなります。高周波の細かい変動を抑え、より滑らかな曲線を抽出したい場合に数値を引き上げます。
  7. trace

    • 内部の計算プロセスを出力するかどうかの論理値。
    • デフォルト: FALSE。

シミュレーションコード

以下に、supsmu 関数の適応的な平滑化能力を確認するためのシミュレーションコードを示します。

本シミュレーションでは、経時とともに周波数が高くなる波形にノイズを付加したデータを生成します。

具体的には、前半部分は低周波で変動し、後半部分は高周波で変動する波形データに、全体的に高めのノイズを付加したサンプルデータとします。

当該データに対し、固定スパンの手法と、span = "cv" による SuperSmoother をそれぞれ適用し、フィッティングの結果を比較します。

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

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


# 1. サンプルデータの生成
# 局所的に曲率が異なるデータとして、チャープ信号を採用します。
n_obs <- 2000
x_vals <- seq(0, 1, length.out = n_obs)

# 前半から後半にかけて周波数が高くなっていくサインカーブ
true_signal <- sin(2 * pi * (1 + 8 * x_vals) * x_vals)

# 観測データ (真の信号 + 正規ノイズ)
y_obs <- true_signal + rnorm(n_obs, mean = 0, sd = 0.5)

df_sim <- data.frame(
  Time = x_vals,
  Observed = y_obs,
  TrueSignal = true_signal
)

# 2. supsmu 関数による平滑化の実行と比較

# パターンA: 固定スパン (狭い)
sm_short <- supsmu(x_vals, y_obs, span = 0.01)

# パターンB: 固定スパン (広い)
sm_long <- supsmu(x_vals, y_obs, span = 0.4)

# パターンC: SuperSmoother (デフォルト: 自動適応)
sm_cv <- supsmu(x_vals, y_obs, span = "cv")

# パターンD: SuperSmoother (自動適応 + 高周波ノイズ抑制)
sm_bass <- supsmu(x_vals, y_obs, span = "cv", bass = 5)

# 結果をデータフレームに格納
df_sim$Short_Span <- sm_short$y
df_sim$Long_Span <- sm_long$y
df_sim$CV_Auto <- sm_cv$y
df_sim$CV_Bass <- sm_bass$y

# 3. データの整形と可視化
# ggplot2 で描画するために縦持ち (Long format) に変換
df_plot <- pivot_longer(
  df_sim,
  cols = c("Short_Span", "Long_Span", "CV_Auto", "CV_Bass"),
  names_to = "Method",
  values_to = "Estimate"
)

# パネルの並び順と日本語ラベルの設定
df_plot$Method <- factor(
  df_plot$Method,
  levels = c("Short_Span", "Long_Span", "CV_Auto", "CV_Bass"),
  labels = c(
    "1. 固定スパン (0.01)",
    "2. 固定スパン (0.40)",
    "3. 自動適応 (span='cv')",
    "4. 自動適応＋高周波抑制 (bass=5)"
  )
)

# プロットの作成
p1 <- ggplot(df_plot, aes(x = Time)) +
  # 観測データ (灰色の点)
  geom_point(aes(y = Observed), color = "gray70", size = 1, alpha = 0.5) +
  # 真の信号 (黒の破線)
  geom_line(aes(y = TrueSignal), color = "black", linetype = "dashed", linewidth = 0.5) +
  # 各手法による平滑化結果 (赤の実線)
  geom_line(aes(y = Estimate), color = "red", linewidth = 1) +
  facet_wrap(~Method, ncol = 2) +
  labs(
    title = "supsmu 関数における適応的平滑化",
    subtitle = "赤線が推定結果、黒破線が真の波形",
    x = "時間 (Time)",
    y = "観測値および推定値"
  ) +
  theme_minimal() +
  theme()

print(p1)
Figure 1

以上です。