In [1]:

    

import sys
sys.path.append('../sample/')
from metropolis_sampler import MetropolisSampler

from random import uniform, gauss
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 [6]:

    

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 [37]:

    

def sampling(iterations, burn_in, dim, step_length, T):
    
    def initialize_state_0():
        return initialize_state(dim)
    
    def markov_process_0(x):
        return markov_process(x, step_length)

    ms = MetropolisSampler(iterations,
                           initialize_state_0,
                           markov_process_0,
                           burn_in,
                           log=False,
                           )
    chain = ms.sampling(log_target_distribution(T))
    print('Accept Ratio: {0}'.format(ms.accept_ratio))
    
    return chain


    

In [39]:

    

dim = 300

## Needs tuning
iterations = int(10 ** 5 * 1)
burn_in = int(0.5 * iterations)
step_length = 0.05
T = 1


chain = sampling(iterations, burn_in, dim, step_length, 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))


    

    




Accept Ratio: 0.89618
Lengh of chain: 50001
Expects: [3.0067161074918647, 6.2320824968266413, 4.986068128271226, 7.1327881876491857, 5.8452986489863026, 6.0829959451377205, 4.868576040734486, 4.5672060695001546, 2.9560042930516257, 5.8602425087711589, 8.0738318383851446, 6.6036025174869986, 2.3889684166855769, 4.9707734440606064, -0.97407990435190694, 3.6818177915237378, 3.7017384242656526, 6.4885749573466054, 5.4438147828511712, 6.9021638635202276, 5.014819486442355, 2.5906733895744836, 2.9245294428785256, 6.0448518131563942, 6.0710505251056333, 4.5038066086773458, 3.9921356480775532, 4.2976243814341704, 6.448576131779526, 8.4715295077643091, 4.9759081171604205, 2.4788079878479072, 12.29599405383243, 3.5538192012603003, 6.4561750564408946, 6.116800895331159, 4.431466625505184, 6.729569687495303, 9.7249446001846014, 3.4638984282780858, 4.2957008856833481, 8.2860009908676098, 4.5480736616982043, 3.5149210797312764, 6.2300795475775246, 4.2804158574898663, 3.0569641521163033, 2.1369476472653073, 6.9660159965772719, 5.7294072902408004, 4.7809522661524406, 5.5938481207297537, 5.9626786422324018, 6.191059801800666, 8.3725761791425022, 5.9404538638253905, 2.6410237907316723, 5.1463895051799611, 5.7481097842293964, 4.8423627973760359, 1.7097601624375893, 7.6376866157521688, 4.4463796322116487, 2.8476138081730786, 7.4437392612233566, 6.5547790490486992, 5.3065051744966922, 5.3276464708249325, 8.7096918355619799, 4.8975037665576657, 7.3015965554339815, 2.3123642311078361, 4.4182693618550157, 3.1376946729740234, 4.9493577608713615, 6.1268616681790107, 4.3224255913833343, 1.3981461419351864, 8.0695571501986976, 3.4685978282533365, 6.4846734180108037, 3.5678888342937034, 6.9115653732047262, 2.679074234636603, 6.1960485716381708, 4.1251916704795217, 3.8534911720004765, 6.1677473464615327, 7.0890174234386203, 6.489948507309184, 6.912734763918384, 2.8715570736804445, 5.9027440047948154, 3.1664837317196386, 3.0866938441286416, 3.3776125223995943, 4.3506418918651635, 4.3526828889067017, 3.8776334496971643, 4.3318563252320708, 5.4987014141603687, 1.8582781643790578, 4.9013947648396856, 6.6429723143323685, 4.82567887818084, 5.9791417916811973, 6.4058535248940993, 6.6556861024869507, 4.2222375457337682, 5.4824403341228507, 6.7624331991376803, 4.2158799803991709, 1.6020434056028003, 2.8502058288733538, 4.47362740753253, 6.0106576476946634, 2.4164348435437581, 7.3578982660726293, 5.3220273943340697, 7.5293166942922776, 6.4546128169248496, 4.0139074268017527, 3.8696384253122229, 2.3567076146976462, 2.7530608885606869, 2.5390719673158548, 8.4221487792991621, 0.71061111600079752, 3.9701339965805995, 6.4138165540612073, 8.6855209974699878, 2.1494492187962395, 3.7282229599609922, 4.3880100838527216, 0.12772388320953207, 5.072633225614279, 5.1689516364577184, 0.79010174539701161, 4.1152232781704754, 5.2020167918201423, 6.0309453918034652, 4.443173949248278, 3.5923920605011506, 2.6830390774764235, 6.2487614112648178, 5.6499400978832401, 5.8398746491297278, 3.3784338438985952, 0.55524964470912364, 3.6795403016485366, 4.5132994432776838, 1.2816608621435057, 10.994629969714721, 3.6347824248240257, 3.4895702284930019, 6.5239386366675101, 6.258743274681529, 4.4702567149611374, 2.1357906734145935, 2.7535704182783758, 4.0900991613314437, 3.865340327984025, 6.0512503430968092, 3.8985690259583268, 0.45278273955076742, 4.5619045924724704, 5.1890762827984309, 3.0001573777216852, 4.2658808976835116, 2.5722762695461605, 5.3756731674615734, 3.6585708377791248, 2.41725835801723, 1.7453800730838192, 3.2354796714287133, 4.4242364969404822, 7.9503468850605348, 3.2931243936261545, 6.063037943667096, 6.4007044634910217, 5.2227437826634047, 7.1913702740434342, 4.564953910970627, 6.6582484641022495, 5.7037587166730264, 5.9203262613139316, 6.6755482249946079, 4.2895828605342485, 2.259993593007481, 4.6363014278562753, 3.8032127195471253, 6.9784560595121494, 2.6131275326412351, 10.043570432175622, 2.7132154400914148, 4.8402554720131574, 4.4923425612067467, 4.5578271592314232, 3.0593172647994695, 5.915410643272967, 4.8798458660402151, 5.3829176994953647, 8.4552480083935126, 4.0112603108039782, 2.7448972072115705, -0.56911055691365575, 3.5162317918467694, 6.4726916462716355, 6.6692413125604419, 5.6450667811513213, 5.0556316444250244, 5.3449451337842424, 1.6838249992668455, 8.2159185452626708, 2.1176710607158693, 6.994516499422601, 7.3280574793453175, 3.6744308981945575, 4.0082065039232964, 4.1990611988248236, 5.3707404204163005, 6.4727373943121851, 3.0954211288762039, 5.4034568407427752, 3.0055210370034802, 2.6241255684005127, 7.4727630730761119, 6.5600834672394903, 4.6782959486293274, 7.7984330152695325, 5.8983822314911798, 5.2934286977859415, 4.31380354300238, 6.7801609626062618, 6.7128753385224984, 0.27747219962902719, 8.0156606159932053, 5.8816485006264445, 6.9637737819312786, 7.5301325101940311, 6.4687680311307689, 7.0358775019985673, 2.8171320605943753, 3.079487674438647, 6.5030825888047117, 5.390272683657809, 3.4131570211855715, 3.2213916326127188, 6.6603378854787945, 3.3181182103235178, 3.2322915517524478, 4.4781003699695745, 3.7698498778564868, 6.6406497940763449, 7.7850739766801622, 2.4934010045427053, 8.325217880701425, 3.334610043494779, 6.4289796908956243, 5.2599572005333233, 3.4519511285314834, 4.1781017529245572, 7.2453445001360457, 2.519913119343665, 5.3638865763238464, 1.9496686231016274, 5.8949893681948655, 4.2640939284557069, 8.453634912868381, 4.4165826980306182, 3.6440034752845825, 3.4309728300938334, 5.7695023649847013, 7.9678275773166503, 7.2424331261099413, 3.2489482958900471, 3.3498791017172955, 8.7396374006245789, 8.1700881271347789, 4.4209241013809182, 5.5180921675460102, 5.8431480899469799, 5.7955714287720212, 4.1917096879724465, 3.9563099050486321, 3.1386107177582843, 3.2842882343441744, 3.084837362722193, 4.0582913710804034, 4.5651889700183581, 6.43770413354697, 1.5408071868934308, 1.6112248128829656, 1.8092141860097732, 1.236478430343346, 6.4169570178102386, 7.580586428365911, 3.8972297955856372, 3.8955701444638726, 4.5054189371945768]

In [11]:

    

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()