Create a permutation with same autocorrelation
<p>My question is similar to <a href="https://stackoverflow.com/questions/33898665/python-generate-array-of-specific-autocorrelation">this one</a>, but with the difference that I need an array of zeros and ones as output. I have an original time series of zeroes and ones with high autocorrelation (i.e., the ones are clustered). For some significance-testing I need to create random arrays with the same number of zeroes and ones. I.e. permutations of the original array, however, also the autocorrelation should stay the same/similar to the original so a simple <code>np.permutation</code> does not help me. </p>

<p>Since I'm doing multiple realizations I would need a solution which is as fast as possible. Any help is much appreciated.   </p>

<p>According to the question to which you refer, you would like to permute <code>x</code> such that </p>

<pre><code>np.corrcoef(x[0: len(x) - 1], x[1: ])[0][1]
</code></pre>

<p>doesn't change. </p>

<p>Say the sequence <em>x</em> is composed of </p>

<p><em>z<sub>1</sub> o<sub>1</sub> z<sub>2</sub>  o<sub>2</sub>  z<sub>3</sub>  o<sub>3</sub> ... z<sub>k</sub> o<sub>k</sub></em>, </p>

<p>where each <em>z<sub>i</sub></em> is a sequence of 0s, and each <em>o<sub>i</sub></em> is a sequence of 1s. (There are four cases, depending on whether the sequence starts with 0s or 1s, and whether it ends with 0s or 1s, but they're all the same in principle).</p>

<p>Suppose <em>p</em> and <em>q</em> are each permutations of <em>{1, ..., k}</em>, and consider the sequence </p>

<p><em>z<sub>p[1]</sub> o<sub>q[1]</sub> z<sub>p[2]</sub>  o<sub>q[2]</sub>  z<sub>p[3]</sub>  o<sub>q[3]</sub> ... z<sub>p[k]</sub> o<sub>q[k]</sub></em>,</p>

<p>that is, each of the run-length sub-sequences of 0s and 1s have been permuted internally. </p>

<p>For example, suppose the original sequence is</p>

<p><em>0, 0, 0, 1, 1, 0, 1</em>.</p>

<p>Then </p>

<p><em>0, 0, 0, 1, 0, 1, 1</em>,  </p>

<p>is such a permutation, as well as</p>

<p><em>0, 1, 1, 0, 0, 0, 1</em>,</p>

<p>and</p>

<p><em>0, 1, 0, 0, 0, 1, 1</em>.</p>

<p>Performing this permutation will not change the correlation:</p>

<ul>
<li>within each run, the differences are the same</li>
<li>the boundaries between the runs are the same as before</li>
</ul>

<p>Therefore, this gives a way to generate permutations which do not affect the correlation. (Also, see at the end another far simpler and more efficient way which can work in many common cases.)</p>

<p>We start with the function <code>preprocess</code>, which takes the sequence, and returns a tuple <code>starts_with_zero, zeros, ones</code>, indicating, respectively,</p>

<ul>
<li>whether <code>x</code> began with 0</li>
<li>The 0 runs</li>
<li>The 1 runs</li>
</ul>

<p>In code, this is</p>

<pre><code>import numpy as np
import itertools

def preprocess(x):
    def find_runs(x, val):
        matches = np.concatenate(([0], np.equal(x, val).view(np.int8), [0]))
        absdiff = np.abs(np.diff(matches))
        ranges = np.where(absdiff == 1)[0].reshape(-1, 2)
        return ranges[:, 1] - ranges[:, 0]

    starts_with_zero = x[0] == 0

    run_lengths_0 = find_runs(x, 0)
    run_lengths_1 = find_runs(x, 1)
    zeros = [np.zeros(l) for l in run_lengths_0]
    ones = [np.ones(l) for l in run_lengths_1]

    return starts_with_zero, zeros, ones
</code></pre>

<p>(This function borrows from an answer to <a href="http://stackoverflow.com/questions/1066758/find-length-of-sequences-of-identical-values-in-a-numpy-array">this question</a>.)</p>

<p>To use this function, you could do, e.g., </p>

<pre><code>x = (np.random.uniform(size=10000) &gt; 0.2).astype(int)

starts_with_zero, zeros, ones = preprocess(x)
</code></pre>

<p>Now we write a function to permute internally the 0 and 1 runs, and concatenate the results:</p>

<pre><code>def get_next_permutation(starts_with_zero, zeros, ones):
    np.random.shuffle(zeros)
    np.random.shuffle(ones)

    if starts_with_zero:
        all_ = itertools.izip_longest(zeros, ones, fillvalue=np.array([]))
    else:
        all_ = itertools.izip_longest(ones, zeros, fillvalue=np.array([]))
    all_ = [e for p in all_ for e in p]

    x_tag = np.concatenate(all_)

    return x_tag
</code></pre>

<p>To generate another permutation (with same correlation), you would use</p>

<pre><code>x_tag = get_next_permutation(starts_with_zero, zeros, ones)
</code></pre>

<p>To generate many permutations, you could do:</p>

<pre><code>starts_with_zero, zeros, ones = preprocess(x)

for i in range(&lt;number of permutations needed):
    x_tag = get_next_permutation(starts_with_zero, zeros, ones)
</code></pre>

<hr>

<p><strong>Example</strong></p>

<p>Suppose we run</p>

<pre><code>x = (np.random.uniform(size=10000) &gt; 0.2).astype(int)
print np.corrcoef(x[0: len(x) - 1], x[1: ])[0][1]

starts_with_zero, zeros, ones = preprocess(x)

for i in range(10):
    x_tag = get_next_permutation(starts_with_zero, zeros, ones)

    print x_tag[: 50]
    print np.corrcoef(x_tag[0: len(x_tag) - 1], x_tag[1: ])[0][1]
</code></pre>

<p>Then we get:</p>

<pre><code>0.00674330566615
[ 1.  1.  1.  1.  1.  0.  0.  1.  1.  1.  1.  1.  1.  1.  1.  0.  1.  0.
  1.  1.  0.  1.  1.  1.  1.  0.  1.  1.  0.  0.  1.  0.  1.  1.  1.  1.
  0.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
0.00674330566615
[ 1.  0.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  0.  1.  1.  0.  1.  1.  1.  1.  1.  1.  0.  0.  1.  0.
  1.  1.  1.  1.  0.  0.  0.  1.  1.  1.  1.  1.  1.  1.]
0.00674330566615
[ 1.  1.  1.  1.  1.  0.  0.  1.  1.  1.  0.  0.  0.  0.  1.  0.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  0.  1.  1.
  1.  1.  1.  1.  1.  1.  0.  1.  0.  0.  1.  1.  1.  0.]
0.00674330566615
[ 1.  1.  1.  1.  0.  1.  0.  1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  0.
  1.  1.  1.  1.  1.  0.  0.  1.  1.  1.  1.  0.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  0.  1.  1.  1.  1.  1.  1.  1.  0.  0.  1.]
0.00674330566615
[ 1.  1.  1.  1.  0.  0.  0.  0.  1.  1.  0.  1.  1.  0.  0.  1.  0.  1.
  1.  1.  0.  1.  0.  1.  1.  0.  1.  1.  1.  1.  1.  1.  1.  0.  0.  1.
  0.  1.  1.  1.  1.  1.  1.  0.  1.  0.  1.  1.  1.  1.]
0.00674330566615
[ 1.  1.  0.  1.  1.  1.  0.  0.  1.  1.  0.  1.  1.  0.  0.  1.  1.  0.
  1.  1.  1.  0.  1.  1.  1.  1.  0.  0.  0.  1.  1.  1.  1.  1.  1.  1.
  0.  1.  1.  1.  1.  0.  1.  1.  0.  1.  0.  0.  1.  1.]
0.00674330566615
[ 1.  1.  0.  0.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  1.  0.  1.  1.  1.  1.  1.
  1.  1.  0.  1.  0.  1.  1.  0.  1.  0.  1.  1.  1.  1.]
0.00674330566615
[ 1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  0.  1.  1.  0.  1.  0.  1.  1.
  1.  1.  1.  0.  1.  0.  1.  1.  0.  1.  1.  1.  0.  1.  1.  1.  1.  0.
  0.  1.  1.  1.  0.  1.  1.  0.  1.  1.  0.  1.  1.  1.]
0.00674330566615
[ 1.  1.  1.  0.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  0.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  1.  0.  1.  1.  1.
  0.  1.  1.  1.  1.  1.  1.  0.  1.  1.  0.  1.  1.  1.]
0.00674330566615
[ 1.  1.  0.  1.  1.  1.  1.  0.  1.  1.  1.  1.  1.  1.  0.  1.  0.  1.
  1.  1.  1.  1.  1.  1.  1.  1.  1.  0.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  0.  1.  0.  1.  0.  1.  1.  1.  1.  1.  1.  0.]
</code></pre>

<hr>

<p>Note that there is a much simpler solution if </p>

<ul>
<li><p>your sequence is of length <em>n</em>, </p></li>
<li><p>some number <em>m</em> has <em>m &lt;&lt; n</em>, and </p></li>
<li><p><em>m!</em> is much larger than the number of permutations you need.</p></li>
</ul>

<p>In this case, simply divide your sequence into <em>m</em> (approximately) equal parts, and permute them randomly. As noted before, only the <em>m - 1</em> boundaries change in a way that potentially affects the correlations. Since <em>m &lt;&lt; n</em>, this is negligible.</p>

<p>For some numbers, say you have a sequence with 10000 elements. It is known that <a href="https://en.wikipedia.org/wiki/Factorial" rel="nofollow">20! = 2432902008176640000</a>, which is far far more permutations than you need, probably. By dividing your sequence into 20 parts and permuting, you're affecting at most 19 / 10000 with might be small enough. For these sizes, this is the method I'd use.</p>

python numpy random permutation