使用多个概率分布在 R 中模拟数据
Simulating data in R with multiple probability distributions
我正在尝试通过自举来模拟数据，以使用漏斗图为我的真实数据创建置信带。我正在建立已接受答案 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) 时，它随机抽取 10 个值（可以是 0.1、0.2、0.3、0.4、0.5、0.6 或 0.7）基于在概率分布 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。
现在我不想对所有模拟都使用一种频率分布，而是根据每个数据点下的观察数量有条件地使用不同的频率分布。
我制作了一个 table、cond_probs，每个真实数据点都有一行。一列包含 total 个观察值，一列给出每个观察值的频率。
cond_probs的例子table:
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 table 并随机 select 一行，其中 total 个观察值的大小与 k 相似：0.9k-1.1k。如果不存在这样的行，则继续。
一旦数据点被 selected，使用 cond_probs 中该行的概率分布就像在原始代码中使用 ps 一样，随机抽样 k 观察次数并输出这些观察结果的平均值。
对于 replicate 的每个 n 迭代，从 cond_probs 的所有行中随机抽样替换一个新数据点 total类似于 k ( 0.9k-1.1k) 的当前值。
这个想法是，对于这个数据集，应该根据数据点下的观察数量来确定要使用的概率分布。这是因为在此数据集中，观察的概率受观察次数的影响（由于遗传连锁和背景 selection，具有更多 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
我没有注意到他们有什么异常。他们都有一个 'total' 的 9 个站点。如果我删除这些行并且 运行 'cond_probs' 文件只包含这些行之前的行，那么代码就可以工作。但是肯定还有其他有问题的行，因为整个 'cond_probs' 仍然不起作用。
我试着将这些有问题的行放回一个较小的 'cond_probs' 文件中，然后这个文件就可以工作了，所以我很困惑，因为看起来这些行本身并没有问题。另一方面，他们总共有 9 个站点的事实表明存在某种因果关系。
如果有帮助，我很乐意私下共享整个文件ps，因为我不知道下一步该怎么做才能进行故障排除。
出现的另一个问题是我不确定代码是否按预期工作。我制作了一个虚拟 cond_probs 文件，其中有两个数据点 'total' 为“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 中随机 select 插入一个新行，每个观察的大小相似，但在这种情况下似乎没有这样做。是我误解了代码还是我的输入不正确仍然是个问题？
整个cond_probs文件可以在以下地址找到：
更新 3
当 运行 模拟解决了这个问题时，将 sapply 更改为 lapply。
我认为保持 cond_probs 不变并选择分布 sampleSize 次的另一个原因可能是最好的解决方案：选择分布的概率应该与其在 cond_probs。如果我们组合分布，选择具有 total 9 或 10 的分布的几率将不再取决于具有这些总数的观察次数。示例：如果有 90 分布 total=10 和 10 分布 total=9 应该有 90% 机会选择 total=10 分布。如果我们组合分布，选择 'total'= 9 或 10（这不理想）的分布的几率不会变成 50/50 吗？
我只是写了一个函数 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
我正在尝试通过自举来模拟数据，以使用漏斗图为我的真实数据创建置信带。我正在建立已接受答案 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) 时，它随机抽取 10 个值（可以是 0.1、0.2、0.3、0.4、0.5、0.6 或 0.7）基于在概率分布 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。
现在我不想对所有模拟都使用一种频率分布，而是根据每个数据点下的观察数量有条件地使用不同的频率分布。
我制作了一个 table、cond_probs，每个真实数据点都有一行。一列包含 total 个观察值，一列给出每个观察值的频率。
cond_probs的例子table:
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 table 并随机 select 一行，其中 total 个观察值的大小与 k 相似：0.9k-1.1k。如果不存在这样的行，则继续。
一旦数据点被 selected，使用 cond_probs 中该行的概率分布就像在原始代码中使用 ps 一样，随机抽样 k 观察次数并输出这些观察结果的平均值。
对于 replicate 的每个 n 迭代，从 cond_probs 的所有行中随机抽样替换一个新数据点 total类似于 k ( 0.9k-1.1k) 的当前值。
这个想法是，对于这个数据集，应该根据数据点下的观察数量来确定要使用的概率分布。这是因为在此数据集中，观察的概率受观察次数的影响（由于遗传连锁和背景 selection，具有更多 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
我没有注意到他们有什么异常。他们都有一个 'total' 的 9 个站点。如果我删除这些行并且 运行 'cond_probs' 文件只包含这些行之前的行，那么代码就可以工作。但是肯定还有其他有问题的行，因为整个 'cond_probs' 仍然不起作用。
我试着将这些有问题的行放回一个较小的 'cond_probs' 文件中，然后这个文件就可以工作了，所以我很困惑，因为看起来这些行本身并没有问题。另一方面，他们总共有 9 个站点的事实表明存在某种因果关系。
如果有帮助，我很乐意私下共享整个文件ps，因为我不知道下一步该怎么做才能进行故障排除。
出现的另一个问题是我不确定代码是否按预期工作。我制作了一个虚拟 cond_probs 文件，其中有两个数据点 'total' 为“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 中随机 select 插入一个新行，每个观察的大小相似，但在这种情况下似乎没有这样做。是我误解了代码还是我的输入不正确仍然是个问题？
整个cond_probs文件可以在以下地址找到：
更新 3
当 运行 模拟解决了这个问题时，将 sapply 更改为 lapply。
我认为保持 cond_probs 不变并选择分布 sampleSize 次的另一个原因可能是最好的解决方案：选择分布的概率应该与其在 cond_probs。如果我们组合分布，选择具有 total 9 或 10 的分布的几率将不再取决于具有这些总数的观察次数。示例：如果有 90 分布 total=10 和 10 分布 total=9 应该有 90% 机会选择 total=10 分布。如果我们组合分布，选择 'total'= 9 或 10（这不理想）的分布的几率不会变成 50/50 吗？
我只是写了一个函数 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