In [ ]:
using PyPlot
srcpath = "/Users/ethananderes/Dropbox/BayesLense/src/"
savepath = "/Users/ethananderes/Google\ Drive/BayesLenseRev1/paper_rev2/"
seed = 2
srand(seed)
In [ ]:
rho = -0.99
sigma1 = [1 rho; rho 1]
sigmainv1 = inv(sigma1)
sigmainv2 = [1 0;-1 1]* sigmainv1 * [1 -1;0 1]
sigma2 = inv(sigmainv2)
pdfimag1(x) = (1/ (2 * pi)) * sqrt(det(sigmainv1)) * exp(-0.5 * (transpose(x) * sigmainv1 * x)[1])
pdfimag2(x) = (1/ (2 * pi)) * sqrt(det(sigmainv2)) * exp(-0.5 * (transpose(x) * sigmainv2 * x)[1])
xlim = 2.0
space = 0.01
gridpoints = Array{Float64,1}[[x,y] for x=(-xlim):space:xlim, y=(-xlim):space:xlim] #|> transpose |> flipud
imag1 = map(pdfimag1, gridpoints)
imag2 = map(pdfimag2, gridpoints)
sum(imag1 * ^(space,2))
function findContours(density, gridarea, probvec)
hts = Float64[]
prb = Float64[]
for dh in density[1:50:end]
push!(hts, dh)
push!(prb, sum( gridarea*density[density .> dh] ))
end
answr = Float64[]
for k in probvec
toss, index = findmin(abs(prb .- k))
push!(answr, hts[index])
end
answr
end
In [ ]:
thetachain1 = [0.0]
phichain1 = [0.0]
for k = 1:20
push!(thetachain1, sigma1[1,2] * phichain1[end] * sqrt(sigma1[1,1] / sigma1[2,2]) + randn() * sqrt(1-sigma1[1,2]^2) * sqrt(sigma1[1,1]) )
push!(phichain1, phichain1[end])
push!(phichain1, sigma1[1,2] * thetachain1[end] * sqrt(sigma1[2,2] / sigma1[1,1]) + randn() * sqrt(1-sigma1[1,2]^2) * sqrt(sigma1[2,2]) )
push!(thetachain1, thetachain1[end])
end
tildethetachain2 = [0.0]
phichain2 = [0.0]
for k = 1:20
push!(tildethetachain2, sigma2[1,2] * phichain2[end] * sqrt(sigma2[1,1] / sigma2[2,2]) + randn() * sqrt(1-sigma2[1,2]^2) * sqrt(sigma2[1,1]) )
push!(phichain2, phichain2[end])
push!(phichain2, sigma2[1,2] * tildethetachain2[end] * sqrt(sigma2[2,2] / sigma2[1,1]) + randn() * sqrt(1-sigma2[1,2]^2) * sqrt(sigma2[2,2]) )
push!(tildethetachain2, tildethetachain2[end])
end
probvec = [.5, .85]
cnt1 = findContours(imag1, space^2, probvec)
cnt2 = findContours(imag2, space^2, probvec)
In [ ]:
contour(imag1,
levels = cnt1,
origin= "lower",
extent=[-2.5,2.5,-2.5,2.5],
colors="k",
linewidth=1.5,
label=L"$P(t, \varphi|data)$"
)
plot(phichain1,thetachain1, linewidth=1.5)
ax = gca()
ax[:spines]["right"][:set_color]("none")
ax[:spines]["top"][:set_color]("none")
ax[:xaxis][:set_ticks_position]("bottom")
ax[:spines]["bottom"][:set_position](("data",0))
ax[:yaxis][:set_ticks_position]("left")
ax[:spines]["left"][:set_position](("data",0))
ax[:text](0.6, 0.6, L"$P$ $(t, \varphi|data)$", fontsize=27)
xticks([-2,-1,1,2], fontsize = 15)
yticks([-2,-1,1,2], fontsize = 15)
savefig(joinpath(savepath, "figure1a.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)
In [ ]:
contour(imag2, levels = cnt2, origin= "lower", extent=[-2.5,2.5,-2.5,2.5], colors="k", linewidth=1.5)
xlabel(L"$\varphi$", fontsize = 28),
ylabel(L"$\widetilde{t}$", fontsize = 28)
plot(phichain2, tildethetachain2, linewidth=1.5)
ax = gca()
ax[:spines]["right"][:set_color]("none")
ax[:spines]["top"][:set_color]("none")
ax[:xaxis][:set_ticks_position]("bottom")
ax[:spines]["bottom"][:set_position](("data",0))
ax[:yaxis][:set_ticks_position]("left")
ax[:spines]["left"][:set_position](("data",0))
ax[:text](0.7, 0.7, L"$P$ $(\widetilde{t}, \varphi|data)$", fontsize=27)
xticks([-2,-1,1,2], fontsize = 15)
yticks([-2,-1,1,2], fontsize = 15)
savefig(joinpath(savepath, "figure1b.pdf"), dpi=300, bbox_inches="tight", pad_inches=0)