Échantillonnage de tranche

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from numpy.random import rand, seed
 
def p_hat(x):
    pi1, pi2 = 0.4, 0.6
    return pi1 * np.exp(-(x + 1)**2) + pi2 * np.exp(-(x - 5) **2  / 2)
 
xrange = np.linspace(-5, 8, 500)
plt.plot(xrange, p_hat(xrange))
plt.title(r"$\tilde p(x)$", fontsize=15);
 
plt.show()
 
z = -2
pz = p_hat(z)
 
print(pz)
 
u = rand() * pz
w = 3
inc = 2
zmin, zmax = z - w, z + w
while p_hat(zmin) > u:
    zmin -= inc
 
while p_hat(zmax) > u:
    zmax += inc
 
    xrange = np.linspace(-5, 8, 500)
plt.plot(xrange, p_hat(xrange))
plt.axhline(y=u, c="tab:gray", linestyle="--", alpha=0.4)
plt.hlines(u, zmin, zmax)
plt.title(r"$\tilde p(x)$", fontsize=15);
 
plt.show()
 
seed(3141)
z_prime = rand() * (zmax - zmin) + zmin
 
if p_hat(z_prime) > u:
    z = z_prime
else:
    zmin_dist = abs(z_prime - zmin)
    zmax_dist = abs(z_prime - zmax)
    if  zmin_dist < zmax_dist:
        zmin = z_prime
    else:
        zmax = z_prime
 
xrange = np.linspace(-5, 8, 500)
plt.plot(xrange, p_hat(xrange))
plt.axhline(y=u, c="tab:gray", linestyle="--", alpha=0.4)
plt.hlines(u, zmin, zmax)
plt.axvline(x=z_prime, c="tab:red")
plt.title(r"$\tilde p(x)$", fontsize=15);
 
plt.show()
 
def one_sample(p_hat, zmin, zmax):        
    z_prime = rand() * (zmax - zmin) + zmin
    if p_hat(z_prime) > u:
        return z_prime
    else:
        zmin_dist = abs(z_prime - zmin)
        zmax_dist = abs(z_prime - zmax)
        if  zmin_dist < zmax_dist:
            zmin = z_prime
        else:
            zmax = z_prime
 
        return one_sample(p_hat, zmin, zmax)
 
def slice_sampling(p_hat, z_init, w, inc=0.2):
    z = z_init
    pz = p_hat(z)
    u = rand() * pz
 
    zmin, zmax = z - w, z + w
    while p_hat(zmin) > u:
        zmin -= inc
 
    while p_hat(zmax) > u:
        zmax += inc
 
    znew = one_sample(p_hat, zmin, zmax)
    return znew
 
z_vec = []
for _ in range(20_000):
    z = slice_sampling(p_hat, z_init=4, w=5, inc=2)
    z_vec.append(z)
 
fig, ax = plt.subplots()
ax2 = ax.twinx()
sns.distplot(z_vec, ax=ax)
ax2.plot(xrange, p_hat(xrange), c="tab:orange")
ax.set_ylabel("Samples", fontsize=13);
ax2.set_ylabel("Ground Truth", fontsize=13);
 
plt.show()
 
z_vec = []
for _ in range(20_000):
    z = slice_sampling(p_hat, z_init=4, w=1, inc=2)
    z_vec.append(z)
 
fig, ax = plt.subplots()
ax2 = ax.twinx()
sns.distplot(z_vec, ax=ax)
ax2.plot(xrange, p_hat(xrange), c="tab:orange")
ax.set_ylabel("Samples", fontsize=13);
ax2.set_ylabel("Ground Truth", fontsize=13);
 
plt.show()