Rの関数:smooth {stats}

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

関数 smooth とは

smooth は、「移動中央値平滑化」を実行するための関数です。

時系列データや連続的な観測データのノイズを除去する際、一般的な移動平均を用いると、外れ値（スパイクノイズ）の影響を強く受けて波形が歪んでしまう場合がありますが、中央値に基づく平滑化は平均値と比較して外れ値に対して頑健（ロバスト）な結果が得られます。

本関数は、単純な幅3の移動中央値にとどまらず、反復処理や分割処理を組み合わせた平滑化アルゴリズムをC言語のコードで提供しています。

関数 smooth の引数

args(smooth)
function (x, kind = c("3RS3R", "3RSS", "3RSR", "3R", "3", "S"), 
    twiceit = FALSE, endrule = c("Tukey", "copy"), do.ends = FALSE) 
NULL
  1. x

    • 平滑化の対象となる数値ベクトル、または時系列オブジェクト。
    • 欠損値（NA）を含めることはできず、実行前に除外または補完する必要があります。
  2. kind

    • 適用する平滑化のアルゴリズムを指定する文字列です。
    • "3": 幅3の移動中央値を1回だけ適用します。
    • "3R": 幅3の移動中央値を、データが変化しなくなる（収束する）まで繰り返し適用します（RはRepeatedを意味します）。
    • "S": 平坦な頂点や谷底を分割（Splitting）して滑らかにします。
    • "3RSS", "3RSR", "3RS3R" (デフォルト): 上記の手法を組み合わせた複合的なアルゴリズムです。デフォルトの "3RS3R" は、「幅3の反復」\(\rightarrow\)「分割処理」\(\rightarrow\)「幅3の反復」を順次行う平滑化手法です。
  3. twiceit

    • 平滑化の残差に対して「Twicing」と呼ばれる手続きを行うかどうかの論理値です。
    • デフォルト: FALSE。
    • TRUE を指定すると、元のデータから1回目の平滑化結果を引いた「残差」に対して再度同じ平滑化アルゴリズムを適用し、当該結果を1回目の結果に加算します。
  4. endrule

    • データの両端（端点）における処理方法を指定する文字列です。
    • "Tukey" (デフォルト): Tukeyの端点処理ルール（端点と内側の値の中央値をとる手法など）を適用します。
    • "copy": 端点のデータをそのままコピーします。
  5. do.ends

    • 端点に対しても分割（Splitting）の手続きを適用するかどうかの論理値です。
    • デフォルト: FALSE。

シミュレーションコード

以下に、smooth 関数の機能を確認するためのシミュレーションコードを示します。

本シミュレーションでは、滑らかなサイン波に「突発的で極端なスパイクノイズ」を混入させたサンプルデータを生成します。

当該データに対し、異なる kind 引数を用いた平滑化を適用し、結果を比較します。

サンプルデータは滑らかな波形に対し、20, 50, 80の3時点に極端な外れ値を付与します。

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

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

# 1. 基礎データの生成
n_time <- 100
time_seq <- 1:n_time

# ベースとなる滑らかな波形 (サイン波)
true_wave <- sin(2 * pi * time_seq / 25)

# 通常の小さな測定ノイズを付加
obs_data <- true_wave + rnorm(n_time, mean = 0, sd = 0.2)

# 極端なスパイクノイズ（外れ値）を3箇所に混入
obs_data[20] <- 5.0
obs_data[50] <- -4.0
obs_data[80] <- 6.0

# 2. smooth 関数を用いた平滑化の実行と手法間の比較

# 手法A: 単純な幅3の移動中央値 (1回のみ)
smooth_3 <- smooth(obs_data, kind = "3")

# 手法B: 幅3の移動中央値を収束するまで反復
smooth_3R <- smooth(obs_data, kind = "3R")

# 手法C: 複合アルゴリズム (デフォルト)
smooth_3RS3R <- smooth(obs_data, kind = "3RS3R")

