用多个概率分布模拟 R 中的数据答案

【问题标题】：Simulating data in R with multiple probability distributions用多个概率分布模拟 R 中的数据
【发布时间】：2016-02-03 13:25:27
【问题描述】：

我正在尝试通过自举来模拟数据，以使用漏斗图为我的真实数据创建置信带。我正在建立已接受答案to a previous question 的策略。我不想使用单个概率分布来模拟我的数据，而是想修改它以根据被模拟的数据部分使用不同的概率分布。

我非常感谢任何可以帮助回答问题或帮助我更清楚地表达问题的人。

我的问题是编写适当的 R 代码来进行更复杂的数据模拟。

目前的代码是：

n <- 1e4
set.seed(42)
sims <- sapply(1:80, 
               function(k) 
                 rowSums(
                   replicate(k, sample((1:7)/10, n, TRUE, ps))) / k)

此代码模拟数据，其中每个数据点都有一个值，该值是1:80 观察值之间的平均值。例如，当数据点的值是 10 个观测值的平均值 (k=10) 时，它会根据概率分布ps，给出了每个值的概率（基于整个经验分布）。

ps 看起来像这样：

ps <- prop.table(table((DF$mean_score)[DF$total_number_snps == 1]))
#        0.1         0.2         0.3         0.4         0.5         0.6         0.7 
#0.582089552 0.194029851 0.124378109 0.059701493 0.029850746 0.004975124 0.004975124

例如，观察值为0.1 的概率是0.582089552。

现在，我希望根据每个数据点的观察数量有条件地使用不同的频率分布，而不是对所有模拟使用一个频率分布。

我制作了一个表格，cond_probs，其中包含我的每个真实数据点的一行。有一列包含total 的观察次数，还有一列给出了每个观察值的频率。

cond_probs 表示例：

gene_name   0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 total
A1  0.664   0.319   0.018   0.000   0.000   0.000   0.000   0.000   0.000   113.000
A2  0.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   1.000

所以对于数据点A2，只有1观察，其值为0.1。因此0.1 观察的频率是1。对于A1，有113 观察值，其中大多数(0.664) 的值为0.1。这个想法是cond_probs 类似于ps，但cond_probs 对每个数据点都有一个概率分布，而不是所有数据的一个概率分布。

我想修改上面的代码，以便将采样修改为使用cond_probs 而不是ps 进行频率分布。并使用观察次数 k 作为选择从 cond_probs 中的哪一行进行采样的标准。所以它会像这样工作：

对于具有k 观察次数的数据点：

查看cond_probs 表并随机选择其中total 观察数与k 大小相似的行：0.9k-1.1k。如果不存在这样的行，请继续。

一旦选择了一个数据点，就使用cond_probs 中该行的概率分布，就像在原始代码中使用ps 一样，随机抽样k 的观察次数并输出这些观察的平均值。

对于replicate 的每个n 迭代，在total 的值与k 的当前值相似的所有行中，随机抽样并替换来自cond_probs 的新数据点（0.9k-1.1k）。

这个想法是，对于这个数据集，应该根据数据点的观察数量来确定要使用的概率分布。这是因为在该数据集中，观察的概率受观察数量的影响（由于遗传连锁和背景选择，具有更多 SNP 的基因在每次观察中的得分往往较低）。

使用下面的答案更新：

我尝试使用下面的答案，它适用于示例中的模拟 cond_probs 数据，但不适用于我的真实 cond_probs 文件。我将我的 cond_probs 文件导入并转换为带有

的矩阵

cond_probs <- read.table("cond_probs.txt", header = TRUE, check.names = FALSE)
cond_probs <- as.matrix(cond_probs)

第一个示例十行（约 20,000 行）如下所示：

>cond_probs
       total   0.1   0.2   0.3   0.4   0.5   0.6   0.7   0.8   0.9   1.0
[1,]     109 0.404 0.174 0.064 0.183 0.165 0.009 0.000 0.000 0.000 0.000
[2,]     181 0.564 0.221 0.144 0.066 0.006 0.000 0.000 0.000 0.000 0.000
[3,]     289 0.388 0.166 0.118 0.114 0.090 0.093 0.028 0.003 0.000 0.000
[4,]     388 0.601 0.214 0.139 0.039 0.008 0.000 0.000 0.000 0.000 0.000
[5,]     133 0.541 0.331 0.113 0.000 0.008 0.008 0.000 0.000 0.000 0.000
[6,]     221 0.525 0.376 0.068 0.032 0.000 0.000 0.000 0.000 0.000 0.000
[7,]     147 0.517 0.190 0.150 0.054 0.034 0.048 0.007 0.000 0.000 0.000
[8,]     107 0.458 0.196 0.252 0.084 0.009 0.000 0.000 0.000 0.000 0.000
[9,]      13 0.846 0.154 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

