Markov Chain Monte Carlo

link管理

链接快照平台

输入网页链接，自动生成快照
标签化管理网页链接

相关文章推荐

高大的豌豆 · 为 Amazon EKS ...· 2 周前 ·

打篮球的铅笔 · 原子能机构：将持续对日本福岛核电站处理水排放 ...· 1 月前 ·

有情有义的菠萝 · Oracle获取一周前，一个月前，一年前, ...· 1 月前 ·

重感情的丝瓜 · Mac Axure8 报错，无法打开 - ...· 2 月前 ·

淡定的啄木鸟 · 操作方法：在多台计算机上自动安装和升级 ...· 2 月前 ·

VIDEO_TAG = """<video controls> <source src="data:video/x-m4v;base64,{0}" type="video/mp4"> Your browser does not support the video tag. </video>""" def anim_to_html ( anim ): if not hasattr ( anim , '_encoded_video' ): with NamedTemporaryFile ( suffix = '.mp4' ) as f : anim . save ( f . name , fps = 20 , extra_args = [ '-vcodec' , 'libx264' ]) video = open ( f . name , "rb" ) . read () anim . _encoded_video = video . encode ( "base64" ) return VIDEO_TAG . format ( anim . _encoded_video )

We start from a simple task, compute the expectation

$$ {\rm I\kern-.3em E}[\phi(x)] , $$

In which

$$ \phi(x) = 0.03 x + 0.2 ,$$$$ x \sim \frac{2}{3}(\mathcal{N}(2, 1) + \frac{1}{2} \mathcal{N}(-1, 1)) . $$

The true value should be 0.23 (as the expected output of target distribution is the 1)

def target_distribution(x):
    return (scipy.stats.norm.pdf(x, loc=2) + scipy.stats.norm.pdf(x, loc=-1) * 0.5) / 1.5
def phi(x):
    return 0.03 * x + 0.2
x_axis = np.arange(-5, 5, 0.01)
target_plot = plt.plot(x_axis, map(target_distribution, x_axis), label="Target distribution")
plt.plot(x_axis, map(phi, x_axis), label="phi(x)", c="r")
plt.legend()
plt.ylim(-0.1, 0.5)
plt.xlim(-6, 6)
plt.title("The plot of target distribution and phi function")
plt.show()
def sampling_distribution(x):
    return scipy.stats.norm.pdf(x, loc=2, scale=1.5)
fig = plt.figure()
ax1 = plt.subplot(211, xlim=(-6, 6), ylim=(-0.1, 0.5))
ax1.plot(x_axis, map(target_distribution, x_axis), label="Target distribution")
ax1.plot(x_axis, map(sampling_distribution, x_axis), label="Sampling distribution")
ax1.set_title("Samples drawed from a wrong normal distribution")
ax1_scat = ax1.scatter([], [], alpha=0.5, label="Samples")
ax1_scat_current = ax1.scatter([], [], alpha=0.5, c="k", s=100, marker="o")
ax1.legend()
ax2 = plt.subplot(212, xlim=(-1, 100), ylim=(0.1, 0.3))
ax2.set_title("The approximation of expectation")
ax2_plot, = ax2.plot([], [], label="Expectation of phi")
ax2.plot([-10, 300], [0.23, 0.23], 'k-', alpha=0.5, label="True value")
ax2.legend()
samples = []
expect_history = []
def animate(i):
    # Sample from normal distribution
    samples.append(np.random.normal(2, 0.8))
    ax1_scat.set_offsets( np.array(zip(samples, np.repeat(-0.05, len(samples)))) )
    ax1_scat_current.set_offsets([(samples[-1], -0.05)])
    # Compute expectation
    expect = sum([phi(x) for x in samples]) / len(samples)
    expect_history.append((i, expect))
    ax2_plot.set_data(map(np.array, zip(*expect_history)))
    return
anim = animation.FuncAnimation(fig, animate, frames=100, interval=20, blit=True)
display_animation(anim)
Rejection sampling works in the following way:
We make a sampling distribution $k \tilde Q(x)$, which is not normalized
Our target distribution $\tilde P(x)$ should be underneath  $k \tilde Q(x)$
Draw a sample $x \sim Q(x)$
Draw a random height $u \sim \mbox{uniform}[0, k \tilde Q(x)]$
We accept the sample if $u \lt \tilde P(x)$
np.random.seed(3)
# Draw the stational things in the animation
def sampling_distribution(x):
    return scipy.stats.norm.pdf(x, loc=0, scale=2) * 3
