Is there an efficient algorithm that allows to obtain a uniform distribution of all possible $4$-bit permutations from a single keyless $4$-bit S-Box?
In Haskell the recursive function

f s x = s (g 16 (x `div` 16) (x `mod` 16))
  where g 1 _ = id
        g n x = t (n-1) (x `mod` n) . g (n-1) (x `div` n)
        t x y z = if x == z then y
                            else if y == z then x else z


runs for any function s that is a permutation on [0..15] through all permutations on [0..15] if you let x run through all numbers 0...20922789888000*16-1.

Of course, s doesn't matter much, so I'd set s = id, and look at the definition of g using t: t x y is simply exchanging the values x and y, keeping all other values the same.

Now g n x is a(n arbitrary!) permutation on [0..n-1] depending on x: x mod n determines which element to exchange with n-1 (all other elements are untouched) and x divided by n is used as parameter to for defining recursively g (n-1).

If you pick n uniformly from 0..20922789887999, the permutation g n will be a uniform random element of the symmetric group on [0..15].