## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)

## ----message=FALSE, warning=FALSE---------------------------------------------
library(RegimeChange)

# Generate example data
set.seed(123)
data <- c(rnorm(100, 0, 1), rnorm(100, 3, 1))

# Basic BOCPD
result <- detect_regimes(data, method = "bocpd")
print(result)

## -----------------------------------------------------------------------------
# Define prior
prior <- normal_gamma(
  mu0 = 0,      # Prior mean for the mean
  kappa0 = 1,   # Prior precision for the mean
  alpha0 = 1,   # Shape for inverse-gamma on variance
  beta0 = 1     # Rate for inverse-gamma on variance
)

print(prior)

## -----------------------------------------------------------------------------
result <- detect_regimes(data, method = "bocpd", prior = prior)

## -----------------------------------------------------------------------------
prior_known <- normal_known_var(
  mu0 = 0,         # Prior mean
  sigma0 = 1,      # Prior standard deviation for mean
  known_var = 1    # Known variance
)

result <- detect_regimes(data, method = "bocpd", prior = prior_known)

## -----------------------------------------------------------------------------
prior_poisson <- poisson_gamma(alpha0 = 1, beta0 = 1)

## -----------------------------------------------------------------------------
# Expected run length = 1/lambda = 100
hazard <- geometric_hazard(lambda = 0.01)
print(hazard)

## -----------------------------------------------------------------------------
hazard_const <- constant_hazard(lambda = 0.05)

## -----------------------------------------------------------------------------
hazard_nb <- negbin_hazard(r = 5, p = 0.1)

## -----------------------------------------------------------------------------
result <- detect_regimes(data, method = "bocpd")

# Access changepoint probabilities
prob_change <- result$prob_change

# Plot posterior probability
plot(prob_change, type = "l", 
     main = "Posterior Probability of Changepoint",
     xlab = "Time", ylab = "P(changepoint)")
abline(v = 100, col = "red", lty = 2)

## -----------------------------------------------------------------------------
plot(result, type = "posterior")

## -----------------------------------------------------------------------------
plot(result, type = "runlength")

## -----------------------------------------------------------------------------
result_sr <- detect_regimes(data, method = "shiryaev",
                            mu0 = 0,      # Pre-change mean
                            mu1 = 3,      # Post-change mean
                            sigma = 1,    # Known standard deviation
                            threshold = 100)
print(result_sr)

## -----------------------------------------------------------------------------
plot(result_sr, type = "diagnostic")

## -----------------------------------------------------------------------------
# Highest posterior density interval
if (length(result$changepoints) > 0) {
  cat("Changepoint estimate:", result$changepoints[1], "\n")
  
  # Find 95% credible interval
  prob <- result$prob_change
  cumprob <- cumsum(prob) / sum(prob)
  lower <- which(cumprob >= 0.025)[1]
  upper <- which(cumprob >= 0.975)[1]
  cat("95% Credible interval: [", lower, ",", upper, "]\n")
}

## -----------------------------------------------------------------------------
existence_prob <- max(result$prob_change)
cat("Maximum changepoint probability:", round(existence_prob, 3), "\n")

## -----------------------------------------------------------------------------
# Create online detector
detector <- regime_detector(
  method = "bocpd",
  prior = normal_gamma(),
  hazard = geometric_hazard(0.01),
  threshold = 0.5
)

# Simulate online processing
set.seed(123)
stream <- c(rnorm(100, 0, 1), rnorm(50, 3, 1))

alarm_time <- NULL
for (i in seq_along(stream)) {
  detector <- update(detector, stream[i])
  
  if (isTRUE(detector$last_result$alarm) && is.null(alarm_time)) {
    alarm_time <- i
    cat("Alarm at observation:", i, "\n")
    cat("Probability of change:", round(detector$last_result$prob_change, 3), "\n")
    break
  }
}

## -----------------------------------------------------------------------------
# Generate bivariate data
set.seed(42)
n <- 200
cp <- 100
data_mv <- rbind(
  matrix(rnorm(cp * 2, 0), ncol = 2),
  matrix(rnorm((n - cp) * 2, 2), ncol = 2)
)

# Normal-Wishart prior for multivariate data
prior_mv <- normal_wishart(
  mu0 = c(0, 0),
  kappa0 = 1,
  nu0 = 4,
  Psi0 = diag(2)
)

result_mv <- detect_regimes(data_mv, method = "bocpd", prior = prior_mv)
print(result_mv)

## -----------------------------------------------------------------------------
# Compare BOCPD with PELT
comparison <- compare_methods(
  data = data,
  methods = c("bocpd", "pelt"),
  true_changepoints = 100
)
print(comparison)