gistfile1.py

#!/usr/bin/python
# -*- coding: utf-8 -*-


from math import sqrt, log, cos, sin, pi


def kai2(values, expected):
    x = 0
    for value in values:
        x += ((value - expected) ** 2.0) / expected
    return x


def fd(count, n, k, randargs=()):
    ''' Frequency Distribution '''
    for i in xrange(k):
        v = count[i]
        print '%3d[%3d]: %s' % (i, v, '*' * (v / 5))
    kai = kai2(count, n / k)
    print 'χ2: %f' % (kai)


seed = 39484


def lcg(setseed=False):
    ''' Linear Congruential Generator '''
    global seed
    if setseed != False:
        seed = setseed
    A = 33
    B = 729
    M = 1024
    seed = (A * seed + B) % M
    return seed * 1.0 / M


def boxmuller(sigma, mu):
    '''
        r1, r2 = 独立した二乱数
        x1 = σ * √(-2 ln r1) * cos(2π r2) + mu
        x2 = σ * √(-2 ln r2) * sin(2π r2) + mu
    '''
    r1, r2 = lcg(), lcg()
    x1 = sigma * sqrt(-2 * log(r1)) * cos(2 * pi * r2) + mu
    x2 = sigma * sqrt(-2 * log(r1)) * sin(2 * pi * r2) + mu
    return (x1, x2)


def main():
    k = 10
    n = 1000

    # 一様乱数（線形合同法）による分布図
    count = [0] * k
    for _ in xrange(n):
        r = int(lcg() * k)
        count[r] += 1
    fd(count, n, k)

    # 正規乱数（ボックス・ミューラー法）による分布図
    count = [0] * k
    for _ in xrange(n):
        x, y = boxmuller(1, 4)
        count[int(x)] += 1
        count[int(y)] += 1
    fd(count, n, k)


if __name__ == '__main__':
    main()