如何在 R 中找到一个永久的方阵(对于一般维度 nxn)?特别是,我正在尝试为独立但非相同分布的人口找到一个顺序统计的pdf,其中包括计算一个矩阵的永久值,该矩阵的元素是原始人口的pdfs和cdfs
谢谢
tl; dr这是一个非平凡的计算问题,似乎没有在 R 中实现,并且计算难度很大,可能需要编译解决方案。最好的办法是编写包装这个开源 C 实现的 R 代码。
根据相关的 Wikipedia 文章,“Ryser”看起来是一个很好的搜索词,可用于查找此计算的实现。library("sos"); findFn("Ryser")
仅找到Spearman 等级相关性的帮助,它说
Spearman 等级相关统计量的精确零分布的计算在 n 中呈指数级困难。这个包使用预先计算的 n <= 22 的精确分布,使用 Ryser 的公式获得,应用于适当的单项式永久。
这甚至不是一般的实现,而是一种特殊情况。 谷歌搜索“永久 Ryser”并没有提出任何实现,直到我们深入到这里,这是 MATLAB 代码。谷歌搜索“永久 Ryser 实现”出现了这个 CodeProject 页面,它提供了相当简单的 C 代码,在相当宽松的Code Project Open License下获得许可。
那么永久的matrix(1:9, 3)
就是:
install.packages("permute"); library(permute)
A<-matrix(1:9, 3)
# Error: sum( apply( allPerms(1:3), 1, function(r) prod( A[1:3, r]) ) )
该allPerms
函数似乎遗漏了原始向量,因此需要 Ben Bolker 的更正之一,我应该用它cbind
来构造 的项目的索引A
:
sum( apply( rbind(1:3,allPerms(1:3)), 1,
function(r) prod( A[cbind(1:3, r)]) ) )
这些值都是正数并且没有减法这一事实表明了为什么不建议使用这种“幼稚”的定义实现的原因。
A <- matrix(1:16,4)
sum( apply( rbind(1:4,allPerms(1:4)), 1,
function(r) prod( A[cbind(1:4, r)]) ) )
#[1] 55456