R/asymptotic_calibration.R

In sstn: Self-Similarity Test for Normality

#' @keywords internal
#' @importFrom MASS mvrnorm

calibrate_mu_sigma <- function(M = 20, t_max = 4, H = 100,
                               weight.type = c("normal", "laplace"),
                               beta = NULL,
                               B = 10000, seed = 18062021) {

  set.seed(seed)

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

  if (is.null(beta)) {
    beta <- if (weight.type == "normal") 1 else 4
  }

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

  t <- seq(-t_max, t_max, length.out = H)
  dt <- t[2] - t[1]
  phi0 <- phi0_vec(t)
  mats <- make_K_C(t)
  K <- mats$K
  C <- mats$C

  cov.Re <- 0.5 * (K + C)
  cov.Im <- 0.5 * (K - C)

  J <- length(t)

  var.Y1 <- 1
  var.Y2 <- 1/2

  cov.Re.Y1 <- rep(0, J)
  cov.Im.Y1 <- t * phi0

  cov.Re.Y2 <- -(t^2/2) * phi0
  cov.Im.Y2 <- rep(0, J)

  p <- 2*J + 2
  SIG <- matrix(0, p, p)

  SIG[1:J, 1:J] <- cov.Re
  SIG[(J+1):(2*J), (J+1):(2*J)] <- cov.Im

  SIG[1:J, 2*J+1] <- cov.Re.Y1
  SIG[J + (1:J), 2*J+1] <- cov.Im.Y1
  SIG[2*J+1, 1:J] <- cov.Re.Y1
  SIG[2*J+1, J + (1:J)] <- cov.Im.Y1

  SIG[1:J, 2*J+2] <- cov.Re.Y2
  SIG[J + (1:J), 2*J+2] <- cov.Im.Y2
  SIG[2*J+2, 1:J] <- cov.Re.Y2
  SIG[2*J+2, J + (1:J)] <- cov.Im.Y2

  SIG[2*J+1, 2*J+1] <- var.Y1
  SIG[2*J+2, 2*J+2] <- var.Y2

  SIG <- (SIG + t(SIG))/2
  SIG <- SIG + diag(1e-10, p)

  Q.collect <- matrix(NA_real_, nrow = B, ncol = M)

  for (b in 1:B) {

    z <- mvrnorm(1, mu = rep(0, p), Sigma = SIG)

    Gre <- z[1:J]
    Gim <- z[(J+1):(2*J)]
    Y1  <- z[2*J + 1]
    Y2  <- z[2*J + 2]
    G   <- Gre + 1i * Gim

    G0 <- G - 1i * t * phi0 * Y1 + (t^2) * phi0 * Y2

    Qm <- numeric(M)

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

      G0.u1 <- interp_complex(t, G0, u1)
      G0.u2 <- interp_complex(t, G0, u2)

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

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

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

    Q.collect[b, ] <- Qm
  }

  mu.hat <- colMeans(Q.collect)
  sd.hat <- apply(Q.collect, 2, sd)

  list(mu = mu.hat, sigma = sd.hat, t = t, dt = dt,
       weight.type = weight.type, beta = beta)
}

simulate_T <- function(mu, sigma, M = 20, t_max = 4, H = 100,
                        weight.type = c("normal", "laplace"),
                        beta = NULL,
                        B = 10000, seed = 18092024) {

  set.seed(seed)

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

  if (is.null(beta)) {
    beta <- if (weight.type == "normal") 1 else 4
  }

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

  t <- seq(-t_max, t_max, length.out = H)
  dt <- t[2] - t[1]
  phi0 <- phi0_vec(t)
  mats <- make_K_C(t)
  K <- mats$K
  C <- mats$C

  cov.Re <- 0.5 * (K + C)
  cov.Im <- 0.5 * (K - C)

  J <- length(t)
  var.Y1 <- 1
  var.Y2 <- 1/2

  cov.Re.Y1 <- rep(0, J)
  cov.Im.Y1 <- t * phi0

  cov.Re.Y2 <- -(t^2/2) * phi0
  cov.Im.Y2 <- rep(0, J)

  p <- 2*J + 2
  SIG <- matrix(0, p, p)

  SIG[1:J, 1:J] <- cov.Re
  SIG[(J+1):(2*J), (J+1):(2*J)] <- cov.Im

  SIG[1:J, 2*J+1] <- cov.Re.Y1
  SIG[J + (1:J), 2*J+1] <- cov.Im.Y1
  SIG[2*J+1, 1:J] <- cov.Re.Y1
  SIG[2*J+1, J + (1:J)] <- cov.Im.Y1

  SIG[1:J, 2*J+2] <- cov.Re.Y2
  SIG[J + (1:J), 2*J+2] <- cov.Im.Y2
  SIG[2*J+2, 1:J] <- cov.Re.Y2
  SIG[2*J+2, J + (1:J)] <- cov.Im.Y2

  SIG[2*J+1, 2*J+1] <- var.Y1
  SIG[2*J+2, 2*J+2] <- var.Y2

  SIG <- (SIG + t(SIG))/2
  SIG <- SIG + diag(1e-10, p)

  T.vals <- numeric(B)

  for (b in 1:B) {

    z <- mvrnorm(1, mu = rep(0, p), Sigma = SIG)

    Gre <- z[1:J]
    Gim <- z[(J+1):(2*J)]
    Y1  <- z[2*J + 1]
    Y2  <- z[2*J + 2]
    G   <- Gre + 1i * Gim

    G0 <- G - 1i * t * phi0 * Y1 + (t^2) * phi0 * Y2

    Qm <- numeric(M)

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

      G0.u1 <- interp_complex(t, G0, u1)
      G0.u2 <- interp_complex(t, G0, u2)

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

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

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

    T.vals[b] <- max(abs((Qm - mu) / sigma))
  }

  T.vals
}

# Example generation (not run):
# calib.normal  <- calibrate_mu_sigma(weight.type = "normal")
# T.normal <- simulate_T(mu = calib.normal$mu, sigma = calib.normal$sigma,
#                          weight.type = "normal")