1
$\begingroup$

Hi,

I have a set of N objects randomly distributed in a 2D physical space. Each object (i) generates a bernoulli random number (0 or 1) based on a marginal probability Pr(xi = 1) = p. These objects a correlated by physical distance. The closer the objects are, the larger their correlation is.

E.g. If objects i and j are co-located, they are expected to generate correlated results. For Example, if P(Xi=1)= 0.6 and P(Xj=1)=0.3 they would produce something like:

Xi= 0 1 0 1 1 1 0 1 0 1

Xj= 0 1 0 0 0 1 0 1 0 0

Such that Pr(Xi|Xj)=1

On the other hand if i and j are distant they would produce uncorrelated results such that Pr(Xi|Xj)=Pr(Xi)

I have tried to use some of the packages in Matlab (Sampling from multivariate correlated binary and poisson random variables) and R (bindata) but I could not produce an acceptable correlation matrix.

Any ideas how I can produce an acceptable correlation matrix?

BTW, I have checked the following earlier posts discrete stochastic process: exponentially correlated Bernoulli?

and

Constructing Bernoulli random variables with prescribed correlation

But I am not sure how I can relate to them.

Thanks

$\endgroup$

1 Answer 1

2
$\begingroup$

Here's a suggestion:

Define a non-negative decreasing function $w(r)$ measuring interaction strength. Given each object its own independent $N(0,1)$ random variable $N_i$. Now set $$ Y_i=\frac{\sum_{j}w(\|x_i-x_j\|)N_j}{\sqrt{\sum_j w(\|x_i-x_j\|)^2}}, $$ where $x_i$ denotes the location of the $i$th object.

Then the $Y_i$ are correlated $N(0,1)$ random variables. If two objects are co-located the normal random variables agree.

Finally set $t_i=\Phi^{-1}(p_i)$ (i.e. $\mathbb P(N < t_i)=p_i$) and set $X_i=1$ if $Y_i < p_i$ and 0 otherwise.

With this setup you can write down the covariance of $Y_i$ and $Y_k$ explicitly: it's just $$ \text{Cov}(Y_i,Y_k)=\frac{\sum_j w(\|x_i-x_j\|)w(\|x_k-x_j\|)} {\sqrt{\sum_j w(\|x_i-x_j\|)^2\sum_j w(\|x_k-x_j\|)^2}}. $$

If you write this as $\cos\theta_{ik}$ then you can write the covariance of $X_i$ and $X_k$ as an integral: $$ 1/(2\pi)\int_{ x < t_1\;,\; cos\theta_{ik}x+\sin\theta_{ik}y < t_2} e^{-(x^2+y^2)/2}\,dxdy-p_ip_k. $$

$\endgroup$
6
  • $\begingroup$ If i understood you correctly, what you are saying is to simply use the Cov(Yi,Yk) above as the covariance matrix and plug it in the bernoulli generator. If so, then Unfortunately, I tried that but I keep getting an Unacceptable correlation matrix, which seams to be not positive definite. $\endgroup$
    – alandalusi
    Apr 26, 2011 at 21:44
  • $\begingroup$ No you didn't understand me correctly. Here is my concrete suggestion. (1) Compute $t_i=\Phi^{-1}(p_i)$; (2) Compute the matrix $Cov(Y_i,Y_k)$ as above (3) Use a multivariate normal generator to build some $Y_i$'s. (4) Set $X_i=1$ if $Y_i<t_i$ and 0 otherwise. You could compute the covariance of the Bernoulli's (using the nasty integral formula I wrote down), but if your purpose is to just generate the random #s, then the procedure I described will work just as well. $\endgroup$ Apr 26, 2011 at 22:59
  • $\begingroup$ Acually, It is positive definite, but does not have a dichotomized Gaussian distribution for the correlation matrix. $\endgroup$
    – alandalusi
    Apr 26, 2011 at 23:03
  • $\begingroup$ OK, I think I understand your suggestion now. I will test it right away. $\endgroup$
    – alandalusi
    Apr 26, 2011 at 23:32
  • $\begingroup$ It works!, but I noticed that the Correlation (or Conditional Probabilities) $Pr(X{_i}=1|X{_k}=1)$ between any fixed pair (i and k) slightly change as we add more objects into the physical space. So, I changed $Cov(Y{_i},Y{_k})$ to be $w(||x{_i}-x{_k}||)$ and it worked with fixed conditional probabilities, even if I add more objects. I am still testing it, but does this make sense to you? For your information, my distance functions is $e^{(-\frac{eucaliandistance}{dcorr})}$ $\endgroup$
    – alandalusi
    Apr 27, 2011 at 16:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.