Sampling from a normal distribution by Metropolis Monte Carlo

ms_normal.py

import sys
import math
import random

def invariant_distribution(x):
    return math.exp(- 0.5 * x ** 2 / (0.3 * 0.3))

def next_state(x, dx):
    return x + dx * (random.uniform(0.0, 1.0) - 0.5)

size = 100000
dx = 1.0
now = 0.0
results = []
for i in range(size):
    next = next_state(now, dx)
    probability = min(1.0,
                      invariant_distribution(next)
                      / invariant_distribution(now))
    if probability > random.uniform(0.0, 1.0):
        now = next
        results.append(now)

step = 0.01
freq = {}
for result in results:
    res = step * round(result / step)
    if res == -0.0:
        res == 0.0
    freq[res] = freq.get(res, 0) + 1

total = 0.0
for res in sorted(freq):
    total += freq[res]
for res in sorted(freq):
    print "%.2lf\t%.6lf" % (res, freq[res] / (total * step))