In [1]:

    

from neatplots.predefined import four as palette
from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec

%pylab inline
%load_ext autoreload
%autoreload 2

gray = (0.3, 0.3, 0.3)
seed(64321)

import latexstyle
latexstyle.setup()


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

def draw_from_proposal_dist(x0, area, proposal_std):
    x = np.inf
    while x < area[0] or x > area[1]:
        x = x0 + proposal_std * randn()
    return x

def sample_with_metropolis_hastings(
        density, x0, area, num_samples, proposal_std):
    positions = np.empty(num_samples)
    values = np.empty(num_samples)

    f = 0
    while f <= 0:
        x = draw_from_proposal_dist(x0, area, proposal_std)
        f = density(x)

    for i in xrange(num_samples):
        x_new = draw_from_proposal_dist(x, area, proposal_std)
        f_new = density(x_new)
        if f <= 0:
            acceptance_ratio = 1
        else:
            acceptance_ratio = f_new / f
        if rand() < acceptance_ratio:
            x = x_new
            f = f_new

        positions[i] = x_new
        values[i] = f

    return positions, values


def gen_probe_locations(area, density):
    area = np.asarray(area)
    source = 0.0

    num_uniform_samples = 0
    num_samples_per_source = 50
    mh_stride = 5

    uniform_samples = (area[1] - area[0]) * rand(num_uniform_samples) + area[0]

    samples = sample_with_metropolis_hastings(
        density, source, area, num_samples_per_source, 2)[0][::mh_stride]

    samples_gauss = []
    for i in xrange(mh_stride):
        samples_gauss.extend(
            2 * randn() + s for s in samples)
    samples_gauss = filter(
        lambda x: np.all(x > area[0]) and np.all(x < area[1]), samples_gauss)

    return np.concatenate([uniform_samples, samples, samples_gauss])

def gen_probe_locations_qrsim(area, density, num):
    samples = []
    while len(samples) < num:
        num_missing = num - len(samples)
        new_samples = (area[1] - area[0]) * rand(num_missing) + area[0]
        samples.extend(filter(lambda x: density(x) > 1e-3, new_samples))
    return samples


    

In [3]:

    

density = lambda x: np.exp(-(x**2) / 2.0)
#def density(x):
#    x = np.atleast_1d(x)
#    mask = x > 0
#    a = 2 * 0.33 * (x[mask] ** 0.86)
#    y = np.zeros_like(x)
#    y[mask] = 1 / np.pi / 6 / a * (np.exp(-(5.0 ** 2) / a) + np.exp(-(45.0 ** 2) / a))
#    return y

x_range = (-7, 13)
xs = np.linspace(x_range[0], x_range[1], 150)
ys = density(xs)


    

In [4]:

    

fig = plt.figure(figsize=(3, 1.5))
ax = fig.add_subplot(1, 1, 1)

loc = gen_probe_locations(x_range, density)
print len(loc)
loc_qrsim = gen_probe_locations_qrsim(x_range, density, len(loc))

ax.plot(xs, ys, c=palette.thin[1])
ax.scatter(loc_qrsim, np.zeros_like(loc_qrsim) - 0.1, c=palette.thick[0], marker='s', label='QRSim sampling')
ax.scatter(loc, np.zeros_like(loc) - 0.2, c=palette.thick[3], label='MH based sampling')

latexstyle.style_axes(ax)
ax.set_xlim(*x_range)
ax.set_ylim(-0.3, 1.0)
ax.legend(frameon=False, columnspacing=1.5, handletextpad=0.2)

ax.set_xlabel(r'$x$', labelpad=0)
ax.set_ylabel(r'$y$', labelpad=0, rotation='horizontal', verticalalignment='center')


    

    




59






    
Out[4]:





<matplotlib.text.Text at 0x111792090>






    

In [5]:

    

fig.savefig('../../thesis/plots/err-sampling.pdf')


    

In [5]: