Interactive Plotting

In [1]:
using Gadfly,Interact,Compose
In [3]:
using Reactive
In [2]:
using Distributions

Interactive Plotting for the Gibbs sampler

simulate data

In [3]:
##These parameter can be modified.
niter=nrow=200
m=[0,0]
v=[1.0 0.6
   0.6 1.0];
In [4]:
ncol=2
ypair=Array(Float64,nrow,ncol);

y=Array(Float64,1,2)

s12 =sqrt( v[1,1] - v[1,2]*v[1,2]/v[2,2]);
s21 = sqrt(v[2,2] - v[1,2]*v[1,2]/v[1,1]);
for (iter in 1:niter)
    m12 = m[1] + v[1,2]/v[2,2]*(y[2] - m[2]);
    y[1] = rand(Normal(m12,s12));
    m21 = m[2] + v[1,2]/v[1,1]*(y[1] - m[1]);
    y[2] = rand(Normal(m21,s21));    
    ypair[iter,:]=y
end


lower_1 = m[1]-3*v[1,1]
upper_1 = m[1]+3*v[1,1]
lower_2 = m[2]-3*v[2,2]
upper_2 = m[2]+3*v[2,2];

##Interact
gady1,gady2=[ypair[1,1]],[ypair[1,2]]
    
for i=2:niter
    push!(gady1,ypair[i-1,1])
    push!(gady2,ypair[i,2])
    push!(gady1,ypair[i,1])
    push!(gady2,ypair[i,2])
end
In [5]:
set_default_plot_size(20cm, 16cm)

@manipulate for n=1:niter        
    Gadfly.plot(
    y=gady1[1:(n-2)], x=gady2[1:(n-2)],
    Geom.point, 
    Guide.ylabel("Trait 2"),
    Guide.xlabel("Trait 1"),
    Guide.title("Bivariate Samples - Gibbs Sampling"),
    Theme(default_color=colorant"orange",default_point_size=4pt),
    Scale.x_continuous(minvalue=-5, maxvalue=5),
    Scale.y_continuous(minvalue=-5, maxvalue=5),

    layer(y=[gady1[n-1]],x=[gady2[n-1]],
          Geom.point,
          Theme(default_color=colorant"blue",default_point_size=4pt)),

    layer(y=[gady1[n]],x=[gady2[n]],
          Geom.point,
          Theme(default_color=colorant"red",default_point_size=4pt))
    )
end
Out[5]:
Trait 1 -20 -15 -10 -5 0 5 10 15 20 -15.0 -14.5 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 -20 -10 0 10 20 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 -20 -15 -10 -5 0 5 10 15 20 -15.0 -14.5 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 -20 -10 0 10 20 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Trait 2 Bivariate Samples - Gibbs Sampling

Interactive Plotting for the Metropolis–Hastings algorithm

In [6]:
#Paramaters below can be modified.
m=[0.0 0.0]
v=[1.0 0.6
   0.6 1.0]

niter = 500;
In [7]:
nrow  = niter
ncol  = 2
ypair = Array(Float64,nrow,ncol);

y = zeros(1,2)
ynew = Array(Float64,1,2)
ycount = 0
xcount = 0
all_accepted_pairs=Array(Float64,nrow,2);
all_rejected_pairs=Array(Float64,nrow,2);
vi = inv(v)


m1 = 0;
m2 = 0;
xx = 0;
y1 = 0;
delta = 1.0;
min1 = -delta*sqrt(v[1,1]);
max1 = +delta*sqrt(v[1,1]);
min2 = -delta*sqrt(v[2,2]);
max2 = +delta*sqrt(v[2,2]);


lower_1 = m[1]-3*sqrt(v[1,1])
upper_1 = m[1]+3*sqrt(v[1,1])
lower_2 = m[2]-3*sqrt(v[2,2])
upper_2 = m[2]+3*sqrt(v[2,2])


ytmp=Array(Float64,nrow,4)
xtmp=Array(Float64,nrow,4)
naccept=0
nreject=0

