# Load parallel computing package
library(parallel)
library(MASS) 

LCR <- function(S, q, r) {
  
  u <- rbinom(1, 1, r)
  v <- rbinom(1, 1, 0.5)
  
  if (u == 1) {
    return(as.numeric(q > S))
  } else {
    return(v)
  }
}
leixing='setting1'
n <- 10000
p <- 5
alpha <- 0.1
num_repeats <- 200  # Number of repetitions
r=0.25



beta <- matrix(0, nrow = 3, ncol = p)
beta[1, ] <- c(1, 2, 1, 0, 0)       
beta[2, ] <- c(0, -1, -2, -1, 0)    
beta[3, ] <- c(0, 0, 1, 2, 1)


# Define a function to run one simulation
run_simulation <- function(repeat_idx) {
  set.seed(202504+repeat_idx)  # Set seed for reproducibility
  # Generate data
  X <- mvrnorm(n, mu = rep(0, 5), Sigma = diag(5))
  epsilon <- rnorm(n)
  Y_vals <- numeric(n)
  Y_vals[1:3500] <- X[1:3500, ] %*% beta[1,] + epsilon[1:3500]
  Y_vals[3501:6500] <- X[3501:6500, ] %*% beta[2,] + epsilon[3501:6500]
  Y_vals[6501:10000] <- X[6501:10000, ] %*% beta[3,] + epsilon[6501:10000]
  
  beta_hat <- rep(0, p)  
  stepsize <- 0.01   
  Y_hat <- numeric(n)  
  s <- numeric(n)  
  W <- numeric(n)  
  lambda <- numeric(n)
  lambda[1] <- 0  
  
  s[1] <- 0
  # Main loop to iterate over all data points
  for (i in 1:n) {
    x_i <- X[i, ]
    y_i <- Y_vals[i]
    
    
    Y_hat[i] <- sum(x_i * beta_hat)
    
    
    r_hat <- y_i - Y_hat[i]
    
    
    grad <- -r_hat * x_i
    beta_hat <- beta_hat - stepsize * grad
    
    
    
    S_t <- abs(Y_vals[i] - Y_hat[i])
    L=LCR(S_t,s[i],r)
    
    if(L==1)
    {g_t=1-(r*(1-alpha)+(1-r)*0.5)}
    else
    {g_t=-(r*(1-alpha)+(1-r)*0.5)}
    
    
    W[i] <- ifelse(i == 1, 1 - g_t * s[1], W[i-1] - g_t * s[i])
    lambda[i + 1] <- (i) / (i + 1) * lambda[i] - 1 / (i + 1) * g_t
    s[i + 1] <- lambda[i + 1] * W[i]
  }
  
  
  
  prediction_intervals <- data.frame(
    Lower = Y_hat[1:n] - s[1:n],
    Upper = Y_hat[1:n] + s[1:n],
    Y_true = Y_vals
  )
  coverage <- ifelse(prediction_intervals$Y_true >= prediction_intervals$Lower & 
                       prediction_intervals$Y_true <= prediction_intervals$Upper, 1, 0)
  return(list(coverage = coverage, s = s[1:n])) 
}


num_cores <- detectCores() - 1 
cl <- makeCluster(num_cores)


clusterExport(cl, c("n", "p", "alpha",  "beta",  "run_simulation",'mvrnorm','LCR','r'))


results <- parLapply(cl, 1:num_repeats, run_simulation)


stopCluster(cl)


coverage_matrix <- sapply(results, function(res) res$coverage)
s_matrix <- sapply(results, function(res) res$s)

average_coverage <- rowMeans(coverage_matrix)
average_s <- 2 * rowMeans(s_matrix)  # 区间长度 = 2*s


burn_in <- 200
cumulative_coverage_matrix <- apply(coverage_matrix, 2, function(x) {
  x_sub <- x[(burn_in+1):length(x)]
  cumsum(x_sub) / seq_along(x_sub)
})
cumulative_coverage <- rowMeans(cumulative_coverage_matrix)

interval_matrix <- 2 * s_matrix
cumulative_interval_matrix <- apply(interval_matrix, 2, function(x) {
  x_sub <- x[(burn_in+1):length(x)]
  cumsum(x_sub) / seq_along(x_sub)
})
cumulative_interval <- rowMeans(cumulative_interval_matrix)


overall_coverage_per_rep <- colMeans(coverage_matrix[(burn_in+1):n, ])
overall_interval_per_rep <- colMeans(interval_matrix[(burn_in+1):n, ])



sd_coverage <- sd(overall_coverage_per_rep)
sd_interval <- sd(overall_interval_per_rep)

cumulative_coverage <- c(rep(0, burn_in), cumulative_coverage)
cumulative_interval <- c(rep(0, burn_in), cumulative_interval)

# 4. 创建包含所有统计量的数据框
tag <- paste0("model", leixing, "r", r)
output_data <- data.frame(
  InstantCoverage = average_coverage,        # 瞬时覆盖率
  InstantIntervalLength = average_s,         # 瞬时区间长度
  CumulativeCoverage = cumulative_coverage,  # 累积覆盖率
  CumulativeIntervalLength = cumulative_interval,  # 累积区间长度
  OverallSDCoverage = sd_coverage,           # 整体覆盖率标准差
  OverallSDInterval = sd_interval            # 整体区间长度标准差
)


setwd('D:/desktop/randomize/Simulations_4.16')
dir.create("results", showWarnings = FALSE)
write.csv(output_data, paste0("results/", tag, ".csv"), row.names = FALSE)
