# 乱数のseed値を指定

srand(9869)

using Distributions
import PyPlot
plt = PyPlot
hist = plt.plt[:hist]
hist2D = plt.plt[:hist2d]

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

mixmodel(a,b) = MixtureModel(Normal, [(0.0,1.0), (b,1.0)], [1-a, a])

log_likelihood(a,b,sample) = loglikelihood(mixmodel(a,b),sample)

# Sample is generated by the model for the parameters a0 and b0
#
function plot_sample_and_posterior(a0,b0,N)
    println("a0 = $a0, b0 = $b0, N = $N")
    
    sample = rand(mixmodel(a0,b0),N)
    
    plt.figure(figsize=(5,3))
    hist(sample, bins=Int(floor(sqrt(N))))
    plt.title("Sample histogram")
    
    a = [linspace(0,1,100);]
    b = [linspace(-1,2*b0+1,100);]
    aa, bb = meshgrid(a,b)
    c = log_likelihood(a0,b0,sample)
    z = ((a,b) -> exp(log_likelihood(a,b,sample) - c)).(aa, bb)

    plt.figure()
    plt.surf(a,b,z, cmap=plt.ColorMap("rainbow"))
    plt.title("Posterior distribution")
    plt.xlabel("a")
    plt.ylabel("b")
    
    plt.figure()
    plt.axes(facecolor="black")
    plt.pcolormesh(a,b,z, cmap="rainbow")
    plt.title("Posterior distribution")
    plt.xlabel("a")
    plt.ylabel("b")
    plt.colorbar()
end

@time plot_sample_and_posterior(0.5, 3.0, 100);

@time plot_sample_and_posterior(0.5, 1.0, 100);

@time plot_sample_and_posterior(0.5, 0.5, 100);

@time plot_sample_and_posterior(0.5, 0.2, 100);

@time plot_sample_and_posterior(0.5, 0.1, 100);

@time plot_sample_and_posterior(0.5, 0.1, 1000);

@time plot_sample_and_posterior(0.5, 0.1, 10000);

@time plot_sample_and_posterior(0.5, 0.1, 100000);

@time plot_sample_and_posterior(0.5, 0.2, 100);

@time plot_sample_and_posterior(0.5, 0.2, 1000);

@time plot_sample_and_posterior(0.5, 0.2, 10000);

@time plot_sample_and_posterior(0.5, 0.5, 100);

@time plot_sample_and_posterior(0.5, 0.5, 1000);

@time plot_sample_and_posterior(0.5, 0.5, 10000);

@time plot_sample_and_posterior(0.5, 0.6, 1000);

@time plot_sample_and_posterior(0.5, 0.7, 1000);

@time plot_sample_and_posterior(0.5, 0.8, 1000);

@time plot_sample_and_posterior(0.5, 0.9, 1000);

@time plot_sample_and_posterior(0.5, 1.0, 1000);

@time plot_sample_and_posterior(0.5, 2.0, 1000);

@time plot_sample_and_posterior(0.5, 3.0, 1000);

@time plot_sample_and_posterior(0.5, 0.1, 10^6);