如何生成相关的二元变量

Question

我需要生成一系列具有给定相关函数的 N 个随机二元变量。令 x = {x_i} 为一系列二进制变量（取值 0 或 1，i 从 1 到 N）。边际概率给定 Pr(x_i = 1) = p，变量应该在方式如下：

更正[x_ix_j ] = const × |i-j|^-α（对于 i!=j）

其中 α 是一个正数。

如果更简单，考虑相关函数：

更正[x_ix_j ] = (|i-j|+1)^-α

重要的部分是我想研究相关函数像幂律一样时的行为。 （不是α^|i-j|）

是否可以生成这样的系列，最好是在 Python 中？

Answer 1

感谢您的所有意见。我在 Chul Gyu Park 等人的可爱小文章中找到了我的问题的答案，所以如果有人遇到同样的问题，请查找：

“一种生成相关二进制变量的简单方法”（jstor.org.stable/2684925）

对于一个简单的算法。如果相关矩阵中的所有元素都是正数，并且对于一般边际分布 Pr(x_i)=p_i，则该算法有效。

j

Answer 2

您正在描述一个随机的过程，这对我来说似乎很难...如果您消除了二进制 (0,1) 要求，而是指定了期望值和方差，可以将其描述为通过 1 极点低通滤波器馈送的白噪声发生器，我认为这会给您 α^|ij| 特性。

这实际上可能符合 mathoverflow.net 的标准，具体取决于它的措辞方式。让我试试问....

---

更新：我为 α^|i-j| 案例做了ask on mathoverflow.net。但也许有一些想法可以适应你的情况。

Answer 3

在RSeek 的快速搜索显示R 有包

• bindata
• binarySimCLF

这样做。

Answer 4

蛮力解决方案是将问题的约束表示为具有2^N 变量pr(w) 的线性程序，其中w 的范围是所有长度为N 的二进制字符串。首先pr是一个概率分布的约束：

for all w: 0 <= pr(w) <= 1
sum_w pr(w) = 1

二、每个变量的期望为p的约束：

for all i: sum_{w such that w[i] = 1} pr(w) = p

三、协方差约束：

for all i < j: sum_{w such that w[i] = w[j] = 1} pr(w) = const * |j - i|^alpha - p^2

这很慢，但粗略的文献搜索没有更好的结果。如果你决定实现它，这里有一些带有 Python 绑定的 LP 求解器：http://wiki.python.org/moin/NumericAndScientific/Libraries

Answer 5

将分布x_i表示为一些独立基分布f_j的线性组合：x_i = a_i1f₁ + a_i2f₂ + 。 .. .让我们将 f_j 约束为均匀分布在 0..1 或 {0,1}（离散）中的自变量。现在让我们用矩阵形式表达我们所知道的一切：

Let X be the vector (x1, x2, .., xn)
Let A be the matrix (a_ij) of dimension (k,n) (n rows, k columns)
Let F be the vector (f1, f2, .., fk) 
Let P be the vector (p1, p2, .., pn)
Let R be the matrix (E[x_i,x_j]) for i,j=1..n
Definition of the X distribution: X = A * F
Constraint on the mean of individual X variables: P = A * (1 ..k times.. 1)
Correlation constraint: AT*A = 3R or 2R in the discrete case (because E[x_i x_j] = 
  E[(a_i1*f_1 + a_i2*f_2 + ...)*(a_j1*f_1 + a_j2*f_2 + ...)] =
  E[sum over p,q: a_ip*f_p*a_jq*f_q] = (since for p/=q holds E[f_p*f_q]=0)
  E[sum over p: a_ip*a_jp*f_p^2] =
  sum over p: a_ip*a_jp*E[f_p^2] = (since E[f_p^2] = 1/3 or 1/2 for the discrete case)
  sum over p: 1/3 or 1/2*a_ip*a_jp
And the vector consisting of those sums over p: a_ip*a_jp is precisely AT*A.

现在你需要解两个方程：

AT*A      = 3R (or 2R in the discrete case)
A*(1...1) = P

第一个方程的解对应于求矩阵 3R 或 2R 的平方根。例如见http://en.wikipedia.org/wiki/Cholesky_factorization 和一般http://en.wikipedia.org/wiki/Square_root_of_a_matrix。 第二个也应该做点什么:)

我请周围的数学家纠正我，因为我很可能将 ATA 与 AAT 混为一谈，或者做的更糟。

要生成 x_i 的值作为基分布的线性混合，请使用两步过程：1）使用均匀随机变量来选择一个的基础分布，用相应的概率加权，2) 使用选择的基础分布生成结果。

Answer 6

这是一种似乎可行的直观/实验方法。

如果 b 是二进制 r.v.， m 是二进制 r.v. 的平均值， c 是你想要的相关性， rand() 生成一个 U(0,1) r.v.，并且 d 是相关的二进制 r.v。你想要的：

d = if(rand()

也就是说，如果一个统一的房车。小于所需的相关性，d = b。否则 d = 另一个随机二进制数。

我为 2000 个二进制 r.v.s 的列运行了 1000 次。 m=.5 和 c = .4 和 c = .5 相关平均值完全符合规定，分布似乎是正态的。 对于 0.4 的相关性，相关性的标准偏差为 0.02。

抱歉 - 我无法证明这始终有效，但您必须承认，这确实很容易。