In [2]:
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
In [3]:
# define distribution function

def gauss(x, mu, sigma):
    return 1./sqrt(2*pi*sigma**2)*exp(-(x-mu)**2/(2*sigma**2))

def P(x):
    mu1, sigma1 = 0.1, 0.1
    mu2, sigma2 = 0.5, 0.07
    return gauss(x, mu1, sigma1)*0.3 + gauss(x, mu2, sigma2)*0.7

x = linspace(-1,1, 200)
plot(x, P(x)) 
Out[3]:
[<matplotlib.lines.Line2D at 0x7f064ec2a7f0>]
In [4]:
def step(x, stepsize = 1):
    return x + random.normal(0, stepsize)
In [5]:
def metropolis(x0, step, nstep, P, stepsize):
    path = zeros((nstep,))
    xold = x0
    path[0] = x0
    accept = 0

    for i in range(1,nstep):
        xnew = step(xold, stepsize)
        r = P(xnew)/P(xold)
        A = min([1,r])
        if A > random.uniform():
            xold = xnew

            accept += 1
        path[i] = xold
    return path, float(accept)/nstep
In [6]:
def step(x, stepsize = 0.1):
    return x + random.normal(0, stepsize)

def plot_path():
    figure()
    ax1 = gcf().add_axes([.15, .15, .55, .8])
    ax2 = gcf().add_axes([.70, .15, .25, .8])
    #ax2 = subplot(122)
    ax2.axes.yaxis.set_ticklabels([])
    ax2.axes.xaxis.set_ticklabels([])
    ax1.plot(path)
    ax2.hist(path, bins=linspace(-1,1,100), normed=1, orientation="horizontal", histtype="stepfilled", color="#DDDDDD")
    ax2.set_ylim([-0.4, 0.8])
    ax1.set_ylim([-0.4, 0.8])
    ax1.set_title("stepsize = %.3f, acceptance = %.3f" % (stepsize, acceptance))

    x = linspace(-1,1, 200)
    ax2.plot(P(x),x,"r", lw=2)

random.seed(12)
for stepsize in [0.01, 0.1, 2.]:
    path, acceptance = metropolis(0., step, 1000, P, stepsize=stepsize)
    plot_path()
pass
In [7]:
def E(x):
    return x**2 - cos((cos(x*4)**8-1)*20)
    
    #return x**2-cos(x*gauss(x, 0, 0.3)*200)*gauss(x, 0, 0.5)*0.4
    #return x**2-cos(x/(0.5+x**2)*200)*gauss(x, 0, 0.5)*0.4

x = linspace(-1,1, 1000)
plot(x, E(x))
Out[7]:
[<matplotlib.lines.Line2D at 0x7f064e9e59b0>]
In [10]:
def greedy(x0, step, nstep, E, stepsize):
    path = zeros((nstep,))
    xold = x0
    path[0] = x0
    accept = 0
    for i in range(1,nstep):
        xnew = step(xold, stepsize)
        Pnew = P(xnew)
        if E(xnew) < E(xold):
            xold = xnew
            accept += 1
        path[i] = xold
    return path, float(accept)/nstep
In [14]:
random.seed(1234)
def step(x, stepsize = 0.1):
    return x + random.normal(0, stepsize)

for stepsize in [0.01, 0.5, 5.]:
    path, acceptance = greedy(0.8, step, 2000, E, stepsize=stepsize)
    figure(figsize=(6,4))
    ax1 = gcf().add_axes([.15, .15, .45, .8])
    ax2 = gcf().add_axes([.70, .15, .35, .8])

    ax1.plot(path, E(path))
    x = linspace(-1,1, 1000)
    ax1.plot(x, E(x))
    ax2.semilogy(E(path)-E(0))
pass
In [15]:
def boltzmann_anneal(x0, step, nstep, E, schedule, stepsize):
    path = zeros((nstep,))
    xold = x0
    path[0] = x0
    accept = 0
    for i in range(1,nstep):
        T = schedule(i)
        xnew = step(xold, i, stepsize)
        r = exp(-(E(xnew)-E(xold))/T)
        A = min(array([1,r]))
        if A > random.uniform():
            xold = xnew
            accept += 1
        path[i] = xold
    return path, float(accept)/nstep
