R的基本样本函数(无替换)忽略小概率

tsm1rwdh  于 2023-05-04  发布在  其他
关注(0)|答案(2)|浏览(159)

在寻找解决另一个问题的方法时,我遇到了一个问题,即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.然而,情节古怪;我本以为它会被分成不同的区域,但它似乎只是跳得到处都是。

这让我彻底糊涂了。有人知道为什么这是,一种迫使它重新计算的方法,或者一种绕过这个问题的方法吗?除了写一个函数来一次绘制一个并重新计算。

bvjxkvbb

bvjxkvbb1#

简而言之,您遇到的是计算中常见的floating point issue
你的概率非常低(最多87个零- i。即1.3e-87)。从纯概率的Angular 来看,在10,000次抽奖中,你仍然极不可能看到它。但这也远远低于大多数机器所能处理的。
不需要太多的细节,您可以检查机器上的最小正浮点数。从R 4.0开始,你可以检查你的机器是否支持长度超过double的C long double类型。

capabilities()["long.double"]
# long.double 
#        TRUE

我的机器支持长双
如果你运行.Machine,你会看到你的机器的规格:

.Machine$double.eps
# [1] 2.220446e-16

.Machine$longdouble.neg.eps
# [1] 5.421011e-20

我可以运行到大约20位小数,这对应于你在概率向量dat$lkhd_wt中的第**~ 10**个位置。
从Limey那里偷了一点代码来快速看看这是如何实现的:

hist(sapply(1:1e4, 
            function(x) sample(1:40, 10, prob = dat$lkhd_wt)), breaks=0:40, freq=FALSE)

hvvq6cgz

hvvq6cgz2#

你观察到的问题可能是由于R中浮点数的有限数值精度造成的。传递给sample()函数的权重非常小(大约为10^-86),因此,它们在计算机内存中的表示可能会被截断或舍入,从而导致数值错误。
解决此问题的一个可行方法是使用rmultinom()函数而不是sample()来生成示例。rmultinom()函数从多项分布生成随机样本,可用于模拟从具有给定概率的总体中绘制样本而不进行替换。以下是如何使用rmultinom()生成示例的示例:

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_multinom <- function(N_) {
  # generate a multinomial sample with given probabilities
  smpl_ <- rmultinom(N_, size = 1, prob = dat$lkhd_wt)
  # convert the multinomial sample to an index vector
  idx_ <- rep(dat$id, N_)[as.logical(smpl_)]
  # randomly permute the index vector to get the final sample
  return(sample(idx_, size = N_, replace = FALSE))
}

dat_sim_draw_multinom <- as.data.frame(matrix(0, nrow = dim(dat)[1], ncol = N_simulations))
for (sim_i in 1:N_simulations) {
  dat_sim_draw_multinom[sample_multinom(N_draw_per_sim), sim_i] <- 1
}

dat_agg_multinom <- data.frame(id = dat$id, cnt = apply(dat_sim_draw_multinom, 1, sum))
dat_agg_multinom$wt <- dat_agg_multinom$cnt / sum(dat_agg_multinom$cnt)

dat_compare <- merge(x = dat_agg_basic[TF_non_zero,],
                     y = dat_agg_manual[TF_non_zero,],
                     z = dat_agg_multinom[TF_non_zero,],
                     by.x = "id", by.y = "id", by.z = "id",
                     all.x = TRUE, all.y = TRUE, all.z = TRUE)

colnames(dat_compare) <- c("id", "cnt_basic", "cnt_manual", "cnt_multinom", "wt_basic", "wt_manual", "wt_multinom")
dat_compare <- dat_compare[, c("id", "cnt_basic", "cnt_manual", "cnt_multinom", "wt_basic", "wt_manual", "wt_multinom")]
dat_compare

在此修改后的代码中,sample_multinom()函数使用rmultinom()函数生成多项式样本,将其转换为索引向量,然后使用sample()函数随机置换索引向量以获得最终样本。代码的其余部分与原始代码类似,不同之处在于它还计算和比较使用rmultinom()获得的权重。
使用rmultinom()而不是sample()应该可以消除您在小概率情况下观察到的数值精度问题,并且在生成样本时也应该比sample_manual()函数更有效。

相关问题