如何在R中实现拒收抽样？
我有一个每行基因的数据集，每行都有它们的基因长度，我希望使用排斥抽样从这些基因中通过它们的基因长度分布进行采样-因为我在这个数据集中有太多太小的基因，无法进入进一步的分析(但我不能自己设置一个截止点来删除它们)。我有一个带有基因长度的基因数据集可供采样，还有一个建议的基因长度分布，我想使用它来对第一个数据集进行拒绝采样。
我的数据示例如下所示：
#df1 data to sample from:
Gene Length
Gene1 5
Gene2 6
Gene3 400000
Gene4 1000
Gene5 25000
Gene6 10
Gene7 50
Gene8 4
Gene9 100
Gene10 2000我的提案数据集：
#df2
Gene Length
Gene1 5000
Gene2 60000
Gene3 400000
Gene4 1000
Gene5 25000
Gene6 10000
Gene7 50000
Gene8 4000
Gene9 1000
Gene10 2000我没有任何统计学背景，我正在尝试进行排斥抽样(我的总体目标是获得长度较小的基因样本，以便进行进一步分析)。
为了做拒绝采样，我从我发现的here教程中尝试了这个
X = df1$Length
U = df2$Length
accept = c()
count = 1
pi_x <- function(x) {
new_x = (3/2)*(x^3)+(11/8)*(x^2)+(1/6)*(x)+(1/12)
return(new_x)
}
while(count <= 50 & length(accept) < 50){
test_u = U[count]
test_x = pi_x(X[count])/(3.125*dunif(X[count],0,1))
if (test_u <= test_x){
accept = rbind(accept, X[count])
count = count + 1
}
count = count + 1
}我的问题是，它只选择了25个基因(我进一步分析的理想采样范围是选择50-100个基因)，并且这25个基因中的大多数在采样后仍然非常小。在运行这个拒绝采样代码之前，我需要以某种方式转换X吗？我的df1的实际数据是800个基因，基因长度呈倾斜/β分布(大多数都很小)。或者我完全忽略了我的理解中的其他东西？任何指导都将不胜感激。
输入数据：
df1 <- structure(list(Gene = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5",
"Gene6", "Gene7", "Gene8", "Gene9", "Gene10"), Length = c(5L,
6L, 400000L, 1000L, 25000L, 10L, 50L, 4L, 100L, 2000L)), row.names = c(NA,
-10L), class = c("data.table", "data.frame"))
df2 <- structure(list(Gene = c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5",
"Gene6", "Gene7", "Gene8", "Gene9", "Gene10"), Length = c(5000L,
60000L, 400000L, 1000L, 25000L, 10000L, 50000L, 40000L, 1000L, 2000L)), row.names = c(NA,
-10L), class = c("data.table", "data.frame"))编辑：
我也尝试过：
sampled <- data.frame(proposal = df2$Length)
sampled$targetDensity <- dbeta(sampled$proposal, 3,6)
maxDens = max(sampled$targetDensity, na.rm = T)
sampled$accepted = ifelse(df1$Length < sampled$targetDensity / maxDens, TRUE, FALSE)
hist(sampled$proposal[sampled$accepted], freq = F, col = "grey", breaks = 100)
curve(dbeta(x, 3,6),0,1, add =T, col = "red")但我确信我没有正确使用dbeta()，因为sampled$targetDensity输出的全是零--有没有办法解决这个问题？我曾尝试更改dbeta()中的值，但没有任何成功。
回答 1
Stack Overflow用户
发布于 2021-05-26 22:57:05
如果您知道要采样的基因数量，则采样函数应该可以很好地工作：
sampled = sample(df$genes, size = n, prob = df$length) 如果您想进一步降低对长度较小的基因进行采样的概率，可以将prob参数的长度平方。参数prob将采样概率关联到每个元素(这里是基于长度)
如果你不知道你想要获得的基因数量，那么你可以定义自己的概率函数：
rejection.prob = function(x){
if (x<too.small) {return(0)} # all genes smaller than too.small won't be sampled
if (x > it.is.ok) {return(1)} # all these ones will be sampled
if (x>too.small & (x<it.is.ok){
# here make any function that is equal to 0 when x == too.small
# and 1 when x == it.is.ok
# it can be a simple linear function
}请注意，您还可以将rejection.prob函数的输出用于sample函数。
根据您的期望，您可能希望您的拒绝函数更加连续(在too.small和it.is.ok中没有这些中断)。如果是这种情况，我会使用sigmoid函数，您可以根据需要的输出调整参数。
https://stackoverflow.com/questions/67671330
复制相似问题
腾讯云开发者
