R/compute_statistics.R

In sstn: Self-Similarity Test for Normality

#' @keywords internal

compute_T_obs <- function(x, mu, sigma, M, t_max, H,
                          weight.type = c("normal", "laplace"),
                          beta) {

  weight.type <- match.arg(weight.type)

  if (weight.type == "normal") {
    w <- function(t) exp(-beta * t^2)
  } else {
    w <- function(t) exp(-beta * abs(t))
  }

  n <- length(x)
  xb <- mean(x)
  Sn <- sqrt(mean((x - xb)^2))

  t <- seq(-t_max, t_max, length.out = H)
  dt <- t[2] - t[1]

  phi.hat.Z <- sapply(t, function(tt)
    exp(-1i * tt * xb / Sn) * mean(exp(1i * (tt / Sn) * x))
  )

  phi0.t <- phi0_vec(t)
  U_n <- sqrt(n) * (phi.hat.Z - phi0.t)

  Qhat_lin <- numeric(M)

  for (m in 1:M) {
    u1 <- t / sqrt(m + 1)
    u2 <- t / sqrt(m)

    U.u1 <- interp_complex(t, U_n, u1)
    U.u2 <- interp_complex(t, U_n, u2)

    phi.u1 <- phi0_vec(u1)
    phi.u2 <- phi0_vec(u2)

    Hm.n <- (m + 1) * (U.u1 / phi.u1) - m * (U.u2 / phi.u2)

    integrand <- w(t) * (phi0.t^2) * Mod(Hm.n)^2
    Qhat_lin[m] <- sum(integrand) * dt
  }

  T_obs <- max(abs((Qhat_lin - mu) / sigma))
  T_obs
}

compute_Qhat_lin_only <- function(x, M, t_max, H, w) {
  n  <- length(x)
  xb <- mean(x)
  Sn <- sqrt(mean((x - xb)^2))

  t  <- seq(-t_max, t_max, length.out = H)
  dt <- t[2] - t[1]

  phi.hat.Z <- sapply(t, function(tt)
    exp(-1i * tt * xb / Sn) * mean(exp(1i * (tt / Sn) * x))
  )

  phi0.t <- phi0_vec(t)
  U_n    <- sqrt(n) * (phi.hat.Z - phi0.t)

  Qhat_lin <- numeric(M)

  for (m in 1:M) {
    u1 <- t / sqrt(m + 1)
    u2 <- t / sqrt(m)

    U.u1 <- interp_complex(t, U_n, u1)
    U.u2 <- interp_complex(t, U_n, u2)

    phi.u1 <- phi0_vec(u1)
    phi.u2 <- phi0_vec(u2)

    Hm.n <- (m + 1) * (U.u1 / phi.u1) - m * (U.u2 / phi.u2)

    integrand   <- w(t) * (phi0.t^2) * Mod(Hm.n)^2
    Qhat_lin[m] <- sum(integrand) * dt
  }

  Qhat_lin
}

compute_exact_T <- function(n, B, M, t_max, H = 100,
                             weight.type = c("normal", "laplace"),
                             beta) {

  weight.type <- match.arg(weight.type)

  if (weight.type == "normal") {
    w <- function(t) exp(-beta * t^2)
  } else {
    w <- function(t) exp(-beta * abs(t))
  }

  Qmat <- matrix(NA_real_, nrow = B, ncol = M)

  for (b in 1:B) {

    x_b <- rnorm(n)

    Qmat[b, ] <- compute_Qhat_lin_only(x = x_b, M = M, t_max = t_max, H = H, w = w)
  }

  mu.hat    <- colMeans(Qmat)
  sigma.hat <- apply(Qmat, 2, sd)

  T.vals <- apply(Qmat, 1, function(Qb)
    max(abs((Qb - mu.hat) / sigma.hat))
  )

  list(T.vals = T.vals, mu = mu.hat, sigma = sigma.hat)
}