2

我需要用 R 生成一个联合概率矩阵（80*80），这样这个矩阵的均值和标准差是预先固定的。是一个对称矩阵。主对角线元素被解释为二进制变量 A i 等于 1 的概率 p ( A i )。非对角线元素是 A i 和 A j 均为 1 的概率 p ( A i A j )。这个矩阵应该响应以下条件：

0 ≤ p A i ≤ 1

max ( 0 , p A i + p A j − 1 ) ≤ p A i A j ≤ min ( p A i , p A j ) , i ≠ j

p A i + p A j + p A k − p A i A j − p A i A k − p A j A k ≤ 1 , i ≠ j , i ≠ k , j ≠ k

使用 check.commonprob 检查这些条件。

这是我尝试过的

# First I need another function to make the matrix symmetric 


makeSymm <- function(m) {
    m[upper.tri(m)] <- t(m)[upper.tri(m)]
    return(m)
}



install.packages("bindata")
library (bindata)

# Function to create the joint probabilites matrix:

joint <- function(b, x, y, k) {
    repeat {
        b[lower.tri(b, diag=TRUE)] <- rnorm(k*(k+1)/2, mean=x,sd=y)
        b <- makeSymm(b)
        if (check.commonprob(b)==TRUE) break
    }
    return(b)
}

c <- matrix(0, 4, 4)
c <- joint(c, 0.5, 0.25, 4)
c
              [,1]      [,2]      [,3]      [,4]
    [1,] 0.3956524 0.3156023 0.3007369 0.2837493
    [2,] 0.3156023 0.4330955 0.2920143 0.4145124
    [3,] 0.3007369 0.2920143 0.7921709 0.3616905
    [4,] 0.2837493 0.4145124 0.3616905 0.4318048

问题是这个函数适用于非常小的矩阵（4 * 4）但是如果我取（6 * 6），运行时间太长所以我停止它。我需要将此函数推广到 80*80 矩阵。有什么建议么 ？

4

0 回答 0