我有一个基因行数据集,每个基因行都有它们的基因长度,我希望使用拒绝抽样通过基因长度分布从这些基因中取样 - 因为我在这个数据集中有太多基因太小而无法进入进一步分析(但是我不能自己设置一个截止点来删除它们)。我有一个基因长度的基因数据集可供采样,还有另一个基因长度的提议分布,我想使用它来对第一个数据集进行拒绝采样。
我的数据示例如下所示:
#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
我没有任何统计背景,我正在尝试进行拒绝抽样(我的总体目标是获取长度极小基因较少的基因样本以进行进一步分析)。
要进行拒绝抽样,我正在从我在这里找到的教程中尝试这个:
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()
但没有任何成功。