z = (y-m)'
denOld = exp(-0.5*z'*vi*z);
In [8]:
nrow=niter
recnum=Array(Int64,nrow,2)
points=Any[]
v3=zeros(3) # a vector as accepted or not / point x axis / point y axis

v3copy=copy(v3)
push!(points,v3copy)

for i=1:niter
    ynew[1] = y[1] + rand(Uniform(min1, max1))
    ynew[2] = y[2] + rand(Uniform(min2, max2))
    denNew  = exp(-0.5*(ynew-m)*vi*(ynew-m)');
    alpha   = denNew/denOld;    
    
    #create proposal interval
    xtmp[i,1]=y[1]+min1
    ytmp[i,1]=y[2]+min2
    xtmp[i,2]=y[1]+max1
    ytmp[i,2]=y[2]+min2    
    xtmp[i,3]=y[1]+min1
    ytmp[i,3]=y[2]+max2    
    xtmp[i,4]=y[1]+max1
    ytmp[i,4]=y[2]+max2
    
    #acceptance probability
    if rand()<alpha[1,1]   #Accept the new pair
       y[1] = ynew[1]
       y[2] = ynew[2] 
       denOld = exp(-0.5*(y-m)*vi*(y-m)'); 
       #naccept = naccept +1
       v3[1] = 1.0 
       #all_accepted_pairs[naccept,:]= y  
    else                  #Reject the new pair
       #nreject = nreject +1
       v3[1] = 0.0
       #all_rejected_pairs[nreject,:]= ynew   
    end
    
    v3[2]=ynew[1]
    v3[3]=ynew[2]
    v3copy=copy(v3)
    push!(points,v3copy)
end

numaccept=int(sum(points)[1])
numreject=niter+1-numaccept
all_accepted_pairs=Array(Float64,numaccept,2);
all_rejected_pairs=Array(Float64,numreject,2);
nreject=0
naccept=0
vecAccept=Int64[]

for i=1:niter
    if points[i][1]==0.0
       nreject += 1 
       all_rejected_pairs[nreject,:]= points[i][2:3]   
    else
       naccept += 1 
       all_accepted_pairs[naccept,:]= points[i][2:3]   
    end
    push!(vecAccept,naccept)
    
end
In [9]:
set_default_plot_size(20cm, 16cm)
In [10]:
@manipulate for n=1:(niter-1)        
    Gadfly.plot(
    x=[points[n+1][2]],y=[points[n+1][3]], #proposal
    Geom.point,
    Scale.x_continuous(minvalue=-3, maxvalue=3),
    Scale.y_continuous(minvalue=-3, maxvalue=3),
    Theme(default_color=colorant"red",default_point_size=4pt),
    Guide.ylabel("Trait 2"),
    Guide.xlabel("Trait 1"),
    Guide.title("Bivariate Samples - MH Sampling"),
    
    #current accepted point
    layer(x=[all_accepted_pairs[vecAccept[n],1]],y=[all_accepted_pairs[vecAccept[n],2]],
          Geom.point,
    Theme(default_color=colorant"black",default_point_size=4pt)),
 
    #poroposal distribution
    layer(y=ytmp[n,1:2],x=xtmp[n,1:2],
          Geom.line,
          Theme(default_color=colorant"blue",default_point_size=20pt)),    
    layer(y=ytmp[n,3:4],x=xtmp[n,3:4],
          Geom.line,
          Theme(default_color=colorant"blue",default_point_size=20pt)),
    layer(y=ytmp[n,[1,3]],x=xtmp[n,[1,3]],
          Geom.line,
          Theme(default_color=colorant"blue",default_point_size=20pt)),
    layer(y=ytmp[n,[2,4]],x=xtmp[n,[2,4]],
          Geom.line,
          Theme(default_color=colorant"blue",default_point_size=20pt)),

    
    #layer(x=[all_accepted_pairs[1:vecAccept[n],1]],y=[all_accepted_pairs[1:vecAccept[n],2]],
    #     Geom.point,
    #Theme(default_color=color("black"),default_point_size=4pt)),
    
    #layer(x=[all_rejected_pairs[1:(n-vecAccept[n]),1]],y=[all_rejected_pairs[1:(n-vecAccept[n]),2]],
    #      Geom.point,
    #Theme(default_color=color("black"),default_point_size=4pt)),
    
    Guide.annotation(
       compose(context(), circle([all_rejected_pairs[1:(n-vecAccept[n-1]),1]], [all_rejected_pairs[1:(n-vecAccept[n-1]),2]],
       [1.0mm]), fill(nothing),
       stroke(colorant"red"))),
    
    Guide.annotation(
       compose(context(), circle([all_accepted_pairs[1:vecAccept[n-1],1]], [all_accepted_pairs[1:vecAccept[n-1],2]],
       [1.0mm]), fill(nothing),
       stroke("black")))
    )
end
Out[10]:
Trait 1 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 -9.0 -8.8 -8.6 -8.4 -8.2 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 -10 -5 0 5 10 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 -20 -10 0 10 20 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 Trait 2 Bivariate Samples - MH Sampling
In [39]:
Gadfly.plot(sin, 0, 2pi,
     Guide.annotation(
      compose(context(), circle([pi/2, 3*pi/2], [1.0, -1.0], [2mm]), fill(nothing),
       stroke(colorant"orange"))))
compose not defined
while loading In[39], in expression starting on line 1

Appendix

Interact.jl works nicely with Jupyter notebook. To make interactive plotting in REPL (terminal), Winston is a better choice.

Winston

In [32]:
using Winston
Warning: using Winston.text in module Main conflicts with an existing identifier.
In [ ]:
p = FramedPlot(
         title="Bivariate Samples - Gibbs Sampling",
         xlabel="Trait 1",
         ylabel="Trait 2",
         xrange=(lower_1,upper_1),
         yrange=(lower_2,upper_2))

point1=Points(ypair[1,1],ypair[1,2],color="green")
point2=Points(ypair[1,1],ypair[1,2])
add(p,point1)

for iter=2:nrow
    readline(STDIN)
    point1=Points(ypair[iter-1,1],ypair[iter,2],color="green")
    add(p,point1,point2)
    display(p)
    point1=Points(ypair[iter-1,1],ypair[iter,2]) #change to black

    readline(STDIN)
    point2=Points(ypair[iter,1],ypair[iter,2],color="green")
    add(p, point1,point2)
    display(p)
    point2=Points(ypair[iter,1],ypair[iter,2]) #change to black
end;