如果我跑：

sampleSize <- 20
set.seed(42)
#replace 1:80 with 1: max number of SNPs in gene in dataset
sims_test <- sapply( 1:50, simulateData, sampleSize )

并查看具有 x 次观察的抽样均值，我只得到一个结果，而应该有 20 个。

例如：

> sims_test[[31]]
[1] 0.1

并且sims_test的排序方式与sims不同：

>sims_test
   [,1] [,2]      [,3]  [,4] [,5]      [,6]      [,7]   [,8]      [,9]
 [1,]  0.1  0.1 0.1666667 0.200 0.14 0.2666667 0.2000000 0.2375 0.1888889
 [2,]  0.1  0.1 0.1333333 0.200 0.14 0.2333333 0.1571429 0.2625 0.1222222
 [3,]  0.1  0.1 0.3333333 0.225 0.14 0.1833333 0.2285714 0.2125 0.1555556
 [4,]  0.1  0.1 0.2666667 0.250 0.10 0.1500000 0.2000000 0.2625 0.2777778
 [5,]  0.1  0.1 0.3000000 0.200 0.16 0.2000000 0.2428571 0.1750 0.1000000
 [6,]  0.1  0.1 0.3666667 0.250 0.16 0.1666667 0.2142857 0.2500 0.2000000
 [7,]  0.1  0.1 0.4000000 0.300 0.12 0.2166667 0.1857143 0.2375 0.1666667
 [8,]  0.1  0.1 0.4000000 0.250 0.10 0.2500000 0.2714286 0.2375 0.2888889
 [9,]  0.1  0.1 0.1333333 0.300 0.14 0.1666667 0.1714286 0.2750 0.2888889

更新 2

使用 cond_probs <- head(cond_probs,n) 我已经确定代码在 n = 517 之前有效，然后对于所有大于此的大小，它都会产生与上述相同的输出。我不确定这是文件本身的问题还是内存问题。我发现如果我删除第 518 行并将之前的行复制几次以制作更大的文件，它可以工作，这表明该行本身导致了问题。第 518 行如下所示：

9.000   0.889   0.000   0.000   0.000   0.111   0.000   0.000   0.000   0.000   0.000

我发现另外 4 条违规行：

9.000   0.444   0.333   0.111   0.111   0.000   0.000   0.000   0.000   0.000   0.000

9.000   0.444   0.333   0.111   0.111   0.000   0.000   0.000   0.000   0.000   0.000

9.000   0.111   0.222   0.222   0.111   0.111   0.222   0.000   0.000   0.000   0.000

9.000   0.667   0.111   0.000   0.000   0.000   0.222   0.000   0.000   0.000   0.000

我没有注意到他们有什么不寻常的地方。他们都有 9 个站点。如果我删除这些行并运行仅包含这些行之前的“cond_probs”文件，那么代码就可以工作。但是肯定还有其他有问题的行，因为整个 'cond_probs' 仍然不起作用。

我尝试将这些有问题的行放回一个较小的“cond_probs”文件中，然后这个文件就可以工作了，所以我很困惑，因为这些行似乎没有本质上的问题。另一方面，它们共有 9 个站点，这表明某种原因模式。

如果有帮助的话，我很乐意私下分享整个文件，因为我不知道下一步该怎么做才能进行故障排除。

出现的另一个问题是我不确定代码是否按预期工作。我制作了一个虚拟 cond_probs 文件，其中有两个数据点的“总”观察值为“1”：

total   0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
1.000   0.000   0.000   0.000   0.000   0.000   1.000   0.000   0.000   0.000   0.000
1.000   0.000   1.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000   0.000

因此，我希望它们都针对具有“1”观察值的数据点进行采样，因此得到大约 50% 的观察值的平均值为“0.2”，而 50% 的观察值的平均值为“0.6”。但是平均值始终为 0.2：

sims_test[[1]]
 [1] 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2

