I am reading a mathematical article that contains the following lemma:
Given any sequence of real numbers a_{1}, ..., a_{n}, let us denote their cumulative sums as s_{1}, ..., s_{n}. There exists a unique cyclic permutation sigma of a_{1}, ..., a_{n} such that
for all k: s_{k}(sigma) <= s_{n}/n * k
(i.e., the sequence of partial sums lies below the line connecting (0, 0) and (n, s_{n}) ).
I decided to write a Python code to illustrate this lemma (that would plot the graph of the sequence of partial sums for a random sequence, together with the graph of the sequence of partial sums for the unique permutation, the existence of which is guaranteed by the lemma). Here is the code:
randarr = np.random.randn(9)
tot_arr = np.array([0, *randarr])
init_cumsum = np.cumsum(tot_arr)
mean_array = np.array([ n * init_cumsum[-1]/(len(init_cumsum) - 1)\
for n in range(len(init_cumsum)) ])
def permuted(k):
return np.cumsum([0, *np.roll(randarr, k)])
fig, ax = plt.subplots(figsize=(5, 3))
ax.plot(np.arange(0, 10, 1), init_cumsum)
ax.plot(np.arange(0, 10, 1), mean_array, '--')
for k in range(1, 10):
if (mean_array >= permuted(k)).all() :
ax.plot(np.arange(0, 10, 1), permuted(k), label="correct permutation")
plt.legend()
break
else:
print(tot_arr)
plt.show()
Since I am a simple man, instead of constructing the right permutation as described by the lemma, I just iterate over all cyclic permutations (except for the initial one), and plot the one that satisfies the condition. However, sometimes Python thinks no cyclic permutation satisfies the condition (i.e., the initial permutation is not the right one, but the code can't find the right one among the rest).
I made it print such random sequences, and then tried inserting them instead of the random ones. When I do that, the code readily finds the right permutation and plots it. I can't see the difference between the data obtained by the rest of the code the first time such sequence is generated, and in the following runs when I input it manually, so, I guess, something in my code works not the way I expect, and this explains why I obtain different results with the same input data. However, I can't understand what could be wrong.
EDIT:
Examples of sequences on which the program breaks:
-0.12254795 -1.18301478 1.9911016 -0.67975716 -0.53877141 0.49170832 -1.08431192 1.03930656 -0.14305318
-0.98415686 -0.19698914 1.25619013 -0.31378433 1.45470612 0.04794845 -2.12723127 -0.69520889 0.28457194
2.11822629 -1.2222578 0.48015473 -0.57059476 0.76872146 0.44806995 -0.76466846 0.404127 -0.35847348
-0.04126626 -2.0829324 -0.5147483 0.13324286 -2.08391787 0.00566639 2.07952122 0.83424157 0.71345859
Also, the program appears to always work when I use np.random.randint(-10, 10, 9) instead of np.radnom.randn(9)
I think this might be due to a rounding error in mean_array, since subtracting a very small term (eg 0.1**15 which is 10 to the power of -15) as per this if statement produces a solution;
for k in range(1, 10):
if (mean_array >= permuted(k)).all() - 0.1**15 :
ax.plot(np.arange(0, 10, 1), permuted(k), label="correct permutation")
plt.legend()
Possibly rounding the mean_array to say (5 to 10) decimal places may also help.