生成具有总和和极大值约束的随机矩阵

机器算法验证 可能性 随机生成 随机矩阵
2022-03-14 02:34:44

我想生成一个随机方阵,使行归一化为一,对角线元素是其列的最大值。如果有一种有效的方法来统一采样这些矩阵?

2×2矩阵很容易创建。第一列是通过从中采样两个独立的制服并将最大值移动到第一个条目而生成的。然后第二列是第一列的补码。我在更高维度上遵循类似过程的尝试感觉非常特别,最麻烦的是,最终补充列的分布与其他列的分布非常不同。某种类型的拒绝抽样是我最好的选择吗?[0,1]

动机是生成似然矩阵,行是世界的状态,列是观察到的信号。每个信号最有可能在其相应的状态下出现,但不一定是该状态下最有可能出现的信号。

编辑:确切的分布并不重要,除了看起来很自然。不过,我想避免对角线对于行也是最大的。

这是我自己以前的 ah-hoc 尝试。

gen.likelihood.matrix <- function(T){
    # T: Number of states

    mat <- matrix(runif(T^2), nrow=T) # Intialize with uniform variates
    diag(mat) <- 1 # To guarantee diag is max in column

    candidate.mat <- matrix(nrow=T, ncol=T)
    # Loop to check for constraint satisfaction
    for(k in 1:100){
        weights <- runif(T-1)
        weights <- weights / sum(weights)

        for(i in 1:T){
            # Rescale columns
            candidate.mat[i,-T] <- mat[i,-T] * weights
            # Replace last column with complement to meet row constraint
            candidate.mat[i,T] <- 1 - sum(candidate.mat[i,-T])
        }

        # Does the last column meet the maximality constraint?
        if(candidate.mat[T,T] > max(candidate.mat[-T,T])){
            # Shuffle to minimize the distortion from last column
            reorder <- sample(1:T,T)
            return(candidate.mat[reorder,reorder])
        }
    }

    # Recurse to get new initial values if necessary
    return(genLikeliMatrix(T))
}

这可以有效地为我提供矩阵,但由于我处理最后一列的方式,我在对角线上得到了一个奇怪的双峰分布。 对角线元素的分布

4个回答

好的,为了推动这些努力(但有些犹豫),我提供了这种方法:首先生成对角线元素。使它们成为大常数根据您想要的任何(非负)分布生成所有非对角元素 iid。规范化行。检查列最大条件。如果违反,请重复。

通过使初始常数足够大,可以使预期的重复次数变小。

显然,对角线元素是独立同分布的,非对角线元素是独立同分布的,但是(当然)这两种分布不同。

这是一些可以玩的代码。

n <- 8                                     # Matrix dimension
y <- rep(1 + 3/sqrt(n),n)                  # A large constant compared to entries in x
x <- matrix(runif(n^2), ncol=n)            # Here, uniform distributions off diagonal
x[cbind(1:n,1:n)] <- y                     # (Paste in the diagonal)
z <- t(apply(x, 1, function(u) u/sum(u)))  # Normalize the rows
which(1:n != apply(z, 2, which.max))       # Find all columns violating the conditions

(希望integer(0)作为输出;否则,将输出最大值不是对角线的列的索引。)我尝试了n从 3 到 300 的范围。

绘制列很有启发性:

plot(z[1,], type="n")
apply(z, 2, function(u) lines(u, col=(256*runif(1))))

首先,我将生成一个矩阵,使得以指数分布,所有 ,否则MMijμij0

然后,计算你的矩阵

Xij=eMijjeMij

为了有效地从分类分布(总概率为 1 的列表)中采样,只需使用每个结果的概率来构建一个霍夫曼树,跟踪每个节点左侧后代的总概率。从均匀分布中采样,并找到其 cdf 最小的叶节点,该叶节点大于您的采样值。(这具有分类分布熵的摊销成本。)

首先生成随机行,然后将每一行除以其总和,以便它们总和为 1。然后以特定方式对行进行排序,在第一行中,您希望在其第一个位置(即第 1 列的最大值)中包含元素的行。因此,在第 n 行中,您想要具有第 n 列最大值的行.

只需从Dirichlet 分布生成每一行,然后在需要时重新排序以将最大值放在对角线上(或选择参数以使最大值可能在对角线上)。