即使我采样了 10000 次，所有观察结果都是 0.2，而不是 0.6。我对代码的理解是，它应该从 cond_probs 中为每个观察随机选择一个具有相似大小的新行，但在这种情况下似乎没有这样做。是我误解了代码还是我的输入不正确仍然存在问题？

整个 cond_probs 文件可以在以下地址找到：

cond_probs

更新 3

在运行模拟时将 sapply 更改为 lapply 可解决此问题。

我认为保留cond_probs 并选择分布sampleSize 次数的另一个原因可能是最好的解决方案：选择分布的概率应该与其在cond_probs 中的频率有关。如果我们结合分布，选择具有total9 或10 的分布的几率将不再取决于这些总数的观察数。示例：如果有90 分布与total=10 和10 与total=9，则应该有90% 机会选择具有total=10 的分布。如果我们结合分布，选择 'total'= 9 或 10 的分布的几率不会变成 50/50（这不是理想的）吗？

【问题讨论】：

我建议您查找bootstrapping 和conditional probability。
我不清楚您希望如何解决您的问题，或者您提出的解决方案是否合理。您似乎将数据中的观察到的频率（条件或其他）与数据的概率分布混淆了......也许这是通过引导程序进行的合理方法，或者可能不是，您应该查看具有信息/非信息先验的条件后验分布....根据您的描述，我个人不清楚
感谢建设性的 cmets。你说得对，我对术语很草率，这本质上是一种引导方法。我会尝试纠正这一点。尽管我仍然对概率分布和经验分布中观察到的频率之间的差异感到困惑。如果我想根据观察到的频率分布以概率从经验分布中重新采样，那么频率分布和概率分布是否指代同一事物？
我很乐意在私聊中进一步讨论
考虑 P(Y|X= x*) ~ N(\mu,\sigma)（或一些 f(\theta)）。你从这样的分布中观察到 y_1 , ..., y_10 。来自 y_1, ... y_10 的引导样本显然与来自 N(\mu, \sigma) 或基于先验 N(\mu, \sigma) 的后验样本以及来自 y_1, ... , y_10。这就是我的观点。同样，也许引导程序适合您的问题，“频率分布和概率分布 [不] 指的是同一件事”。

标签： r simulation probability resampling replicate

【解决方案1】：

我只是写了一个函数ps，它从cond_probs中选择一个合适的分布：

N <- 10  # The sampled values are 0.1, 0.2, ... , N/10
M <- 8   # number of distributions in "cond_probs"

#-------------------------------------------------------------------
# Example data:

set.seed(1)

cond_probs <- matrix(0,M,N)

is.numeric(cond_probs)

for(i in 1:nrow(cond_probs)){ cond_probs[i,] <- dnorm((1:N)/M,i/M,0.01*N) }

is.numeric(cond_probs)

total <- sort( sample(1:80,nrow(cond_probs)) )
cond_probs <- cbind( total, cond_probs/rowSums(cond_probs) )

colnames(cond_probs) <- c( "total", paste("P",1:N,sep="") )

#---------------------------------------------------------------------
# A function that chooses an appropiate distribution from "cond_prob",
# depending on the number of observations "numObs":

ps <- function( numObs,
                similarityLimit = 0.1 )
{
  similar <- which( abs(cond_probs[,"total"] - numObs) / numObs < similarityLimit )

  if ( length(similar) == 0 )
  { 
    return(NA)
  }
  else
  {
    return( cond_probs[similar[sample(1:length(similar),1)],-1] )
  }
}

#-----------------------------------------------------------------
# A function that simulates data using a distribution that is
# appropriate to the number of observations, if possible:

simulateData <- function( numObs, sampleSize )
{
  if (any(is.na(ps(numObs))))
  {
    return (NA)
  }
  else
  {
    return( rowSums(
              replicate(
                numObs,
                replicate( sampleSize, sample((1:N)/10, 1, prob = ps(numObs))))
            ) / numObs )
  }
}

#-----------------------------------------------------------------
# Test:

sampleSize <- 30
set.seed(42)
sims <- lapply( 1:80, simulateData, sampleSize )

cond_probs中的分布：

    total           P1           P2           P3           P4           P5           P6           P7           P8           P9          P10
