Rの関数loessの出力 one.delta、two.delta、trace.hatおよびenpをスクラッチで算出

build_hat_matrix <- function(x, weights = NULL, span = 0.75, degree = 2L) {

  x   <- as.matrix(x)
  N   <- nrow(x)
  if (is.null(weights)) weights <- rep(1.0, N)

  # 近傍点数 nf の計算
  nf <- min(N, as.integer(floor(N * span + 1e-5)))

  L <- matrix(0.0, N, N)

  for (j in seq_len(N)) {
    q <- x[j, , drop = TRUE]

    # 二乗距離の計算
    dist2 <- rowSums(sweep(x, 2, q, "-")^2)

    # 近傍 nf 点の選択（partial sort）
    psi <- order(dist2)[seq_len(nf)]

    # バンド幅 rho = nf 番目二乗距離 × max(1, span)
    rho <- dist2[psi[nf]] * max(1.0, span)
    if (rho == 0) rho <- .Machine$double.eps

    # tricube カーネル重みの計算（ehg127 の 2 ループを忠実に再現）
    u      <- sqrt(dist2[psi] / rho)                    # u_i = dist_i / sqrt(rho)
    w_kern <- sqrt(weights[psi] * pmax(0, 1 - u^3)^3)  # sqrt(rw * tricube)

    # 重み付きデザイン行列 B = W^{1/2} X (nf × p) の構築
    dx    <- x[psi, 1, drop = TRUE] - q[1]
    X_raw <- if (degree == 1L) cbind(1, dx) else cbind(1, dx, dx^2)
    B     <- X_raw * w_kern   # 各行 i に w_kern[i] を乗算

    # WLS による hat matrix 行 j の計算
    # B^T B = X^T W X  (W = diag(w_kern^2))
    BtB <- crossprod(B)

    BtB_inv <- tryCatch(
      solve(BtB),
      error = function(e) {
        # ランク落ち時は Moore-Penrose 疑似逆行列
        s   <- svd(BtB)
        tol <- max(dim(BtB)) * .Machine$double.eps * max(s$d)
        s$v %*% diag(ifelse(s$d > tol, 1 / s$d, 0), length(s$d)) %*% t(s$u)
      }
    )

    # 評価点 j は psi の中で距離 0 の点（必ず含まれる）
    j_idx <- which(psi == j)
    xj    <- matrix(X_raw[j_idx, ], nrow = 1)   # 1 × p

    # L[j, psi] = xj %*% (X^T W X)^{-1} X^T W
    #           = xj %*% BtB_inv %*% t(B) * w_kern
    # t(B) * w_kern は各列 i に w_kern[i] を乗算 → X^T W に等しい
    L[j, psi] <- as.numeric(xj %*% BtB_inv %*% t(B) * w_kern)
  }

  L
}

lowesc_R <- function(L) {
  N <- nrow(L)

  # M = L - I = -(I - L)
  M <- L - diag(N)

  # LL = M M^T = (I-L)(I-L)^T
  LL <- M %*% t(M)

  # trace(L) → trace.hat
  trl <- sum(diag(L))

  # trace(LL) → one.delta = Σ_i ||(I-L)[i,:]||²
  delta1 <- sum(diag(LL))

  # trace(LL^2) → two.delta
  delta2 <- sum(diag(LL %*% LL))

  list(trl = trl, delta1 = delta1, delta2 = delta2)
}

compute_loess_stats_scratch <- function(y, x, weights = NULL,
                                        span = 0.75, degree = 2L) {
  N   <- length(y)
  L   <- build_hat_matrix(x, weights, span, degree)
  res <- lowesc_R(L)

  # enp の計算
  enp <- res$delta1 + 2 * res$trl - N

  list(
    one_delta = res$delta1,
    two_delta = res$delta2,
    trace_hat = res$trl,
    enp       = enp,
    L         = L
  )
}

seed <- 20260611
set.seed(seed)
n <- 250
x <- seq(0, 2 * pi, length.out = n)
y <- sin(x) + rnorm(n, 0, 0.2)
w <- rep(1, n)

res <- compute_loess_stats_scratch(y, x, w, span = 0.75, degree = 2L)
cat("スクラッチ one.delta =", res$one_delta, "\n")
cat("スクラッチ two.delta =", res$two_delta, "\n")
cat("スクラッチ trace.hat =", res$trace_hat, "\n")
cat("スクラッチ enp       =", res$enp,       "\n")
cat("-----------------\n")

fit <- loess(y ~ x, span = 0.75, degree = 2, weights = w,
             control = loess.control(surface = "direct",
                                     statistics = "exact",
                                     trace.hat = "exact"))
cat("公式   one.delta =", fit$one.delta, "\n")
cat("公式   two.delta =", fit$two.delta, "\n")
cat("公式   trace.hat =", fit$trace.hat, "\n")
cat("公式   enp       =", fit$enp,       "\n")