example1.py

import copy
import math
import random

def random_x():
    return 2*random.randint(0, 1) - 1

def modify(x):
    i = random.randint(0, 2)
    xr = copy.copy(x)
    xr[i] = -xr[i]
    return xr

def p_ratio(x, theta):
    return math.exp(-2.0*theta*x[0]*(x[1] + x[2]))

theta = 1.0
range_max = 10000

x = [random_x(), random_x(), random_x()]
x1_sum = 0
x1x2_sum = 0
for i in range(range_max):
    r = p_ratio(x, theta)
    R = random.random()
    if R < r:
        x = modify(x)
    print ('x = %s'%x)
    x1_sum += x[0]
    x1x2_sum += x[0]*x[1]

print('average x1 = %f' % (x1_sum / range_max))
print('average x1*x2 = %f' % (x1x2_sum / range_max))