# 数据初始化
last_year <- c(33, 36, 71, 48, 19, 69, 48)  # 去年底的商店数
current_year <- c(33, 39, 63, 51, 23, 77, 45)  # 现有商店数
n <- 2  # 抽取的样本数

# 计算累计和
cumsum_last_year <- cumsum(last_year)  # 去年底商店数的累计和

# 给定的随机数
r <- c(0.54463226, 0.15389852)

# 使用累计和方法抽取样本单元
sample_units <- sapply(r, function(r_val) {
  which(cumsum_last_year >= r_val * sum(last_year))[1]
})

# 输出抽取的样本单元
cat("抽取的样本单元:", sample_units, "\n")

# 获取抽取的样本的现有商店数和去年底商店数
sample_current <- current_year[sample_units]
sample_last <- last_year[sample_units]

# 计算HT估计值（Horvitz-Thompson Estimate）
ht_estimate <- sum(sample_current / (sample_last / sum(last_year)))
cat("HT估计值 Ŷ_HT:", ht_estimate, "\n")

# 计算样本的权重（即每个样本的抽样概率的倒数）
weights <- sum(last_year) / sample_last

# 计算方差的估计
# 首先计算每个样本的贡献
contribution <- (sample_current / sample_last - ht_estimate / sum(last_year))^2

# 估计的方差
variance <- sum(contribution * weights) / (n - 1)

# 计算标准差的估计值
std_dev_hat <- sqrt(variance)
cat("标准差的估计值:", std_dev_hat, "\n")

# 计算均方偏差（MSE）
mse_estimate <- variance
cat("均方偏差（MSE）:", mse_estimate, "\n")