In [ ]:
"""
This module implements Metropolis-Hastings algorithm for random variable
generation. The algorithm generates random variables from a desired
distribution (which may be unnormalized).
"""
def metropolis( desiredPDF, initValue, computableRVS, skipIterations = 1500 ):
"""
This function returns a generator, which generates random variables
from some space S with a desired distribution using Metropolis-Hastings
algorithm.
Args:
desiredPDF (func) : PDF of desired distribution p( T ), where T from S
initValue : an object from S to initialize the starting point
of iterative proccess
computableRVS (func) : a generator of random value from space S
with given parameter T, which is also from S
skipIterations (int) : number of iterations to skip
(skipping more iterations leads to better accuracy,
but more time consuming)
Returns: generator, which produce some values from S
and their denisity according to distribution desiredPDF
"""
random_variable = initValue
random_variableDensityValue = desiredPDF( random_variable )
"""
A state of MCMC
"""
#ignore first iterations to let the iterative proccess
#converge to some distribution, which is close to desired
for i in xrange( skipIterations ):
candidate = computableRVS( random_variable )
candidateDensityValue = desiredPDF( candidate )
"""
next candidate for sample, generated by computableRVS
"""
#acceptanceProb = min( 1, candidateDensityValue / random_variableDensityValue )
#logp is returnd by desiredPDF_, so here is the change
acceptanceProb = min(0, candidateDensityValue - random_variableDensityValue )
"""
probability to accept candidate to sample
"""
if math.log(random.random()) < acceptanceProb:
random_variable = candidate
random_variableDensityValue = candidateDensityValue
#now when the procces is converged to desired distribution,
#return acceptable candidates
while True:
candidate = computableRVS( random_variable )
candidateDensityValue = desiredPDF( candidate )
"""
next candidate for sample, generated by computableRVS
"""
#acceptanceProb = min( 1, candidateDensityValue / random_variableDensityValue )
# logp is returnd by desiredPDF_, so here is the change
acceptanceProb = min( 0, candidateDensityValue - random_variableDensityValue )
"""
probability to accept candidate to sample
"""
if math.log(random.random()) < acceptanceProb:
random_variable = candidate
random_variableDensityValue = candidateDensityValue
yield random_variable, random_variableDensityValue
In [ ]:
"""
This module provides some functions
that generate random permutations with different distributions.
There are a uniform distribution and a symmetric distribution,
which depends on some other permutation.
"""
def uniform( n ):
"""
Generates random permutation using Knuth algorithm.
Args:
n (int) : length of permutation
Returns: random permutation of length n from uniform distribution
"""
#initialize permutation with identical
permutation = [ i for i in xrange( n ) ]
#swap ith object with random onject from i to n - 1 enclusively
for i in xrange( n ):
j = random.randint( i, n - 1 )
permutation[ i ], permutation[ j ] = permutation[ j ], permutation[ i ]
permutation.append(26)
return permutation
def applyTransposition( basePermutation ):
"""
This function returns random permutation by applying random transposition
to given permutation.
The result distribution is not uniform and summetric assuming parameter.
Args:
basePermutation (array) : parameter of distribution
Returns: random permutation generated from basePermutation
"""
n = len( basePermutation )
"""
length of permutation
"""
permutation = copy( basePermutation )
"""
permutation to return after some modifications
we use a copy method, because initial arguments must be left unchanged
"""
#apply n random transpositions (including identical) to base permutation
#for i in xrange( n ):
k, l = random.randint( 0, n - 1 ), random.randint( 0, n - 1 )
permutation[ k ], permutation[ l ] = permutation[ l ], permutation[ k ]
return permutation