In [1]:
import sys
sys.path.append('../sample/')
from metropolis_sampler import MetropolisSampler
from random import uniform, gauss
import random
import numpy as np
import matplotlib.pyplot as plt
In [2]:
def initialize_state(dim):
return np.array([uniform(-10, 10) for i in range(dim)])
def markov_process(x, step_length):
result = x.copy()
for i, item in enumerate(result):
result[i] = item + gauss(0, 1) * step_length
return result
In [19]:
def N(mu, sigma):
""" float * float -> ([float] -> float)
"""
return lambda x: np.exp(- np.sum(np.square((x - mu) / sigma)))
## Recall SimulatedAnnealing is searching the argmin, instead of argmax.
def target_function(x):
""" [float] -> float
"""
return 1 * N(-5, 5)(x) + 10 * N(5, 5)(x)
def log_target_distribution(T):
""" float -> ([float] -> float)
"""
return lambda x: np.log(target_function(x))
In [94]:
def sampling_0(iterations_0, burn_in_ratio, dim, step_length_0, T):
def initialize_state_0():
return initialize_state(dim)
def markov_process_0(x):
return markov_process(x, step_length_0)
ms = MetropolisSampler(iterations_0,
initialize_state_0,
markov_process_0,
int(burn_in_ratio * iterations_0),
log=False,
)
chain_0 = ms.sampling(log_target_distribution(T))
print('Pre- Accept Ratio: {0}'.format(ms.accept_ratio))
return chain_0
def sampling_1(chian_0, iterations_1, burn_in_ratio, dim, step_length_1, T):
def initialize_state_1():
return random.choice(chain_0)
def markov_process_1(x):
return markov_process(x, step_length_1)
chain_1 = []
accept_ratios = []
for i in range(1000):
ms = MetropolisSampler(iterations_1,
initialize_state_1,
markov_process_1,
int(burn_in_ratio * iterations_1),
log=False,
)
chain_1 += ms.sampling(log_target_distribution(T))
accept_ratios.append(ms.accept_ratio)
print('Accept Ratio: {0}'.format(np.mean(accept_ratios)))
return chain_1
In [99]:
dim = 300
## Needs tuning
iterations_0 = int(10 ** 5 * 1)
burn_in_ratio = 0.5
step_length_0 = 0.1
T = 1
chain_0 = sampling_0(iterations_0, burn_in_ratio, dim, step_length_0, T)
In [101]:
iterations_1 = int(1000)
step_length_1 = 0.02
chain_1 = sampling_1(chain_0, iterations_1, burn_in_ratio, dim, step_length_1, T)
print('Lengh of chain: {0}'.format(len(chain)))
expect = [np.mean([state[axis] for state in chain])
for axis in range(dim)
]
print('Expects: {0}'.format(expect))
In [ ]:
steps = np.arange(len(chain))
targets = [target_function(state) for state in chain]
plt.plot(steps, targets)
plt.xlabel('step')
plt.ylabel('target value')
plt.show()
for axis in range(dim):
xs = [state[axis] for state in chain]
plt.plot(steps, xs)
plt.xlabel('step')
plt.ylabel('x[{0}]'.format(axis))
plt.show()
Splendid, even for dim = 10.