[1,]    16 6.654875e-01 3.046824e-01 2.923948e-02 5.881753e-04 2.480041e-06 2.191926e-09 4.060763e-13 1.576900e-17 1.283559e-22 2.189990e-28
[2,]    22 2.335299e-01 5.100762e-01 2.335299e-01 2.241119e-02 4.508188e-04 1.900877e-06 1.680045e-09 3.112453e-13 1.208647e-17 9.838095e-23
[3,]    30 2.191993e-02 2.284110e-01 4.988954e-01 2.284110e-01 2.191993e-02 4.409369e-04 1.859210e-06 1.643219e-09 3.044228e-13 1.182153e-17
[4,]    45 4.407425e-04 2.191027e-02 2.283103e-01 4.986755e-01 2.283103e-01 2.191027e-02 4.407425e-04 1.858391e-06 1.642495e-09 3.042886e-13
[5,]    49 1.858387e-06 4.407417e-04 2.191023e-02 2.283099e-01 4.986746e-01 2.283099e-01 2.191023e-02 4.407417e-04 1.858387e-06 1.642492e-09
[6,]    68 1.642492e-09 1.858387e-06 4.407417e-04 2.191023e-02 2.283099e-01 4.986746e-01 2.283099e-01 2.191023e-02 4.407417e-04 1.858387e-06
[7,]    70 3.042886e-13 1.642495e-09 1.858391e-06 4.407425e-04 2.191027e-02 2.283103e-01 4.986755e-01 2.283103e-01 2.191027e-02 4.407425e-04
[8,]    77 1.182153e-17 3.044228e-13 1.643219e-09 1.859210e-06 4.409369e-04 2.191993e-02 2.284110e-01 4.988954e-01 2.284110e-01 2.191993e-02

分布的均值：

> cond_probs[,-1] %*% (1:10)/10
          [,1]
[1,] 0.1364936
[2,] 0.2046182
[3,] 0.3001330
[4,] 0.4000007
[5,] 0.5000000
[6,] 0.6000000
[7,] 0.6999993
[8,] 0.7998670

31 次观测的模拟数据均值：

> sims[[31]]
 [1] 0.2838710 0.3000000 0.2935484 0.3193548 0.3064516 0.2903226 0.3096774 0.2741935 0.3161290 0.3193548 0.3032258 0.2967742 0.2903226 0.3032258 0.2967742
[16] 0.3129032 0.2967742 0.2806452 0.3129032 0.3032258 0.2935484 0.2935484 0.2903226 0.3096774 0.3161290 0.2741935 0.3161290 0.3193548 0.2935484 0.3032258

合适的分布是第三种：

> ps(31)
          P1           P2           P3           P4           P5           P6           P7           P8           P9          P10 
2.191993e-02 2.284110e-01 4.988954e-01 2.284110e-01 2.191993e-02 4.409369e-04 1.859210e-06 1.643219e-09 3.044228e-13 1.182153e-17

【讨论】：

非常感谢。我会试试看，让你知道。
你的答案对我来说非常适合你的 cond_probs 文件，但当我将它应用到我的实际 cond_probs 文件时却不行。我更新了答案以显示尝试时会发生什么。你知道为什么我的文件可能会产生不同的结果吗？再次感谢您的帮助，我认为这一定是一个小文件格式问题导致它无法正常工作。
cond_probs 的结构完全符合我们的要求：total 必须是第一列。我无法重现错误并通过使用问题中的cond_probs 矩阵的测试结果来增强我的答案。正如预期的那样，生成的 sims_test 是一个长度为 50 的列表。但问题中显示的 sims_test 看起来像一个 9×9 矩阵。我建议您将脚本复制并粘贴到我的答案中，看看是否得到相同的结果。如果是这样，请从文件中再次读取cond_probs，并将其替换为head(cond_probs,n)，并使用不断增加的数字n。 n=9 就是上面的例子。
感谢您的 cmets 和有用的建议。我遵循了您的建议，发现它在第 518 行停止工作。如果我从前 600 行中删除此行，它可以工作，但整个 'cond_probs' 文件仍然不起作用，所以我假设还有其他例外情况。我在问题中添加了违规行。知道是什么导致了问题吗？
我发现了更多违规行（请参阅更新的答案），但是当我将这些行放回一个较小的“cond_probs”文件时，它并没有停止工作，所以我不确定潜在的问题是什么。你知道它可能是什么吗？如果有帮助，我很乐意将整个文件（21813 行）发送给您（它太大而无法粘贴到问题中）。