%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
# 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))
[<matplotlib.lines.Line2D at 0x7fa3dfbcd510>]
def step(x, stepsize = 1):
return x + random.normal(0, stepsize)
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
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
/usr/local/lib/python2.7/dist-packages/matplotlib/axes/_axes.py:6448: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg. warnings.warn("The 'normed' kwarg is deprecated, and has been "
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))
[<matplotlib.lines.Line2D at 0x7fa3dfad0ad0>]
x = linspace(-2,2, 1000)
plot(x, E(x*1e-8)+1)
[<matplotlib.lines.Line2D at 0x7fa3df870750>]
x = linspace(-2,2, 1000)
plot(x, ((x*1e-8)**2-1)+1)
[<matplotlib.lines.Line2D at 0x7fa3df768150>]
def greedy(x0, step, nstep, E, schedule, 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
def plot_optim(optimizer, x0, step, nstep, E, schedule, stepsizes):
for stepsize in stepsizes:
path, acceptance = optimizer(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)))
random.seed(1234)
def step(x, stepsize = 0.1):
return x + random.normal(0, stepsize)
plot_optim(greedy, 0.8, step, 2000, E, lambda x:1+0*x, [0.01, 0.5, 5.])
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
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
plot_optim(boltzmann_anneal, x0, step, nstep, E, schedule, [0.01, 0.5, 5.])
-c:3: RuntimeWarning: divide by zero encountered in true_divide
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
plot_optim(boltzmann_anneal, x0, step, nstep, E, schedule, [0.1, 1., 10.])
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
plot_optim(boltzmann_anneal, x0, step, nstep, E, schedule, [0.1, 1., 10.])
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
plot_optim(boltzmann_anneal, x0, step, nstep, E, schedule, [0.1, 1., 10., 100., 1., 1., 1.])
-c:9: RuntimeWarning: overflow encountered in exp