【问题标题】:How to make a probability simulator in R?如何在 R 中制作概率模拟器?
【发布时间】:2017-11-21 04:28:15
【问题描述】:

使用以下数据框:

A1  A2  EFF       FRQ      
A   G   0.0125    0.4578  
T   C   0.0143    0.1293    
T   C   -0.017    0.8984  
A   G   -0.018    0.8945   
A   G   -0.009    0.8652   
A   G   0.0001    0.3931   

我想根据FRQ 列从效果大小中“抽取”两个概率。我想创建一个名为 sim_1 的新列,其中 45.78% 的时间,EFF 保持它的符号,而 54.22% 的时间,EFF 切换它的符号。然后,我想为每一行总结其中两个随机事件。例如,假设生成了两个随机数 0-100。 78.33 和 32.16。我会将任何 EFF 不变的指示。由于我随机掷出 78 和 32,总和将为 -0.0125(对于 78.33 掷骰)和 0.0125 对于(32.16)掷骰,等于 0。

在第二行,假设我滚动了两个随机数 88.22 和 67.10。因为这两个数字都不低于 12.93,所以EFF 符号将在 88.22 和 67.10 滚动中翻转,使我们得到 -0.0286 (-0.0143 + -0.0143) 的总和。

我想以这种方式做 500 个模拟列,以便最终输出如下所示:

A1  A2  EFF       FRQ      Sim_1   Sim_2   Sim_3...
A   G   0.0125    0.4578   0       -       -
T   C   0.0143    0.1293   -0.0286 -       -
T   C   -0.017    0.8984  -        -       -
A   G   -0.018    0.8945  -        -       -
A   G   -0.009    0.8652  -        -       -
A   G   0.0001    0.3931  -        -       -

注意:如果你生成一个输出文件,它可能与我的不匹配,因为它是基于随机性的。

【问题讨论】:

    标签: r simulation probability


    【解决方案1】:

    使用您的数据:

    tmp_df <- structure(list(A1 = structure(c(1L, 2L, 2L, 1L, 1L, 1L), 
                                            .Label = c("A", "T"), class = "factor"), 
                             A2 = structure(c(2L, 1L, 1L, 2L, 2L, 2L),
                                            .Label = c("C", "G"), class = "factor"), 
                             EFF = c(0.0125, 0.0143, -0.017, -0.018, -0.009, 1e-04), 
                             FRQ = c(0.4578, 0.1293, 0.8984, 0.8945, 0.8652, 0.3931)),
                        .Names = c("A1", "A2", "EFF", "FRQ"), class = "data.frame", row.names = c(NA, -6L))
    

    执行以下操作

    set.seed(0)
    
    tmp_results <- lapply(1:500, function(i) rowSums(2 * (0.5 - (matrix(runif(nrow(tmp_df) * 2), ncol = 2) >= tmp_df$FRQ)) * tmp_df$EFF))
    
    
    
    tmp_out <- as.data.frame(tmp_results)
    names(tmp_out) <- paste("Sim", 1:500)
    
    tmp_out <- cbind(tmp_df, tmp_out)
    

    制作:

    > tmp_out[, 1:10]
      A1 A2     EFF    FRQ   Sim 1   Sim 2   Sim 3   Sim 4   Sim 5   Sim 6
    1  A  G  0.0125 0.4578 -0.0250  0.0000  0.0250 -0.0250  0.0000  0.0250
    2  T  C  0.0143 0.1293 -0.0286 -0.0286 -0.0286 -0.0286  0.0000 -0.0286
    3  T  C -0.0170 0.8984 -0.0340 -0.0340 -0.0340 -0.0340 -0.0340 -0.0340
    4  A  G -0.0180 0.8945 -0.0360  0.0000 -0.0360 -0.0360 -0.0360 -0.0360
    5  A  G -0.0090 0.8652  0.0000 -0.0180 -0.0180 -0.0180 -0.0180  0.0000
    6  A  G  0.0001 0.3931  0.0002 -0.0002 -0.0002  0.0000 -0.0002  0.0000
    

    lapply 步骤说明:

    1) matrix(runif(nrow(tmp_df) * 2)
    Draw two columns filled with random numbers drawn uniformly in the interval [0, 1].
    Alternatively, you can look into using `rbinom`.
    
    2) 2 * (... >= tmp_df$FRQ) * tmp_df$EFF
    Create (-1, 1) indicator to see whether `EFF` should be fliped, then multiply, exploiting conformability rules.
    
    3) lapply(...) 
    Do the above 500 times.
    

    其余的只是标记,并将模拟结果绑定到您的原始数据。

    【讨论】:

    • tmp_results &lt;- lapply(1:500, function(i) rowSums(2 * (0.5 - matrix(runif(nrow(tmp_df) * 2), ncol = 2) &gt;= tmp_df$FRQ) * tmp_df$EFF)) (tmp_results) 行正在向我的 R 输出数百万行并导致我的程序崩溃
    • 对此感到抱歉,不确定(tmp_results) 是如何出现的。只需删除它,因为这是打印所有结果。
    • 第二行有一些模拟,我得到 0.05758(基本上是效果列的 4 倍)。你知道为什么会这样吗?模拟的大小不应超过效果大小值的两倍。
    • 对不起,我在尝试使指示器功能时将括号放在错误的位置。现在修好了,我希望。
    猜你喜欢
    • 2019-07-26
    • 2021-11-18
    • 2013-11-10
    • 1970-01-01
    • 1970-01-01
    • 2018-08-29
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多