hw2 for cogley's course

hw2_plotting.jl

# load in the model and associated tools
require("../nk_model.jl")
require("../nk_model_tools.jl")

# load packages we use
import Gadfly
import PyPlot
using Distributions
using DataFrames


# function plot_hpd(d::Distribution)
#     fig, ax = subplots()
#     x = 0.0:0.01:2.0
#     a, b = hpd_cs(d)
#     ax[:plot](x, pdf(d, x))
#     ax[:vlines](a, 0.0, 3.5)
#     ax[:vlines](b, 0.0, 3.5)
# end


function plot_post_sample(chains::Vector{MetropolisChain},
                          true_vals::ModelParams,
                          ti="Posterior sample from Metropolis algorithm")
    fig, axis = PyPlot.subplots(4, 3)
    axis = axis[:]  # flatten

    n_j = size(lots_of_chains, 1)  # number of lots_of_chains
    max_kde_j = min(n_j, 4) # only do kde for 4 lots_of_chains
    n_params = lots_of_chains[1].task.model.size  # number of parameters

    # set up colors for fill_between
    colors = PyPlot.plt.rcParams["axes.color_cycle"][1:max_kde_j]

    # extract latex names for axis title
    nms = tex_names(true_vals)
    true_vals = vec(true_vals)  # treat it as an array below

    monster_chain = sum(lots_of_chains)

    for i=1:n_params  # assumes each column is a samples
        # super ghetto, but I'm tired of wasting time on this...
        x_max = i in [2, 3, 6, 7, 8] ? 1.1 : 2.1
        x_max = i == 7 ? 4.1 : x_max

        ax = axis[i]
        ax[:set_title](nms[i])
        data_max = 0.0
        true_val = true_vals[i]

        for j=1:max_kde_j  # accross lots_of_chains
            k = KernelDensity.kde(lots_of_chains[j].samples[:, i])
            ax[:plot](k.x, k.density, alpha=0.4, lw=1., color=colors[j])
            ax[:fill_between](k.x, 0.0, k.density, alpha=0.05, color=colors[j])
            data_max = max(data_max, maximum(k.density))
        end
        ax[:vlines](true_val, 0, data_max, color="k", linestyle="--")
        ax[:set_xlim](0.001, x_max)

        big_k = KernelDensity.kde(monster_chain.samples[:, i])
        ax[:plot](big_k.x, big_k.density, lw=2, color="k", linestyle=":",
                 label="Overall")
        ax[:legend](loc=0, fontsize=7)
        ax[:set_ylim](bottom=0.0)

    end
    fig[:suptitle](ti * " ($(n_j) chains)")
end


function plot_post_mode(J::MultivariateNormal, nms::Vector, true_vals)
    data = rand(J, 10000)

    fig, axis = PyPlot.subplots(4, 3)
    axis = axis[:]  # flatten
    color = PyPlot.plt.rcParams["axes.color_cycle"][1]

    for i=1:size(data, 1)
        # super ghetto, but I'm tired of wasting time on this...
        x_max = i in [2, 3, 6, 7, 8] ? 1.1 : 2.1
        x_max = i == 7 ? 4.1 : x_max

        ax = axis[i]
        ax[:set_title](nms[i])
        data_max = 0.0
        true_val = true_vals[i]
        k = kde(data[i, :][:])
        ax[:plot](k.x, k.density, alpha=0.5, color=color)
        ax[:fill_between](k.x, 0.0, k.density, color=color, alpha=0.2)
        ax[:vlines](true_val, 0, maximum(k.density), color="k",
                    linestyle="--")
        ax[:set_xlim](0.001, x_max)
    end
    fig[:suptitle]("Normal approximation to Posterior Mode")
    return fig
end


function mpl_plot(p::Prior, θ::ModelParams=true_θ)
    fig, axis = PyPlot.subplots(4, 3)
    axis = axis[:]  # flatten

    N = 150

    xd(d::Beta) = linspace(0.001, 1.1, N)
    xd(d::InverseGamma) = linspace(0.001, 2.1, N)

    dist_name(d::InverseGamma) = "IG($(d.shape), $(d.scale))"
    dist_name(d::Beta) = "Beta" * "($(d.alpha), $(d.beta))"

    color = PyPlot.plt.rcParams["axes.color_cycle"][1]

    tex_θ_names = tex_names(θ)

    for (i, nm) in enumerate(names(p))
        true_val = getfield(θ, nm)
        dist = getfield(p, nm)
        std_dist = std(dist)

        x = xd(dist)
        if nm == :σ_i
            x = linspace(0.001, 4.1, N)
        end

        pdf_x = pdf(dist, x)
        axis[i][:plot](x, pdf_x)     

        axis[i][:vlines](true_val, 0, maximum(pdf_x), color="k",
                         linestyle="--")
        axis[i][:fill_between](x, 0.0, pdf_x, color=color, alpha=0.1)
        ti = tex_θ_names[i] * ": " * dist_name(dist)
        axis[i][:set_title](ti)
    end
    fig[:suptitle]("Prior Distributions")
    return fig
end

function gadfly_plot(p::Prior, θ::ModelParams=true_θ)
    N = 150

    df = DataFrame();
    x = linspace(0.001, 4.1, N)
    df[:x] = x
    for nm in names(the_prior)
        true_val = getfield(true_θ, nm)
        dist = getfield(the_prior, nm)

        pdf_x = pdf(dist, x);
        df[nm] = pdf_x
    end

    new_df = stack(df, [2:size(df, 2)])

    p = Gadfly.plot(new_df, x=:x, y=:value, color=:variable, Gadfly.Geom.line,
                    Gadfly.Guide.ylabel(""),
                    Gadfly.Guide.xlabel(""),
                    Gadfly.Guide.title("Prior Distributions"))

    return p
end

hw2_present.ipynb

nk_model.jl

nk_model_tools.jl