{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Interactive Plotting"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using Gadfly,Interact,Compose"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using Distributions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interactive Plotting for the Gibbs sampler"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###simulate data"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"##These parameter can be modified.\n",
"niter=nrow=200\n",
"m=[0,0]\n",
"v=[1.0 0.6\n",
" 0.6 1.0];"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ncol=2\n",
"ypair=Array(Float64,nrow,ncol);\n",
"\n",
"y=Array(Float64,1,2)\n",
"\n",
"s12 =sqrt( v[1,1] - v[1,2]*v[1,2]/v[2,2]);\n",
"s21 = sqrt(v[2,2] - v[1,2]*v[1,2]/v[1,1]);\n",
"for (iter in 1:niter)\n",
" m12 = m[1] + v[1,2]/v[2,2]*(y[2] - m[2]);\n",
" y[1] = rand(Normal(m12,s12));\n",
" m21 = m[2] + v[1,2]/v[1,1]*(y[1] - m[1]);\n",
" y[2] = rand(Normal(m21,s21)); \n",
" ypair[iter,:]=y\n",
"end\n",
"\n",
"\n",
"lower_1 = m[1]-3*v[1,1]\n",
"upper_1 = m[1]+3*v[1,1]\n",
"lower_2 = m[2]-3*v[2,2]\n",
"upper_2 = m[2]+3*v[2,2];\n",
"\n",
"##Interact\n",
"gady1,gady2=[ypair[1,1]],[ypair[1,2]]\n",
" \n",
"for i=2:niter\n",
" push!(gady1,ypair[i-1,1])\n",
" push!(gady2,ypair[i,2])\n",
" push!(gady1,ypair[i,1])\n",
" push!(gady2,ypair[i,2])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [],
"text/plain": [
"Interact.Slider{Int64}([Reactive.Input{Int64}] 100,\"n\",100,1:200)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 15,
"metadata": {
"comm_id": "05ae4472-350c-45ae-84ca-9d09c4e10142",
"reactive": true
},
"output_type": "execute_result"
}
],
"source": [
"set_default_plot_size(20cm, 16cm)\n",
"\n",
"@manipulate for n=1:niter \n",
" Gadfly.plot(\n",
" y=gady1[1:(n-2)], x=gady2[1:(n-2)],\n",
" Geom.point, \n",
" Guide.ylabel(\"Trait 2\"),\n",
" Guide.xlabel(\"Trait 1\"),\n",
" Guide.title(\"Bivariate Samples - Gibbs Sampling\"),\n",
" Theme(default_color=colorant\"orange\",default_point_size=4pt),\n",
" Scale.x_continuous(minvalue=-5, maxvalue=5),\n",
" Scale.y_continuous(minvalue=-5, maxvalue=5),\n",
"\n",
" layer(y=[gady1[n-1]],x=[gady2[n-1]],\n",
" Geom.point,\n",
" Theme(default_color=colorant\"blue\",default_point_size=4pt)),\n",
"\n",
" layer(y=[gady1[n]],x=[gady2[n]],\n",
" Geom.point,\n",
" Theme(default_color=colorant\"red\",default_point_size=4pt))\n",
" )\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Interactive Plotting for the Metropolis–Hastings algorithm "
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#Paramaters below can be modified.\n",
"m=[0.0 0.0]\n",
"v=[1.0 0.6\n",
" 0.6 1.0]\n",
"\n",
"niter = 500;"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"nrow = niter\n",
"ncol = 2\n",
"ypair = Array(Float64,nrow,ncol);\n",
"\n",
"y = zeros(1,2)\n",
"ynew = Array(Float64,1,2)\n",
"ycount = 0\n",
"xcount = 0\n",
"all_accepted_pairs=Array(Float64,nrow,2);\n",
"all_rejected_pairs=Array(Float64,nrow,2);\n",
"vi = inv(v)\n",
"\n",
"\n",
"m1 = 0;\n",
"m2 = 0;\n",
"xx = 0;\n",
"y1 = 0;\n",
"delta = 1.0;\n",
"min1 = -delta*sqrt(v[1,1]);\n",
"max1 = +delta*sqrt(v[1,1]);\n",
"min2 = -delta*sqrt(v[2,2]);\n",
"max2 = +delta*sqrt(v[2,2]);\n",
"\n",
"\n",
"lower_1 = m[1]-3*sqrt(v[1,1])\n",
"upper_1 = m[1]+3*sqrt(v[1,1])\n",
"lower_2 = m[2]-3*sqrt(v[2,2])\n",
"upper_2 = m[2]+3*sqrt(v[2,2])\n",
"\n",
"\n",
"ytmp=Array(Float64,nrow,4)\n",
"xtmp=Array(Float64,nrow,4)\n",
"naccept=0\n",
"nreject=0\n",
"\n",
"z = (y-m)'\n",
"denOld = exp(-0.5*z'*vi*z); "
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nrow=niter\n",
"recnum=Array(Int64,nrow,2)\n",
"points=Any[]\n",
"v3=zeros(3) # a vector as accepted or not / point x axis / point y axis\n",
"\n",
"v3copy=copy(v3)\n",
"push!(points,v3copy)\n",
"\n",
"for i=1:niter\n",
" ynew[1] = y[1] + rand(Uniform(min1, max1))\n",
" ynew[2] = y[2] + rand(Uniform(min2, max2))\n",
" denNew = exp(-0.5*(ynew-m)*vi*(ynew-m)');\n",
" alpha = denNew/denOld; \n",
" \n",
" #create proposal interval\n",
" xtmp[i,1]=y[1]+min1\n",
" ytmp[i,1]=y[2]+min2\n",
" xtmp[i,2]=y[1]+max1\n",
" ytmp[i,2]=y[2]+min2 \n",
" xtmp[i,3]=y[1]+min1\n",
" ytmp[i,3]=y[2]+max2 \n",
" xtmp[i,4]=y[1]+max1\n",
" ytmp[i,4]=y[2]+max2\n",
" \n",
" #acceptance probability\n",
" if rand()\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 20,
"metadata": {
"comm_id": "2dadc49d-648d-4399-a316-a767248d2f42",
"reactive": true
},
"output_type": "execute_result"
}
],
"source": [
"@manipulate for n=1:(niter-1) \n",
" Gadfly.plot(\n",
" x=[points[n+1][2]],y=[points[n+1][3]], #proposal\n",
" Geom.point,\n",
" Scale.x_continuous(minvalue=-3, maxvalue=3),\n",
" Scale.y_continuous(minvalue=-3, maxvalue=3),\n",
" Theme(default_color=colorant\"red\",default_point_size=4pt),\n",
" Guide.ylabel(\"Trait 2\"),\n",
" Guide.xlabel(\"Trait 1\"),\n",
" Guide.title(\"Bivariate Samples - MH Sampling\"),\n",
" \n",
" #current accepted point\n",
" layer(x=[all_accepted_pairs[vecAccept[n],1]],y=[all_accepted_pairs[vecAccept[n],2]],\n",
" Geom.point,\n",
" Theme(default_color=colorant\"black\",default_point_size=4pt)),\n",
" \n",
" #poroposal distribution\n",
" layer(y=ytmp[n,1:2],x=xtmp[n,1:2],\n",
" Geom.line,\n",
" Theme(default_color=colorant\"blue\",default_point_size=20pt)), \n",
" layer(y=ytmp[n,3:4],x=xtmp[n,3:4],\n",
" Geom.line,\n",
" Theme(default_color=colorant\"blue\",default_point_size=20pt)),\n",
" layer(y=ytmp[n,[1,3]],x=xtmp[n,[1,3]],\n",
" Geom.line,\n",
" Theme(default_color=colorant\"blue\",default_point_size=20pt)),\n",
" layer(y=ytmp[n,[2,4]],x=xtmp[n,[2,4]],\n",
" Geom.line,\n",
" Theme(default_color=colorant\"blue\",default_point_size=20pt)),\n",
"\n",
" \n",
" #layer(x=[all_accepted_pairs[1:vecAccept[n],1]],y=[all_accepted_pairs[1:vecAccept[n],2]],\n",
" # Geom.point,\n",
" #Theme(default_color=color(\"black\"),default_point_size=4pt)),\n",
" \n",
" #layer(x=[all_rejected_pairs[1:(n-vecAccept[n]),1]],y=[all_rejected_pairs[1:(n-vecAccept[n]),2]],\n",
" # Geom.point,\n",
" #Theme(default_color=color(\"black\"),default_point_size=4pt)),\n",
" \n",
" Guide.annotation(\n",
" compose(context(), circle([all_rejected_pairs[1:(n-vecAccept[n-1]),1]], [all_rejected_pairs[1:(n-vecAccept[n-1]),2]],\n",
" [1.0mm]), fill(nothing),\n",
" stroke(colorant\"red\"))),\n",
" \n",
" Guide.annotation(\n",
" compose(context(), circle([all_accepted_pairs[1:vecAccept[n-1],1]], [all_accepted_pairs[1:vecAccept[n-1],2]],\n",
" [1.0mm]), fill(nothing),\n",
" stroke(\"black\")))\n",
" )\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.0",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.4.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}