fig = plt.figure()
ax1 = plt.subplot(211, xlim=(-6, 6), ylim=(-0.1, 0.5))
ax1.plot(x_axis, map(target_distribution, x_axis), label="Target distribution")
ax1.plot(x_axis, map(sampling_distribution, x_axis), label="Sampling distribution")
ax1.set_title("An example of rejection sampling")
ax1_scat = ax1.scatter([], [], alpha=0.5, label="Accepted samples")
ax1_discard_scat = ax1.scatter([], [], alpha=0.2, label="Discarded samples", c="r")
ax1_scat_current = ax1.scatter([], [], alpha=0.5, c="k", s=100, marker="o")
ax1.legend()
ax2 = plt.subplot(212, xlim=(-1, 300), ylim=(0.1, 0.3))
ax2.set_title("The approximation of expectation")
ax2_plot, = ax2.plot([], [], label="Expectation of phi")
ax2.plot([-10, 300], [0.23, 0.23], 'k-', alpha=0.5, label="True value")
ax2.legend()
samples = []
discard_samples = []
expect_history = []
def animate(i):
    # Draw a sample
    sample = np.random.normal(0, 2)
    u = np.random.uniform(0, sampling_distribution(sample))
    discard_sample = u > target_distribution(sample)
    if discard_sample:
        discard_samples.append(sample)
    else:
        samples.append(sample)
    if samples:
        ax1_scat.set_offsets( np.array(zip(samples, np.repeat(-0.05, len(samples)))) )
    if discard_samples:
        ax1_discard_scat.set_offsets( np.array(zip(discard_samples, np.repeat(-0.05, len(discard_samples)))) )
    # Compute expectation
    if samples:
        expect = sum([phi(x) for x in samples]) / len(samples)
        expect_history.append((i, expect))
        ax2_plot.set_data(map(np.array, zip(*expect_history)))
    ax1_scat_current.set_offsets([(sample, -0.05)])
    return
anim = animation.FuncAnimation(fig, animate, frames=300, interval=20, blit=True)
display_animation(anim)
np.random.seed(3)
# Draw the stational things in the animation
def sampling_distribution(x):
    return scipy.stats.norm.pdf(x, loc=1, scale=3)
fig = plt.figure()
ax1 = plt.subplot(211, xlim=(-6, 6), ylim=(-0.1, 0.5))
ax1.plot(x_axis, map(target_distribution, x_axis), label="Target distribution")
ax1.plot(x_axis, map(sampling_distribution, x_axis), label="Sampling distribution")
ax1.set_title("An example of importance sampling")
ax1_scat = ax1.scatter([], [], alpha=0.5, label="Samples")
ax1_scat_current = ax1.scatter([], [], alpha=0.5, c="k", s=100, marker="o")
ax1.legend()
ax2 = plt.subplot(212, xlim=(-1, 300), ylim=(0.1, 0.3))
ax2.set_title("The approximation of expectation")
ax2_plot, = ax2.plot([], [], label="Expectation of phi")
ax2.plot([-10, 300], [0.23, 0.23], 'k-', alpha=0.5, label="True value")
ax2.legend()
samples = []
expect_history = []
def animate(i):
    # Draw a sample
    sample = np.random.normal(1, 3)
    samples.append(sample)
    ax1_scat.set_offsets( np.array(zip(samples, np.repeat(-0.05, len(samples)))) )
    ax1_scat_current.set_offsets([(sample, -0.05)])
    # Compute expectation
    expect = sum([phi(x) * target_distribution(x) / sampling_distribution(x) for x in samples]) / len(samples)
    expect_history.append((i, expect))
    ax2_plot.set_data(map(np.array, zip(*expect_history)))
    return ax1_scat
