from __future__ import division
import os
import sys
import glob
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline
%precision 4
plt.style.use('ggplot')
Most of the notebook is taken from http://twiecki.github.io/blog/2014/01/02/visualizing-mcmc/
! pip install -U git+git://github.com/jakevdp/JSAnimation &> /dev/null
You will aslo need to install ffmpeg
- on Ubunutu this can be done
with
sudo add-apt-repository ppa:jon-severinsson/ffmpeg
sudo apt-get update
sudo apt-get install ffmpeg
sudo apt-get install frei0r-plugins
import warnings
warnings.filterwarnings("ignore")
from JSAnimation import IPython_display
from matplotlib import animation
import pymc3 as pm
# Generate some data
np.random.seed(124)
size = 50
true_intercept = 1
true_slope = 2
x = np.linspace(0, 1, size)
y = true_intercept + x*true_slope + np.random.normal(scale=.5, size=size)
data = dict(x=x, y=y)
# Quickly hacked plotting code
samples = 600
fig = plt.figure(figsize=(6, 6))
i_width = (true_intercept-.7, true_intercept+.7)
s_width = (true_slope-.7, true_slope+.7)
samples_width = (0, samples)
ax1 = fig.add_subplot(221, xlim=i_width, ylim=samples_width)
ax2 = fig.add_subplot(224, xlim=samples_width, ylim=s_width)
ax3 = fig.add_subplot(223, xlim=i_width, ylim=s_width,
xlabel='intercept',
ylabel='slope')
fig.subplots_adjust(wspace=0.0, hspace=0.0)
line1, = ax1.plot([], [], lw=1)
line2, = ax2.plot([], [], lw=1)
line3, = ax3.plot([], [], 'o', lw=2, alpha=.1)
line4, = ax3.plot([], [], lw=1, alpha=.3)
line5, = ax3.plot([], [], 'k', lw=1)
line6, = ax3.plot([], [], 'k', lw=1)
ax1.set_xticklabels([])
ax2.set_yticklabels([])
#path = plt.scatter([], [])
lines = [line1, line2, line3, line4, line5, line6]
def init():
for line in lines:
line.set_data([], [])
return lines
def animate(i):
with model:
if i == samples * .75:
for j in range(500): iter_sample.next()
trace = iter_sample.next()
line1.set_data(trace['Intercept'][::-1], range(len(trace['Intercept'])))
line2.set_data(range(len(trace['x'])), trace['x'][::-1])
line3.set_data(trace['Intercept'], trace['x'])
line4.set_data(trace['Intercept'], trace['x'])
intercept = trace['Intercept'][-1]
x = trace['x'][-1]
line5.set_data([intercept, intercept], [x, s_width[1]])
line6.set_data([intercept, i_width[1]], [x, x])
return lines
with pm.Model() as model:
pm.glm.glm('y ~ x', data)
step = pm.Metropolis()
iter_sample = pm.iter_sample(samples+1000, step)
animation.FuncAnimation(fig, animate, init_func=init,
frames=samples, interval=5, blit=True)