如何生成对非对角元素有约束的伪随机正定矩阵？

Question

用户希望对 var/covar 矩阵中每对变量之间的相关性施加一个唯一的、非平凡的上限/下限。

例如：我想要一个方差矩阵，其中所有变量都有 0.9 > |rho(x_i,x_j)| > 0.6，rho(x_i,x_j) 为变量 x_i 和 x_j 之间的相关性。

谢谢。

---

好的，已经找到了一种快速而肮脏的解决方案，但如果有人知道更准确的方法可以到达那里，我们将不胜感激。

---

我失去了原来的登录名，所以我在新的登录名下重新发布了这个问题。 The previous iteration得到以下答案

*你的意思是伪随机，这是 semi 随机的正确术语 – Robert Gould

*好点，但我认为他的意思是半伪随机（在谈论计算机随机性时假定伪：-p） - fortran

*“相关”是指“协方差”吗？ ——斯万特

*不，我的意思是相关性。我想生成一个正定矩阵，使得所有相关性都比平凡的界限更紧密。 – 瓦克

*看我的回答。您是否坚持样本相关性位于指定范围内，或者只是生成样本的总体相关性？如果您的问题是前者，我确实提出了一个可行的想法。 – 木片

*woodship：不，我担心你的解决方案不起作用，请在原始威胁中查看我的答案（上面的链接）。谢谢。

Answer 1

这是您在原帖中对我的回答的回复：

“来人，一定有更简单的”

对不起，没有。想中彩票是不够的。要求小熊队赢得系列赛是不够的。你也不能只是要求一个数学问题的解决方案，然后突然发现它很容易。

使用指定范围内的样本参数生成伪随机偏差的问题并非微不足道，至少如果偏差在任何意义上都是真正的伪随机。根据范围，一个人可能是幸运的。我提出了一个拒绝方案，但也表示这不太可能是一个好的解决方案。如果相关性有很多维度和狭窄的范围，那么成功的概率就很低。样本量也很重要，因为这将推动所得相关性的样本方差。

如果您真的想要一个解决方案，您需要坐下来明确、准确地指定您的目标。您是否想要一个具有名义上指定相关结构但对相关性有严格限制的随机样本？满足目标界限的样本相关矩阵是否令人满意？是否也给出了方差？

Answer 2

您可以创建一组大小为 M 和单位方差的 N 个随机向量。并向它们添加一个随机向量（大小 N 和单位方差）乘以某个数字 k。 然后你取所有这些向量之间的相关性，这将是一个正定矩阵。如果 M 非常大，则相关分布将没有方差，相关性将为：k^2/(1+k^2)。 M 越小，非对角线元素的分布就越广。 或者，您可以让 M 非常大，并将“公共向量”分别乘以不同的 k。如果您正确使用这些参数，您可能会获得更严格的控制。这里有一些 Matlab 代码来做到这一点：

clear all;
vecLarg=10;
theDim=1000;
corrDist=0*randn(theDim,1);
Baux=randn(vecLarg,theDim)+  (corrDist*randn(1,vecLarg))'+(k*ones(theDim,1)*randn(1,vecLarg))'  ;
A=corrcoef(Baux);
hist(A(:),100);

Answer 3

也许这个答案将有助于实现它：

具有这种非负确定性的一类矩阵是Wishart Distribution。并且来自 ~W() 的样本使得所有非对角线条目都在某些范围 [l,u] 之间将适合您的问题。但是，我不认为这与 [l,u] 中所有非对角线的正定矩阵的分布相同。

在维基百科页面上，有一个从 ~W() 计算的算法。

一个更简单的hackish解决方案（可能近似于此）是：

（假设 u>l 和 l>0）

1. 从多元正态图中提取，其中 Sigma = mean(l,u)。
2. 然后取样本，计算其相关矩阵 => C
3. 这个矩阵会有一些随机性（模糊），但它有多少模糊的数学有点超出我的计算范围。此 C 矩阵中的非对角线的值以 [-1,1] 为界，均值为 mean(l,u)。通过眼球，我猜测某种贝塔/指数。在任何情况下，除非 (l,u) = [-1,1]，否则 C 中关闭诊断的连续分布保证它不会表现并位于边界 (l,u) 内。
4. 您可以通过增加/减少步骤 1 中样本的长度来调整“模糊”量。我敢打赌（未经证实）C 的奇数对数的方差量与平方根成正比样本数。

因此，要真正回答这似乎并非易事！

正如其他海报所建议的那样，您可以从 Wishart 生成，然后保留您想要的属性为 true 的那些，但您可能会采样很长时间！如果你排除那些 0-确定的（那是一个词吗？），那么这应该可以很好地生成好的矩阵。然而，这并不是所有 pos-def 矩阵的真实分布，其 off-diags 在 [l,u] 中。

上面提出的哑采样方案的代码（R 中）

sigma1 <- function(n,sigma) {
    out <- matrix(sigma,n,n)
    diag(out) <- 1
    return (out)
}

library(mvtnorm)
sample_around_sigma <- function(size, upper,lower, tight=500) {
    #  size:  size of matrix
    #  upper, lower:  bounds on the corr, should be > 0
    #  tight:  number of samples to use.  ideally this
    #     would be calcuated such that the odd-diags will
    #     be "pretty likely" to fall in [lower,upper]
    sigma <- sigma1(size,mean(c(upper,lower)))
    means <- 0*1:size
    samples <- rmvnorm(n=tight, mean=means,sigma=sigma)
    return (cor(samples))
}

> A <- sample_around_sigma(5, .3,.5)
> A
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 1.0000000 0.3806354 0.3878336 0.3926565 0.4080125
[2,] 0.3806354 1.0000000 0.4028188 0.4366342 0.3801593
[3,] 0.3878336 0.4028188 1.0000000 0.4085453 0.3814716
[4,] 0.3926565 0.4366342 0.4085453 1.0000000 0.3677547
[5,] 0.4080125 0.3801593 0.3814716 0.3677547 1.0000000
> 
> summary(A[lower.tri(A)]); var(A[lower.tri(A)])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3678  0.3808  0.3902  0.3947  0.4067  0.4366 
[1] 0.0003949876

Answer 4

好的，出色的 Gregg：我们正在取得进展。将您的想法与木片的想法结合起来，就产生了这种替代方法。它在数学上非常很脏，但它似乎有效：

library(MCMCpack)
library(MASS)
p<-10
lb<-.6
ub<-.8
zupa<-function(theta){
    ac<-matrix(theta,p,p)
    fe<-rwish(100*p,ac%*%t(ac))
    det(fe)
}
ba<-optim(runif(p^2,-10,-5),zupa,control=list(maxit=10))
ac<-matrix(ba$par,p,p)
fe<-rwish(100*p,ac%*%t(ac))
me<-mvrnorm(p+1,rep(0,p),fe)
A<-cor(me)
bofi<-sqrt(diag(var(me)))%*%t(sqrt((diag(var(me)))))
va<-A[lower.tri(A)]
l1=100
while(l1>0){
    r1<-which(va>ub)
    l1<-length(r1)
    va[r1]<-va[r1]*.9
}
A[lower.tri(A)]<-va
A[upper.tri(A)]<-va
vari<-bofi*A
mk<-mvrnorm(10*p,rep(0,p),vari)
pc<-sign(runif(p,-1,1))
mf<-sweep(mk,2,pc,"*")
B<-cor(mf)
summary(abs(B[lower.tri(B)]))

基本上，这就是这个想法（比如上限 =.8 和下限 =.6），它有足够好的接受率，虽然不是 100%，但它会在项目的这个阶段完成.