cedric's challenge

demo.py

"""
Challenge: http://beust.com/weblog/2012/02/16/a-new-coding-challenge-2/
Author: Praveen Kumar T

The Test:
--------
The probability distribution of the number of fixed points of a 
uniformly distributed random permutation approaches a Poisson 
distribution with expected value 1 as n grows.
he first n moments of this distribution are exactly those of 
the Poisson distribution.

Source: http://en.wikipedia.org/wiki/Random_permutation#Statistics_on_random_permutations
"""

import scipy
import scipy.stats

def seq_producer(n):
    """
    Produces a random shuffled array of [1..n]
    """
    seq = scipy.arange(n)
    scipy.random.shuffle(seq)
    return seq

def count_fixed_points(n, seq):
    """
    Count fixed points.
    """
    orig_seq = scipy.arange(n)
    return sum(p == q for p, q in zip(orig_seq, seq))

seq_size = 100 # sequence size to generate
sample_size = 100 # number of sequences to generate

# raw data of fixed points of the generated sequences
fp_data = [count_fixed_points(seq_size, seq_producer(seq_size)) for _ in range(sample_size)]

# number of moments to verify.
# Higher this value means more the sample data needed for accurate check
cmp_num_moments = 6

# central moments of the fixed point data
seq_moments = [scipy.stats.moment(fp_data, moment=i) for i in range(1, cmp_num_moments)]

# central moments of poisson distribution with mean = 1
# data courtesy: http://www.wolframalpha.com/input/?i=central+moments+of+poisson+distribution+with+mean+1
poisson_moments = [0, 1, 1, 4, 11, 41, 162, 715, 3425, 17722,
                   98253, 580317, 3633280, 24011157, 166888165,
                   1216070380, 9264071767, 73600798037,
                   608476008122, 5224266196935]

# sum the differences between corresponding moments
# Ideally this should be zero for a random sequence
print(sum(abs(x - y) for x, y in zip(seq_moments, poisson_moments)))