public
Created

cedric's challenge

  • Download Gist
demo.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
"""
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)))

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.