2010年12月12日星期日

序贯蒙特卡洛(二)代码

在实验过程中发现,M2并没有比M3、M4差太多。
M3和M4等价是因为转置不改变矩阵的积和式。物理上可以理解:每一步(一个时间)找一个地方走,和每个地方找一步(一个时间)走,其实都是一种置换。
M1虽然计算复杂度低,但是方差实在太大了,没法要。

import math
import random
import itertools
import pylab
import numpy

def weighted_choice(weights):
    cum = numpy.cumsum(weights)
    return numpy.searchsorted(cum, numpy.random.rand() * cum[-1]), cum[-1]

def restrict(A,p):
    return all([A[i][p[i]] for i in xrange(len(A))])

def make_zero_one_matrix_with_given_fillrate(dimension, fillrate):
    return numpy.random.random((dimension, dimension)) < fillrate

def perm(A):
    return sum([restrict(A, p) for p in itertools.permutations(range(len(A)))])

def eperm1(A, m):
    c=0
    p = range(len(A))
    for _ in xrange(m):
        random.shuffle(p)
        c += restrict(A, p)
    return math.factorial(len(A))*float(c)/m

def eperm2(A, m):
    Backup=A.copy()
    w=numpy.ones(m)
    for trip in xrange(m):
        A=Backup.copy()
        for step in xrange(len(A)):
            candidates = filter(lambda k : A[step][k], xrange(len(A)))
            if len(candidates) == 0:
                w[trip] = 0
                break
            k = random.choice(candidates)
            A[:,k]=False
            w[trip] *= len(candidates)
    return sum(w)/m

def eperm3(A, m):
    Backup=A.copy()
    w=numpy.ones(m)
    for trip in xrange(m):
        A=Backup.copy()
        csum=A.sum(0)
        for step in xrange(len(A)):
            for k, s in enumerate(csum):
                if s == 1:
                    A[:,k]=False
                    break #you have no other choice but goto k, w is kept
            else:
                candidates = filter(lambda k : A[step][k], xrange(len(A)))
                if len(candidates) == 0:
                    w[trip] = 0
                    break
                p = 1.0 / (csum[candidates] - 1)
                k, s = weighted_choice(p)
                w[trip] *= s/p[k]
                A[:,candidates[k]]=False
    return sum(w)/m

def eperm4(A, m):
    Backup=A.copy()
    w=numpy.ones(m)
    for trip in xrange(m):
        A=Backup.copy()
        rsum=A.sum(1)
        for pos in xrange(len(A)):
            for step, s in enumerate(rsum):
                if s == 1:
                    A[step,:]=False
                    break
            else:
                candidates = filter(lambda i : A[i][pos], xrange(len(A)))
                if len(candidates) == 0:
                    w[trip] = 0
                    break
                p = 1.0 / (rsum[candidates] - 1)
                k, s = weighted_choice(p)
                w[trip] *= s/p[k]
                A[candidates[k], :]=False
    return sum(w)/m
   

def main():
    n=15
    K=10
    m=100
    A=make_zero_one_matrix_with_given_fillrate(n,0.7)
    print(A)
    print(n,math.factorial(n),m)
#    print(perm(A))
    x=numpy.zeros((K,4))
    for i in range(K):
        x[i][0]=eperm1(A.copy(), m)
        x[i][1]=eperm2(A.copy(), m)
        x[i][2]=eperm3(A.copy(), m)
        x[i][3]=eperm4(A.copy(), m)
    print(pylab.average(x,0))
    print(pylab.var(x,0))
#    pylab.figure()
#    pylab.boxplot(x)
#    pylab.show()
   
main()

没有评论: