学步园

参考资料：http://www.52nlp.cn/lda-math-mcmc-%e5%92%8c-gibbs-sampling1

python脚本：

import math
import time
import numpy as np
import matplotlib.pylab as plt
import random
class Samples:
    def __init__(self):
        pass
    def rand(self, num, seed = 1):
        m = math.pow(2, 32)
        a = 214013
        c = 2531011
        i = 1
        x = np.zeros(num)
        x[0] = seed
        while(i < num):
            x[i] = (a * x[i-1] + c) % m
            i += 1
        return x
    def uniform(self, num, seed = 1):
        m = math.pow(2, 32)
        x = self.rand(num, seed)
        return x / m
    def normal(self, num):
        t = int(time.time()) 
        u1 = self.uniform(num, t)
        t1 = int(time.time())+1
        u2 = self.uniform(num, t1)
        z0 = np.sqrt(-2 * np.log(u2)) * np.sin(2.0 * np.pi * u1)
        #plt.hist(z0)
        #plt.show()
        return z0
if __name__=='__main__':
    s = Samples()
    x00 = s.normal(1000)
    plt.hist(x00,100)
    plt.show()

样本的直方图如下：