What's wrong with the following shuffle algorithm?

Question

Suppose I have an array A of size n, with the initial state of:

A[0] == 1, A[1] = 2, ... A[n-1] = n

I know that one way to get a uniformly-distributed permutation of that array is to use Fisher-Yates Algorithm, but I am more interested in what's flawed with the following naive approach:

for i = 0 to n-1:
   rand_i = random(0, n-1, UNIFORM_DIST)
   swap(A[i], A[rand_i])

When I tried to use that with n = 3, it seems like some permutations are more likely to be drawn than others.
Here's the output of running the above algorithm 10,000,000 times, and then averaging on the results (shown in percentage of likelihood):

#> swap.out --iteartions 10000000
00: (1,2,3): 14.8237
01: (1,3,2): 18.5105
02: (2,1,3): 18.504
03: (2,3,1): 18.5216
04: (3,1,2): 14.7975
05: (3,2,1): 14.8427

To rule out the possibility that the Pseudo-Random-Generator is not uniform, I compared that with an algorithm that builds db of all permutations and then, using the same PRG - chooses one permutation:

#> choose_permutation.out --iterations 10000000
00: (1,2,3): 16.6847
01: (1,3,2): 16.6649
02: (2,1,3): 16.6731
03: (2,3,1): 16.6706
04: (3,1,2): 16.6516
05: (3,2,1): 16.655

This pattern is consistent; permutations 1,2,3 always come up with a higher likelihood than 0,4,5.
With n = 4:

#> swap.out --iteartions 10000000
00: (1,2,3,4): 3.90774
01: (1,2,4,3): 3.90958
02: (1,3,2,4): 3.91321
03: (1,3,4,2): 5.46167
04: (1,4,2,3): 4.29965
05: (1,4,3,2): 3.51932
06: (2,1,3,4): 3.89793
07: (2,1,4,3): 5.83975
08: (2,3,1,4): 5.45905
09: (2,3,4,1): 5.47537
10: (2,4,1,3): 4.30715
11: (2,4,3,1): 4.3
12: (3,1,2,4): 4.29691
13: (3,1,4,2): 4.3075
14: (3,2,1,4): 3.51411
15: (3,2,4,1): 4.29534
16: (3,4,1,2): 4.30452
17: (3,4,2,1): 3.90896
18: (4,1,2,3): 3.12621
19: (4,1,3,2): 3.50413
20: (4,2,1,3): 3.52526
21: (4,2,3,1): 3.12843
22: (4,3,1,2): 3.89564
23: (4,3,2,1): 3.90257

And comparing to randomly choosing from db of permutations:

#> choose_permutation.out --iterations 10000000  
00: (1,2,3,4): 4.16284
01: (1,2,4,3): 4.16721
02: (1,3,2,4): 4.15893
03: (1,3,4,2): 4.17306
04: (1,4,2,3): 4.15853
05: (1,4,3,2): 4.16169
06: (2,1,3,4): 4.16584
07: (2,1,4,3): 4.17245
08: (2,3,1,4): 4.17309
09: (2,3,4,1): 4.15519
10: (2,4,1,3): 4.17007
11: (2,4,3,1): 4.17163
12: (3,1,2,4): 4.16276
13: (3,1,4,2): 4.17367
14: (3,2,1,4): 4.17147
15: (3,2,4,1): 4.16955
16: (3,4,1,2): 4.16576
17: (3,4,2,1): 4.1659
18: (4,1,2,3): 4.16653
19: (4,1,3,2): 4.16749
20: (4,2,1,3): 4.17225
21: (4,2,3,1): 4.1668
22: (4,3,1,2): 4.16797
23: (4,3,2,1): 4.15932

What's flawed with that method and why are those specific permutations come up more often?

Answer 1

This approach cannot work, for the following simple reason. The probability to obtain any permutation is of the form $A/n^n$, for integer $A$. However, we need it to be $1/n!$, so we need $A = n^n/n!$. Unfortunately, for $n \geq 3$ this is not an integer.

---

One can still ask why we might expect this algorithm to work, and where is the gap between our intuition and reality.

After the $i$'th step of the algorithm, $A[i]$ is a uniformly random element. So if we randomize $A[0]$, then $A[1]$, then $A[2]$, and so on until $A[n-1]$, won't we get a random permutation?

While it is true that $A[n-1]$ will be completely random, this is not the case for other entries. Indeed, let us try to analyze the distribution of $A[0]$ as the algorithm progresses.

After the first iteration, $A[0]$ is uniformly random.

After the second iteration, $A[1]$ is uniformly random, but this comes at a cost: we might have swapped $A[0]$ and $A[1]$. This means that $A[0]$ will have a slight preference towards $1$ (the original value of $A[1]$) and $0$ (the value of $A[1]$ in case it was swapped with $A[0]$).

Quantitatively, after the second iteration, if $i > 1$ then $$ \Pr[A[0] = i] = \frac{1}{n} \cdot \left(1 - \frac{1}{n}\right) = \frac{1}{n} - \frac{1}{n^2}, $$ since this happens when $A[0]$ was swapped with $A[i]$, and $A[1]$ was not swapped with $A[0]$. Similarly, $$ \Pr[A[0] = 0] = \frac{1}{n} \cdot \left(1 - \frac{1}{n}\right) + \frac{1}{n} \cdot \frac{1}{n} = \frac{1}{n}, $$ accounting for the additionally possibility in which $A[1]$ is swapped with $A[0]$ in the second round. Finally, $$ \Pr[A[0] = 1] = \frac{1}{n} \cdot \left(1 - \frac{1}{n}\right) + \left(1 - \frac{1}{n}\right) \cdot \frac{1}{n} = \frac{2}{n} - \frac{2}{n^2}, $$ accounting for the possibility that $A[0]$ was not swapped with $A[1]$, but $A[1]$ was swapped with $A[0]$.

Similar phenomena persist throughout the algorithm.