# 手法D: 複合アルゴリズム + Twicing (残差の再適用)
smooth_Twice <- smooth(obs_data, kind = "3RS3R", twiceit = TRUE)

# 3. データの整形と可視化
df_smooth <- data.frame(
  Time = time_seq,
  Observed = obs_data,
  TrueWave = true_wave,
  Method_3 = as.numeric(smooth_3),
  Method_3R = as.numeric(smooth_3R),
  Method_3RS3R = as.numeric(smooth_3RS3R),
  Method_Twice = as.numeric(smooth_Twice)
)

# ggplot 用にデータを縦持ちに変換
df_plot <- pivot_longer(
  df_smooth,
  cols = starts_with("Method"),
  names_to = "Method",
  values_to = "Smoothed"
)

# ファセット表示用のラベルを設定
df_plot$Method <- factor(
  df_plot$Method,
  levels = c("Method_3", "Method_3R", "Method_3RS3R", "Method_Twice"),
  labels = c(
    "1. kind = '3' (1回のみ)",
    "2. kind = '3R' (収束まで反復)",
    "3. kind = '3RS3R' (デフォルト・複合)",
    "4. kind = '3RS3R' + twiceit=TRUE"
  )
)

# プロットの作成
p1 <- ggplot(df_plot, aes(x = Time)) +
  # 観測データ
  geom_line(aes(y = Observed), color = "gray80", linewidth = 0.5) +
  geom_point(aes(y = Observed), color = "gray50", size = 2) +

  # 真の波形
  geom_line(aes(y = TrueWave), color = "blue", linewidth = 2.5, alpha = 0.5) +

  # 平滑化結果
  geom_line(aes(y = Smoothed), color = "red", linewidth = 1.0) +
  facet_wrap(~Method, ncol = 2) +
  labs(
    title = "Tukeyの平滑化アルゴリズム (smooth関数) の比較",
    x = "時間",
    y = "値"
  ) +
  theme_minimal() +
  theme(strip.text = element_text(size = 11, face = "bold"))

print(p1)
Figure 1

Figure 1 の青の実践は真の波形、赤の実線は平滑化結果、灰色の点は観測データです。

kind='3' と kind='3R' および kind='3RS3R' と kind='3RS3R' + twiceit=TRUE' とでは、20点目のスパイクノイズに対する反応の違いが目視で確認できますが、今回のサンプルデータでは4つのアルゴリズムいずれも真の波形に概ね追随しています。

純Rコードでの平滑化

関数smoothはラッパーであり、計算本体は以下のリンク先のC言語関数が担っています。

smooth.c の中では主に以下の処理が実装されています。

  • Tukey のスムージング手法（running median）
    • 3-point median（“3”）
    • repeated smoothing（“3R” 系）
    • split smoothing（“S”）
  • 端点処理（Tukey / copy）
  • 反復処理（収束まで）

そこで、上記のsmooth.cに頼らず、純Rコードによる移動中央値平滑化をシミュレーションします。

なお、平滑化のアルゴリズムは3Rとし、端点処理はTukeyとします。

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

# ------------------------------------------------------------------
# 1. スクラッチコード
# ------------------------------------------------------------------

# 3つの値の中央値を返す関数
med3 <- function(a, b, c) {
  return(sort(c(a, b, c))[2])
}

# 「3」の処理: 幅3の移動中央値とTukeyの端点処理
sm_3_R <- function(x, n, iend) {
  y <- numeric(n)
  if (n < 3) {
    return(list(y = x, iter = 1))
  }

  # 中央部分の処理
  for (i in 2:(n - 1)) {
    y[i] <- med3(x[i - 1], x[i], x[i + 1])
  }

  # 端点処理 (iend = 2 は "Tukey" ルール)
  if (iend == 1) {
    y[1] <- x[1]
    y[n] <- x[n]
  } else if (iend == 2) {
    # 外挿予測（内側の傾きから端の値を推測し、それと実データの中央値をとる手法）
    y[1] <- med3(x[1], y[2], 3.0 * y[2] - 2.0 * y[3])
    y[n] <- med3(x[n], y[n - 1], 3.0 * y[n - 1] - 2.0 * y[n - 2])
  } else {
    y[1] <- x[1]
    y[n] <- x[n]
  }

  return(list(y = y, iter = 1))
}