anim = animation.FuncAnimation(fig, animate, frames=300, interval=20, blit=True)
display_animation(anim)
In Metropolis-Hastings, we propose a trasitional probability distribution $Q(x'; x)$, say $\mathcal N(x, \sigma^2)$.




    

Given previous accepted sample, we draw the next sample from $Q(x'; x)$.
Then we accept it with probability $\min (1, \frac{P(x')Q(x; x')}{P(x)Q(x'; x)})$.
np.random.seed(3)
# Draw the stational things in the animation
def sampling_distribution(x, loc=0):
    return scipy.stats.norm.pdf(x, loc=loc, scale=0.5)
fig = plt.figure()
ax1 = plt.subplot(211, xlim=(-6, 6), ylim=(-0.1, 0.5))
ax1.plot(x_axis, map(target_distribution, x_axis), label="Target distribution")
cond_plot, = ax1.plot([], [], label="Transitional distribution", lw=1)
ax1.set_title("An example of Metropolis-Hastings sampling")
ax1_scat = ax1.scatter([], [], alpha=0.5, label="Accepted samples")
ax1_discard_scat = ax1.scatter([], [], alpha=0.2, label="Discarded samples", c="r")
ax1_scat_current = ax1.scatter([], [], alpha=0.5, c="k", s=100, marker="o")
ax1.legend()
ax2 = plt.subplot(212, xlim=(-1, 300), ylim=(0.1, 0.3))
ax2.set_title("The approximation of expectation")
ax2_plot, = ax2.plot([], [], label="Expectation of phi")
ax2.plot([-10, 300], [0.23, 0.23], 'k-', alpha=0.5, label="True value")
ax2.legend()
start_point = 0.
samples = [start_point]
discard_samples = []
expect_history = []
def animate(i):
    # Draw a sample
    sample = samples[-1] + np.random.normal(0, 0.5)
    cond_plot.set_data([x_axis, np.array([sampling_distribution(x, sample) for x in x_axis])])
    accept_probability = min([1., target_distribution(sample) / target_distribution(samples[-1])])
    discard_sample = np.random.uniform(0, 1) > accept_probability
    if discard_sample:
        discard_samples.append(sample)
    else:
        samples.append(sample)
    if samples:
        ax1_scat.set_offsets( np.array(zip(samples, np.repeat(-0.05, len(samples)))) )
    if discard_samples:
        ax1_discard_scat.set_offsets( np.array(zip(discard_samples, np.repeat(-0.05, len(discard_samples)))) )
    ax1_scat_current.set_offsets([(sample, -0.05)])
    # Compute expectation
    if samples:
        expect = sum([phi(x) for x in samples]) / len(samples)
        expect_history.append((i, expect))
        ax2_plot.set_data(map(np.array, zip(*expect_history)))
    return ax1_scat, cond_plot
anim = animation.FuncAnimation(fig, animate, frames=300, interval=20, blit=True)
display_animation(anim)
Gibbs sampling is a technique for sampling a joint distribution. To use it, we MUST know the true transitional probability distribution $P(x_i|\mathrm x_{j \neq i})$ of the target tribution (but not to propose one).
Algorithm:
Intitialize $\mathrm x$
Sample $x_i \sim P(x_i|\mathrm x_{j \neq i})$ in turn
Gibbs sampling does not produce rejections. But in the begining the samples can not be trust, so we just throw them away (burn-in).
ax = plt.axes(xlim=(-2, 4), ylim=(-2, 4))
scat = ax.scatter([], [], c="b", s=40, alpha=0.3)
scat2 = ax.scatter([], [], c="k", s=100, marker="o", alpha=0.5)
line_plot, = ax.plot([], [])
rho = 0.9
x, y = np.mgrid[-4:4:.01, -4:4:.01]
pos = np.dstack((x, y))
rv = multivariate_normal([1, 0], [[1., rho], [rho, 1.]])
ax.contourf(x, y, rv.pdf(pos), 30, alpha=0.3)
m1 = -1.5
m2 = 3.5
samples = [(m1, m2)]
def animate(i):
    global m1, m2
    new_m1 = np.random.normal(1 + rho * (m2 - 0), (1 - rho**2))
    new_m2 = np.random.normal(rho * (new_m1 - 1), (1 - rho**2))
    m1 = new_m1
    m2 = new_m2
    samples.append((new_m1, new_m2))
    samples_array = np.array(samples)
    scat.set_offsets(samples_array)
    scat2.set_offsets(samples_array[-1:])
    line_plot.set_data(zip(*samples[-2:]))
    return scat, scat2
anim = animation.FuncAnimation(fig, animate, frames=300, interval=20, blit=True)
display_animation(anim)