In [16]:
T0 = 1.
def schedule(k):
    return T0/log(k+1)

def step(x, k, stepsize = 0.1):
    return x + random.normal(0, stepsize*sqrt(schedule(k)))

random.seed(1234)
nstep = 10000
x0 = 0.8
for stepsize in [0.01, 0.1, 5.]:
    path, acceptance = boltzmann_anneal(x0, step, nstep, E, schedule, stepsize=stepsize)
    figure(figsize=(12,4))
    ax1 = gcf().add_axes([.10, .15, .40, .8])
    ax2 = gcf().add_axes([.55, .15, .20, .8])
    ax3 = gcf().add_axes([.80, .15, .20, .8])

    ax1.plot(path, E(path))
    x = linspace(-1,1, 1000)
    ax1.plot(x, E(x))
    ax2.semilogy(E(path)-E(0))
    ax3.plot(schedule(arange(nstep)))
pass
-c:3: RuntimeWarning: divide by zero encountered in true_divide
In [17]:
T0 = 10.
def schedule(k):
    return T0/(k+10)

def step(x, k, stepsize = 0.1):
    return x + random.normal(0, stepsize*sqrt(schedule(k)))

random.seed(1234)
nstep = 10000
x0 = 0.8
for stepsize in [0.1, 1, 10.]:
    path, acceptance = boltzmann_anneal(x0, step, nstep, E, schedule, stepsize=stepsize)
    print(acceptance)
    figure(figsize=(12,4))
    ax1 = gcf().add_axes([.10, .15, .40, .8])
    ax2 = gcf().add_axes([.55, .15, .20, .8])
    ax3 = gcf().add_axes([.80, .15, .20, .8])

    ax1.plot(path, E(path))
    x = linspace(-1,1, 1000)
    ax1.plot(x, E(x))
    ax2.semilogy(E(path)-E(0))
    ax3.plot(schedule(arange(nstep)))
pass
0.4687
0.1205
0.013
In [18]:
T0 = 10.
def schedule(k):
    return T0/(k+1)

def step(x, k, stepsize = 0.1):
    return x + random.standard_cauchy() * stepsize*(schedule(k))

random.seed(1234)
nstep = 10000
x0 = 0.8
for stepsize in [0.1, 1, 10.]:
    path, acceptance = boltzmann_anneal(x0, step, nstep, E, schedule, stepsize=stepsize)
    print(acceptance)
    figure(figsize=(12,4))
    ax1 = gcf().add_axes([.10, .15, .40, .8])
    ax2 = gcf().add_axes([.55, .15, .20, .8])
    ax3 = gcf().add_axes([.80, .15, .20, .8])

    ax1.plot(path, E(path))
    x = linspace(-1,1, 1000)
    ax1.plot(x, E(x))
    ax2.semilogy(E(path)-E(0))
    ax3.plot(schedule(arange(nstep)))
pass
0.8643
0.6665
0.1934
In [19]:
T0 = 10.
def schedule(k):
    eps=0.1
    return T0*exp(-k*eps)

def step(x, k, stepsize = 0.1):
    return x + random.standard_cauchy() * stepsize*sqrt(schedule(k))

random.seed(1234)
nstep = 2000
x0 = 0.8
for stepsize in [0.1, 1., 10., 100., 1., 1., 1.]:
    path, acceptance = boltzmann_anneal(x0, step, nstep, E, schedule, stepsize=stepsize)
    print(acceptance)
    figure(figsize=(12,4))
    ax1 = gcf().add_axes([.10, .15, .40, .8])
    ax2 = gcf().add_axes([.55, .15, .20, .8])
    ax3 = gcf().add_axes([.80, .15, .20, .8])

    ax1.plot(path, E(path))
    x = linspace(-1,1, 1000)
    ax1.plot(x, E(x))
    ax2.semilogy(E(path)-E(0))
    ax3.plot(schedule(arange(nstep)))
pass
0.7115
0.863
0.7985
0.7595
0.866
0.874
0.8615
-c:9: RuntimeWarning: overflow encountered in exp