在R中为每个ID抽取单条记录，匹配目标date分布的方法

我有一份数据da，其date字段存在特定的频率分布（如下所示）。另有一份数据db，其中每个id可能包含一条或多条记录。请问是否存在可行方法，为每个id抽取且仅抽取一条记录，使得db中date字段的抽样分布尽可能接近da中的date分布？

示例数据代码 #

library(data.table)
library(dplyr)
library(lubridate)

# 目标分布数据da
da = data.table(id = paste0('a', 1:7),
                date = ymd(c('2021-1-10', 
                             rep('2021-1-11', 2), 
                             rep('2021-1-12', 3), 
                             '2021-1-13')))

da
# 查看da的date字段频率分布
da[,.N,(date)][,.(date, N, perc = N/sum(N))]
#          date N      perc
# 1: 2021-01-10 1 0.1428571
# 2: 2021-01-11 2 0.2857143
# 3: 2021-01-12 3 0.4285714
# 4: 2021-01-13 1 0.1428571


# 需要为每个id仅抽取一条记录，使得date分布匹配da的分布
set.seed(123)
db = structure(list(id = c(1L, 2L, 3L, 3L, 3L, 4L, 5L, 6L, 6L, 8L), 
               date = structure(c(18638, 18639, 18639, 18640, 18640, 18637, 
                                  18640, 18637, 18638, 18639), class = "Date")), 
          class = c("data.table", "data.frame"))

# 查看db数据
db
#    id       date
# 1:  1 2021-01-11
# 2:  2 2021-01-12
# 3:  3 2021-01-12
# 4:  3 2021-01-13
# 5:  3 2021-01-13
# 6:  4 2021-01-10
# 7:  5 2021-01-13
# 8:  6 2021-01-10
# 9:  6 2021-01-11
#10:  8 2021-01-12

这是一个约束优化下的抽样匹配问题，根据对匹配精度的需求，可选择以下两种方法：

1. 加权随机抽样（近似匹配） #

这种方法简单易实现，通过给每个id的候选日期分配与目标分布成正比的权重，随机抽取符合权重偏向的记录，能快速得到近似目标分布的结果。

# 提取da的目标日期分布权重
target_dist = da[, .(target_perc = .N / nrow(da)), by = date]

# 将目标权重合并到db中，db中存在但da中没有的日期权重设为0
db_with_weights = merge(db, target_dist, by = "date", all.x = TRUE)
db_with_weights[is.na(target_perc), target_perc := 0]

# 按id分组，基于目标权重抽样
set.seed(456)
sampled_db = db_with_weights[, .SD[sample(.N, 1, prob = target_perc)], by = id]

# 查看抽样后的分布
sampled_db[, .N, by = date][, .(date, N, perc = N / nrow(sampled_db))]

2. 整数规划（精确匹配） #

如果需要尽可能精准地匹配目标分布，可将问题转化为整数规划，通过约束条件最小化实际分布与目标分布的差异，适合对匹配精度要求高的场景。

library(lpSolve)

# 整理候选记录的id、date信息
records = db[, .(record_id = .I, id, date)]
# 计算每个日期需要选中的目标数量（基于db总记录数和da的分布比例取整）
target_counts = target_dist[, .(target_count = round(nrow(db)*target_perc)), by = date]

# 构建约束矩阵：每个id仅选1条；每个日期选中数接近目标数
id_constraints = model.matrix(~ factor(id) - 1, data = records)
date_constraints = model.matrix(~ factor(date) - 1, data = records)
constraint_matrix = rbind(id_constraints, date_constraints)

# 设置约束方向和约束值
constraint_dir = c(rep("=", nrow(unique(records[,.(id)]))), rep("=", nrow(target_counts)))
constraint_rhs = c(rep(1, nrow(unique(records[,.(id)]))), target_counts$target_count)

# 定义目标函数：最小化实际分布与目标分布的偏差平方
objective = rep(0, nrow(records))
for (d in target_counts$date) {
  idx = records$date == d
  objective[idx] = (1/target_counts[date == d, target_count])^2
}

# 求解整数规划
lp_result = lp("min", objective, constraint_matrix, constraint_dir, constraint_rhs, binary.vec = rep(1, nrow(records)))

# 提取选中的记录
sampled_exact = records[lp_result$solution == 1, ]
sampled_exact = merge(sampled_exact, db, by = c("id", "date"))

# 查看精确匹配后的分布
sampled_exact[, .N, by = date][, .(date, N, perc = N / nrow(sampled_exact))]

内容的提问来源于stack exchange，提问作者Miao Cai