# 「3R」の処理: 収束するまで「3」を反復
sm_3R_R <- function(x, n, iend) {
  y_current <- x
  iter <- 0
  max_iter <- 100

  repeat {
    iter <- iter + 1
    res <- sm_3_R(y_current, n, iend)
    y_next <- res$y

    # 一致（収束）したら終了
    if (identical(y_current, y_next) || iter >= max_iter) {
      break
    }
    y_current <- y_next
  }

  return(list(y = y_current, iter = iter))
}


# ------------------------------------------------------------------
# 2. 検証用データの生成
# ------------------------------------------------------------------
n_time <- 80
time_seq <- 1:n_time

# サイン波にノイズとスパイク（外れ値）を混入
true_wave <- sin(2 * pi * time_seq / 20)
obs_data <- true_wave + rnorm(n_time, mean = 0, sd = 0.3)
obs_data[15] <- 4.0 # 強いスパイク
obs_data[65] <- -3.5 # 強いスパイク

# ------------------------------------------------------------------
# 3. 両手法による平滑化の実行
# ------------------------------------------------------------------

# 手法A: 関数 smooth
# kind = "3R", 端点処理はデフォルトの "Tukey" (iend = 2 に相当)
res_smooth <- smooth(obs_data, kind = "3R", endrule = "Tukey")

# 手法B: スクラッチ実装
res_scratch <- sm_3R_R(obs_data, n = n_time, iend = 2)

# ------------------------------------------------------------------
# 4. 結果の数値的比較
# ------------------------------------------------------------------

# ベクトルとして値が一致するかを確認
val_builtin <- as.numeric(res_smooth)
val_scratch <- res_scratch$y

is_exact_match <- identical(val_builtin, val_scratch)

cat(sprintf("結果は一致しているか: %s\n", is_exact_match))
cat(sprintf("収束までの反復回数 (関数smooth): %d 回\n", attr(res_smooth, "iter")))
cat(sprintf("収束までの反復回数 (スクラッチ実装): %d 回\n\n", res_scratch$iter))
結果は一致しているか: TRUE
収束までの反復回数 (関数smooth): 3 回
収束までの反復回数 (スクラッチ実装): 4 回
# ------------------------------------------------------------------
# 5. 視覚的確認
# ------------------------------------------------------------------
df_plot <- data.frame(
  Time = time_seq,
  Observed = obs_data,
  Builtin = val_builtin,
  Scratch = val_scratch
)

p1 <- ggplot(df_plot, aes(x = Time)) +
  # 観測データ
  geom_point(aes(y = Observed), color = "gray70", size = 2.5, alpha = 0.8) +
  geom_line(aes(y = Observed), color = "gray80", linewidth = 0.5) +

  # R標準関数の結果 (太い青線)
  geom_line(aes(y = Builtin, color = "関数smooth (C言語)"), linewidth = 1.5, alpha = 0.5) +

  # スクラッチ実装の結果 (細い赤の破線)
  geom_line(aes(y = Scratch, color = "スクラッチ実装 (R言語)"), linewidth = 1, linetype = "dashed") +
  scale_color_manual(
    name = "推計手法",
    values = c("関数smooth (C言語)" = "blue", "スクラッチ実装 (R言語)" = "red")
  ) +
  labs(
    title = "平滑化アルゴリズム (kind='3R') の一致確認",
    x = "時間",
    y = "値"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p1)
Figure 2

Figure 2 からも関数smoothの結果とスクラッチ実装の結果が一致していることを確認できます。

以上です。