Rの関数:pwilcox {stats}

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

# 1. パラメータ設定
# 小規模なサンプルサイズを設定（計算過程を理解しやすくするため）
m_size <- 5 # グループAのサンプル数
n_size <- 4 # グループBのサンプル数

# 統計量 W のとりうる最大値は m * n
w_max <- m_size * n_size

cat("--- 設定パラメータ ---\n")
cat(sprintf("グループAのサイズ (m): %d\n", m_size))
cat(sprintf("グループBのサイズ (n): %d\n", n_size))
cat(sprintf("統計量 W の範囲: 0 ~ %d\n\n", w_max))

# 2. 理論的な累積分布の計算 (pwilcoxを使用)
# 0 から w_max までの全ての q について確率を計算
q_vals <- 0:w_max
prob_vals <- pwilcox(q = q_vals, m = m_size, n = n_size)

# データフレーム化
df_theoretical <- data.frame(
  q = q_vals,
  Probability = prob_vals,
  Type = "理論値 (pwilcox)"
)

cat("--- pwilcoxによる計算例 ---\n")
cat(sprintf("W <= 5 となる確率: %.4f\n", pwilcox(5, m_size, n_size)))
cat(sprintf("W <= 10 となる確率: %.4f (分布の中央)\n", pwilcox(10, m_size, n_size)))
cat(sprintf("W <= 15 となる確率: %.4f\n\n", pwilcox(15, m_size, n_size)))

# stats::pwilcox のスクラッチ実装
scratch_pwilcox <- function(q, m, n, lower.tail = TRUE, log.p = FALSE) {
  # 1. 引数のバリデーション
  if (length(m) != 1 || length(n) != 1) stop("m と n は単一の整数でなければなりません")
  m <- as.integer(m)
  n <- as.integer(n)

  # エッジケース: サンプルサイズが0の場合
  if (m == 0 || n == 0) {
    # 統計量は常に0になるため、q >= 0 なら確率は1
    p <- ifelse(q >= 0, 1, 0)
    if (!lower.tail) p <- 1 - p
    if (log.p) {
      return(log(p))
    } else {
      return(p)
    }
  }

  # 2. 動的計画法 (DP) による分布の計算
  # 目的: サイズ(m, n)の時に、統計量Wがkになる組み合わせの数を計算する
  # 漸化式: count(k, m, n) = count(k - n, m - 1, n) + count(k, m, n - 1)
  # 意味: 「最大のランクを持つ要素」がグループmに属するか、グループnに属するかで場合分け

  # 初期化: i=0 (m=0) の状態。W=0 となる組み合わせが1通りだけある。
  # current_dists[[j]] には、サイズ(0, j-1) における分布のベクトルを格納
  # リストのインデックスはRの仕様上 1 始まりのため、nに対応するデータは index=n+1 に格納

  current_dists <- rep(list(c(1)), n + 1)

  # i (グループmのサイズ) を 1 から m まで増やしていく
  for (i in 1:m) {
    next_dists <- list()

    # j=0 (グループnのサイズが0) の場合、W=0 の1通りのみ
    next_dists[[1]] <- c(1)

    # j (グループnのサイズ) を 1 から n まで増やしていく
    for (j in 1:n) {
      # 漸化式の適用

      # 1. count(k - j, i - 1, j) に対応
      # 直前のステップ (i-1, j) の分布を取得し、スコアを j だけシフト（加算）する
      # (Rのリストインデックスは j+1 を参照)
      dist_prev_i <- current_dists[[j + 1]]
      shifted_dist <- c(rep(0, j), dist_prev_i) # 先頭に0をj個追加してシフト

      # 2. count(k, i, j - 1) に対応
      # 同じiのステップでの (i, j-1) の分布を取得
      dist_prev_j <- next_dists[[j]] # こちらは既に計算済み

      # 3. 2つの分布（ベクトル）を足し合わせる
      # ベクトルの長さが異なるため、短い方の末尾を0で埋めて長さを合わせる
      len_diff <- length(shifted_dist) - length(dist_prev_j)

      if (len_diff > 0) {
        dist_prev_j <- c(dist_prev_j, rep(0, len_diff))
      } else if (len_diff < 0) {
        shifted_dist <- c(shifted_dist, rep(0, -len_diff))
      }

      # 加算して新しい分布とする
      next_dists[[j + 1]] <- shifted_dist + dist_prev_j
    }

    # 次のイテレーションのために更新
    current_dists <- next_dists
  }

  # 最終的な頻度分布 (サイズ m, n)
  freqs <- current_dists[[n + 1]]

  # 3. 累積分布関数 (CDF) の計算
  # 総組み合わせ数 choose(m + n, m)
  total_comb <- sum(freqs)

  # 確率質量分布 (PMF) -> 累積分布 (CDF)
  probs <- cumsum(freqs) / total_comb

  # 4. 入力 q に対応する確率の取得
  # q はベクトル対応する必要がある
  # 統計量は離散値なので、floor(q) 以下の累積確率を参照する
  q_idx <- floor(q)
  res_p <- numeric(length(q))

  max_w <- length(probs) - 1 # 最大スコア m*n

  for (k in seq_along(q)) {
    val <- q_idx[k]
    if (val < 0) {
      res_p[k] <- 0
    } else if (val >= max_w) {
      res_p[k] <- 1
    } else {
      # freqs[1] が W=0 に対応するため、インデックスは val + 1
      res_p[k] <- probs[val + 1]
    }
  }

  # 5. オプション処理 (lower.tail, log.p)
  if (!lower.tail) res_p <- 1 - res_p
  if (log.p) res_p <- log(res_p)

  return(res_p)
}

# テストパラメータ
m <- 5
n <- 4
q_vals <- c(-1, 0, 5, 10, 15, 20, 25) # 境界値を含むテスト

# stats::pwilcox の結果
res_stats <- pwilcox(q_vals, m, n)

# 自作関数の結果
res_my <- scratch_pwilcox(q_vals, m, n)

# 結果の表示
df_compare <- data.frame(
  q = q_vals,
  stats_pkg = res_stats,
  my_func = res_my,
  Is_Equal = abs(res_stats - res_my) < 1e-10 # 誤差許容
)

print(df_compare)

cat("\n全ての値が一致しているか確認: ", all(df_compare$Is_Equal), "\n")