在寻找解决另一个问题的方法时,我遇到了一个问题,即base::sample()
函数在不进行替换的情况下进行采样时,似乎会截断小概率(尽管不是零概率)。
使用下面的代码,我从40 - 10,000次中提取了10个值。注意:当用100,000运行时,每次需要3分钟,而不是〈3秒。
N_simulations = 10000
N_draw_per_sim = 10
dat = data.frame(id = 1:40,
log_likelihood = seq(from = 550, to = 350, length.out = 40))
dat$lkhd_wt = (function(x) {x / sum(x)}) (exp(dat$log_likelihood - max(dat$log_likelihood)))
sample_base = function(N_) {sample(x = dat$id, size = N_, prob = dat$lkhd_wt, replace = FALSE)}
sample_manual = function(N_) {
# manually sample once, remove drawn value, re-weight
smpl_ = sample(x = dat$id, size = 1, prob = (function(x) {x / sum(x)}) (dat$lkhd_wt), replace = FALSE)
for(i in 2:N_) {
prv_idx = which(dat$id %in% smpl_)
smpl_ = c(smpl_, sample(x = dat$id[-prv_idx], size = 1, prob = (function(x) {x / sum(x)}) (dat$lkhd_wt[-prv_idx]), replace = FALSE))
}
return(smpl_)
}
# 10k < 2 sec
dat_sim_draw_basic = as.data.frame(matrix(0, nrow = dim(dat)[1], ncol = N_simulations))
for(sim_i in 1:N_simulations) {dat_sim_draw_basic[sample_base(N_draw_per_sim), sim_i] = 1}
# 10k < 3 sec
dat_sim_draw_manual = as.data.frame(matrix(0, nrow = dim(dat)[1], ncol = N_simulations))
for(sim_i in 1:N_simulations) {dat_sim_draw_manual[sample_manual(N_draw_per_sim), sim_i] = 1}
# Calculate weights from the count
dat_agg_basic = data.frame(id = dat$id, cnt = apply(dat_sim_draw_basic, 1, sum))
dat_agg_basic$wt = dat_agg_basic$cnt / sum(dat_agg_basic$cnt)
dat_agg_manual = data.frame(id = dat$id, cnt = apply(dat_sim_draw_manual, 1, sum))
dat_agg_manual$wt = dat_agg_manual$cnt / sum(dat_agg_manual$cnt)
# Merge for comparison
TF_non_zero = dat_agg_basic$wt != 0 | dat_agg_manual$wt != 0
dat_compare = merge(x = dat_agg_basic[TF_non_zero,],
y = dat_agg_manual[TF_non_zero,],
by.x = c("id"),
by.y = c("id"),
all.x = TRUE,
all.y = TRUE)
colnames(dat_compare) = c("id", "cnt_basic", "wt_basic", "cnt_manual", "wt_manual")
dat_compare = dat_compare[,c("id", "cnt_basic", "cnt_manual", "wt_basic", "wt_manual")]
dat_compare
输出(i.e. dat_compare
)
id cnt_basic cnt_manual wt_basic wt_manual
1 1 10000 10000 0.1 0.10000
2 2 10000 10000 0.1 0.10000
3 3 10000 10000 0.1 0.10000
4 4 10000 10000 0.1 0.10000
5 5 10000 10000 0.1 0.10000
6 6 10000 10000 0.1 0.10000
7 7 10000 10000 0.1 0.10000
8 8 10000 10000 0.1 0.10000
9 9 10000 10000 0.1 0.10000
10 10 10000 9958 0.1 0.09958
11 11 0 42 0.0 0.00042
这似乎表明传递给sample()
的权重没有全部使用。我的猜测是,小的权重被截断为零,并且在绘制后从不重新计算。
我尝试将权重乘以不同的幂(10^i,其中i为1:86),并计算等于0的权重的数量。1. sample()
函数有10个权重等于0。1.然而,情节古怪;我本以为它会被分成不同的区域,但它似乎只是跳得到处都是。
这让我彻底糊涂了。有人知道为什么这是,一种迫使它重新计算的方法,或者一种绕过这个问题的方法吗?除了写一个函数来一次绘制一个并重新计算。
2条答案
按热度按时间bvjxkvbb1#
简而言之,您遇到的是计算中常见的floating point issue。
你的概率非常低(最多87个零- i。即
1.3e-87
)。从纯概率的Angular 来看,在10,000次抽奖中,你仍然极不可能看到它。但这也远远低于大多数机器所能处理的。不需要太多的细节,您可以检查机器上的最小正浮点数。从R 4.0开始,你可以检查你的机器是否支持长度超过double的
C long double
类型。我的机器支持长双
如果你运行
.Machine
,你会看到你的机器的规格:我可以运行到大约20位小数,这对应于你在概率向量
dat$lkhd_wt
中的第**~ 10**个位置。从Limey那里偷了一点代码来快速看看这是如何实现的:
hvvq6cgz2#
你观察到的问题可能是由于R中浮点数的有限数值精度造成的。传递给
sample()
函数的权重非常小(大约为10^-86),因此,它们在计算机内存中的表示可能会被截断或舍入,从而导致数值错误。解决此问题的一个可行方法是使用
rmultinom()
函数而不是sample()
来生成示例。rmultinom()
函数从多项分布生成随机样本,可用于模拟从具有给定概率的总体中绘制样本而不进行替换。以下是如何使用rmultinom()
生成示例的示例:在此修改后的代码中,
sample_multinom()
函数使用rmultinom()
函数生成多项式样本,将其转换为索引向量,然后使用sample()
函数随机置换索引向量以获得最终样本。代码的其余部分与原始代码类似,不同之处在于它还计算和比较使用rmultinom()
获得的权重。使用
rmultinom()
而不是sample()
应该可以消除您在小概率情况下观察到的数值精度问题,并且在生成样本时也应该比sample_manual()
函数更有效。