using PyPlot
using Distributions
using Optim
linspace(a, b, l) = range(a, b, length=l)
using Printf

using Base64
displayfile(mime, file; tag="img") = open(file) do f
    display("text/html", """<$tag src="data:$mime;base64,$(base64encode(f))"/>""")
end

displayfile("image/png", "Seven scientists.png")

Y = [-27.020, 3.570, 8.191, 9.808, 9.603, 9.945, 10.056]
n = length(Y);

f(nu, x) = pdf(TDist(nu),x)
x = collect(linspace(-10,10,400))

figure(figsize=(6,3.6))
for nu in [0.25,0.5,1.0,2.0,10.0]
    plot(x, f.(nu,x), label="\$\\nu = $nu\$", linewidth=0.8)
end
grid("on")
legend()
title("Probability density functions of t-distributions");

f(nu, x) = pdf(TDist(nu),x)
g(sigma, x) = pdf(Normal(0,sigma),x)

nu = 1
sigma = 1.5

x = collect(linspace(-10,10,400))
figure(figsize=(6,3.6))
plot(x, g.(sigma,x), label="Normal(0,\$\\sigma\$=$sigma)", linewidth=0.8, linestyle="-")
plot(x, f.(nu,x), label="TDist(\$\\nu\$=$nu)", linewidth=0.8, linestyle="--")
grid("on")
legend()
title("Probability density functions");

logpdf_TDist1(mu,nu,x) = logpdf(TDist(nu), x-mu)
logpdf_TDist2(mu,nu,sigma,x) = logpdf(TDist(nu), (x-mu)/sigma) - log(sigma)

negloglik_Normal(mu,sigma)    = sigma > zero(sigma) ? -sum(logpdf(Normal(mu,sigma),Y[i])   for i in 1:n) : Inf
negloglik_Cauchy(mu,sigma)    = sigma > zero(sigma) ? -sum(logpdf(Cauchy(mu,sigma),Y[i])   for i in 1:n) : Inf
negloglik_TDist1(mu,nu)       = nu > zero(nu)       ? -sum(logpdf_TDist1(mu,nu,Y[i])       for i in 1:n) : Inf
negloglik_TDist2(mu,nu,sigma) = 
             (nu > zero(nu) && sigma > zero(sigma)) ? -sum(logpdf_TDist2(mu,nu,sigma,Y[i]) for i in 1:n) : Inf

negloglik_Normal(x) = negloglik_Normal(x[1],x[2])
negloglik_Cauchy(x) = negloglik_Cauchy(x[1],x[2])
negloglik_TDist1(x) = negloglik_TDist1(x[1],x[2])
negloglik_TDist2(x) = negloglik_TDist2(x[1],x[2],x[3])

lik_Normal(mu,sigma)    = exp(-negloglik_Normal(mu,sigma))
lik_Cauchy(mu,sigma)    = exp(-negloglik_Cauchy(mu,sigma))
lik_TDist1(mu,nu)       = exp(-negloglik_TDist1(mu,nu))
lik_TDist2(mu,nu,sigma) = exp(-negloglik_TDist2(mu,nu,sigma))
;

fit_Normal = fit_mle(Normal,Y)

opt_Normal = optimize(negloglik_Normal, [10.0, 1.0])
@show opt_Normal

@show nmu    = opt_Normal.minimizer[1]
@show nsigma = opt_Normal.minimizer[2]
@show nmu_str = @sprintf("%.1f", nmu)
@show nsigma_str = @sprintf("%.0f", nsigma)

pred_Normal(x) = pdf(Normal(nmu, nsigma), x)
@show maxloglik_Normal = - negloglik_Normal(nmu, nsigma)
@show AIC_Normal = -2*maxloglik_Normal + 2*2

sleep(0.2)

mu = collect(linspace(-10,15,100))
sigma = collect(linspace(6,24,100))
L = lik_Normal.(mu',sigma)
cmap = "OrRd"
figure()
pcolormesh(mu, sigma, L, cmap=cmap, vmin=0.0)
colorbar(label="likelihood")
xlabel("\$\\mu\$")
ylabel("\$\\sigma\$")
grid("on", linestyle=":")
title("Model:  \$y \\sim \\mathrm{Normal}(\\mu,\\sigma)\$")
text(-9, 6.3, "Sample: $Y", fontsize=9)
scatter([nmu],[nsigma], label="MLE")
legend()

opt_Cauchy = optimize(negloglik_Cauchy, [10.0, 1.0])
@show opt_Cauchy

@show cmu    = opt_Cauchy.minimizer[1]
@show csigma = opt_Cauchy.minimizer[2]
@show cmu_str = @sprintf("%.1f", cmu)
@show csigma_str = @sprintf("%.2f", csigma)

pred_Cauchy(x) = pdf(Cauchy(cmu, csigma), x)
@show maxloglik_Cauchy = - negloglik_Cauchy(cmu, csigma)
@show AIC_Cauchy = -2*maxloglik_Cauchy + 2*2

sleep(0.2)

mu = collect(linspace(8.5,11,100))
sigma = collect(linspace(eps(),2,100))
L = lik_Cauchy.(mu',sigma)
cmap = "OrRd"
figure()
pcolormesh(mu, sigma, L, cmap=cmap, vmin=0.0)
colorbar(label="likelihood")
xlabel("\$\\mu\$")
ylabel("\$\\sigma\$")
grid("on", linestyle=":")
title("Model:  \$y \\sim \\mu + \\sigma\\mathrm{Cauchy}()\$")
text(8.6, 0.05, "Sample: $Y", fontsize=9)
scatter([cmu],[csigma], label="MLE")
legend()

opt_TDist1 = optimize(negloglik_TDist1, [10.0, 1.0])
@show opt_TDist1

@show t1mu = opt_TDist1.minimizer[1]
@show t1nu = opt_TDist1.minimizer[2]
@show t1mu_str = @sprintf("%.1f", t1mu)
@show t1nu_str = @sprintf("%.2f", t1nu)

pred_TDist1(x) = pdf(TDist(t1nu), x-t1mu)
@show maxloglik_TDist1 = - negloglik_TDist1(t1mu, t1nu)
@show AIC_TDist1 = -2*maxloglik_TDist1 + 2*2

sleep(0.2)

mu = collect(linspace(8.5,11,100))
nu = collect(linspace(eps(),2,100))
L = lik_TDist1.(mu',nu)
cmap = "OrRd"
figure()
pcolormesh(mu, nu, L, cmap=cmap, vmin=0.0)
colorbar(label="likelihood")
xlabel("\$\\mu\$")
ylabel("\$\\nu\$")
grid("on", linestyle=":")
title("Model:  \$y \\sim \\mu + \\mathrm{TDist}(\\nu)\$")
text(8.6, 0.05, "Sample: $Y", fontsize=9)
scatter([t1mu],[t1nu], label="MLE")
legend()

opt_TDist2 = optimize(negloglik_TDist2, [10.0, 1.0, 1.0])
@show opt_TDist2

@show t2mu = opt_TDist2.minimizer[1]
@show t2nu = opt_TDist2.minimizer[2]
@show t2sigma = opt_TDist2.minimizer[3]
@show t2mu_str = @sprintf("%.1f", t2mu)
@show t2nu_str = @sprintf("%.2f", t2nu)
@show t2sigma_str = @sprintf("%.2f", t2sigma)

pred_TDist2(x) = pdf(TDist(t2nu), (x-t2mu)/t2sigma)/t2sigma
@show maxloglik_TDist2 = - negloglik_TDist2(t2mu, t2nu, t2sigma)
@show AIC_TDist2 = -2*maxloglik_TDist2 + 2*3

mu = collect(linspace(8.5,11,100))
nu = collect(linspace(eps(),2,100))
L = lik_TDist2.(mu', nu, t2sigma)
cmap = "OrRd"
figure()
pcolormesh(mu, nu, L, cmap=cmap, vmin=0.0)
colorbar(label="likelihood")
xlabel("\$\\mu\$")
ylabel("\$\\nu\$")
grid("on", linestyle=":")
title("Model:  \$y \\sim \\mu + \\sigma\\;\\mathrm{TDist}(\\nu)\$")
text(8.6, 1.65, "Sample: $Y", fontsize=9)
text(8.6, 1.85, "\$\\sigma = $t2sigma_str\$")
scatter([t2mu],[t2nu], label="MLE")
legend()

sleep(0.2)

mu = collect(linspace(8.5,11,100))
sigma = collect(linspace(eps(),1,100))
L = lik_TDist2.(mu', t2nu, sigma)
cmap = "OrRd"
figure()
pcolormesh(mu, sigma, L, cmap=cmap, vmin=0.0)
colorbar(label="likelihood")
xlabel("\$\\mu\$")
ylabel("\$\\sigma\$")
grid("on", linestyle=":")
title("Model:  \$y \\sim \\mu + \\sigma\\;\\mathrm{TDist}(\\nu)\$")
text(8.6, 0.85, "Sample: $Y", fontsize=9)
text(8.6, 0.92, "\$\\nu = $t2nu_str\$")
scatter([t2mu],[t2sigma], label="MLE")
legend()

@show AIC_Normal
@show AIC_Cauchy
@show AIC_TDist1
@show AIC_TDist2;

figure(figsize=(6,3.6))

x = collect(linspace(8.5,11.5,1000))
plot(x, pred_Normal.(x), label="Normal(\$\\mu=$nmu_str\$, \$\\sigma=$nsigma_str\$)", lw=0.8, "-")
plot(x, pred_Cauchy.(x), label="Cauchy(\$\\mu=$cmu_str\$, \$\\sigma=$csigma_str\$)", lw=1.2, ls=":")
plot(x, pred_TDist1.(x), label="$t1mu_str + TDist(\$\\nu=$t1nu_str\$)", lw=0.8, ls="--")
plot(x, pred_TDist2.(x), label="$t2mu_str + $t2sigma_str TDist(\$\\nu=$t2nu_str\$)", lw=0.8, ls="-.")
xlabel("y")
grid("on", linestyle=":")

y = zero(Y)
scatter(Y[4:n], y[4:n], s=15, marker="o", alpha=0.5, color="red", label="Sample")

legend(fontsize=9)

title("PDFs of the predictive distributions");

figure(figsize=(6,3.6))

x = collect(linspace(2,9,1000))
plot(x, pred_Normal.(x), label="Normal(\$\\mu=$nmu_str\$, \$\\sigma=$nsigma_str\$)", lw=0.8, "-")
plot(x, pred_Cauchy.(x), label="Cauchy(\$\\mu=$cmu_str\$, \$\\sigma=$csigma_str\$)", lw=1.2, ls=":")
plot(x, pred_TDist1.(x), label="$t1mu_str + TDist(\$\\nu=$t1nu_str\$)", lw=0.8, ls="--")
plot(x, pred_TDist2.(x), label="$t2mu_str + $t2sigma_str TDist(\$\\nu=$t2nu_str\$)", lw=0.8, ls="-.")
xlabel("y")
grid("on", linestyle=":")

y = zero(Y)
scatter(Y[2:3], y[2:3], s=15, marker="o", alpha=0.5, color="red", label="Sample")

legend(fontsize=9)

title("PDFs of the predictive distributions");

figure(figsize=(6,3.6))

x = collect(linspace(-40,-20,1000))
plot(x, pred_Normal.(x), label="Normal(\$\\mu=$nmu_str\$, \$\\sigma=$nsigma_str\$)", lw=0.8, "-")
plot(x, pred_Cauchy.(x), label="Cauchy(\$\\mu=$cmu_str\$, \$\\sigma=$csigma_str\$)", lw=1.2, ls=":")
plot(x, pred_TDist1.(x), label="$t1mu_str + TDist(\$\\nu=$t1nu_str\$)", lw=0.8, ls="--")
plot(x, pred_TDist2.(x), label="$t2mu_str + $t2sigma_str TDist(\$\\nu=$t2nu_str\$)", lw=0.8, ls="-.")
xlabel("y")
grid("on", linestyle=":")

y = zero(Y)
scatter(Y[1:1], y[1:1], s=15, marker="o", alpha=0.5, color="red", label="Sample")

legend(fontsize=9)

title("PDFs of the predictive distributions");

figure(figsize=(6,3.6))

x = collect(linspace(-40,40,1000))
plot(x, pred_Normal.(x), label="Normal(\$\\mu=$nmu_str\$, \$\\sigma=$nsigma_str\$)", lw=0.8, "-")
#plot(x, pred_Cauchy.(x), label="Cauchy(\$\\mu=$cmu_str\$, \$\\sigma=$csigma_str\$)", lw=1.2, ls=":")
#plot(x, pred_TDist1.(x), label="$t1mu_str + TDist(\$\\nu=$t1nu_str\$)", lw=0.8, ls="--")
#plot(x, pred_TDist2.(x), label="$t2mu_str + $t2sigma_str TDist(\$\\nu=$t2nu_str\$)", lw=0.8, ls="-.")
xlabel("y")
grid("on", linestyle=":")

y = zero(Y)
y[5] = 0.001
y[6] = 0.002
y[7] = 0.003
scatter(Y, y, s=15, marker="o", alpha=0.5, color="red", label="Sample")

legend(fontsize=9)

title("PDFs of the predictive distribution");

@show csigma
@show t1nu
@show t2nu, t2sigma

MyInvGamma(nu, sigma) = InverseGamma(nu/2,nu*sigma^2/2)
InvGamma_Cauchy = MyInvGamma(1.0, csigma)
InvGamma_TDist1 = MyInvGamma(t1nu, 1.0)
InvGamma_TDist2 = MyInvGamma(t2nu, t2sigma)

@show InvGamma_Cauchy
@show InvGamma_TDist1
@show InvGamma_TDist2

@show mean(InvGamma_Cauchy)
@show mean(InvGamma_TDist1)
@show mean(InvGamma_TDist2)
println()
@show median(InvGamma_Cauchy)
@show median(InvGamma_TDist1)
@show median(InvGamma_TDist2)
println()
@show mode(InvGamma_Cauchy)
@show mode(InvGamma_TDist1)
@show mode(InvGamma_TDist2)
;

pdfvar_Cauchy(x) = pdf(InvGamma_Cauchy, x)
pdfvar_TDist1(x) = pdf(InvGamma_TDist1, x)
pdfvar_TDist2(x) = pdf(InvGamma_TDist2, x)

figure(figsize=(12,4))

subplot2grid((1,20), (0,0), rowspan=1, colspan=10)
xmax = 2.0
x = collect(linspace(eps(),xmax,1000))
plot([0],[0], label="Normal: \$\\delta_{$nsigma_str^2}\$", lw=1)
plot(x, pdfvar_Cauchy.(x), label="Cauchy: $csigma inv\$\\Gamma(1/2,$csigma_str^2/2)\$", lw=1.2, ls=":")
plot(x, pdfvar_TDist1.(x), label="TDist1: inv\$\\Gamma($t1nu_str/2, $t1nu_str/2)\$", lw=0.8, ls="--")
plot(x, pdfvar_TDist2.(x), label="TDist2: inv\$\\Gamma($t2nu_str/2, $t2nu_str\\cdot$t2sigma_str^2/2)\$", lw=0.8, ls="-.")
xlabel("variance of the normal distribution")
ylabel("probability density")
grid("on", ls=":")
legend(fontsize=9)
title("The predictive distributions of variances", fontsize=10)

subplot2grid((1,20), (0,10), rowspan=1, colspan=10)
x = collect(linspace(8.5,11.5,1000))
plot(x, pred_Normal.(x), label="Normal(\$\\mu=$nmu_str\$, \$\\sigma=$nsigma_str\$)", lw=0.8, "-")
plot(x, pred_Cauchy.(x), label="Cauchy(\$\\mu=$cmu_str\$, \$\\sigma=$csigma_str\$)", lw=1.2, ls=":")
plot(x, pred_TDist1.(x), label="$t1mu_str + TDist(\$\\nu=$t1nu_str\$)", lw=0.8, ls="--")
plot(x, pred_TDist2.(x), label="$t2mu_str + $t2sigma_str TDist(\$\\nu=$t2nu_str\$)", lw=0.8, ls="-.")
xlabel("y")
grid("on", linestyle=":")
legend(fontsize=9)
title("PDFs of the predictive distributions")

tight_layout();

InvGamma(nu, sigma) = InverseGamma(nu/2,nu*sigma^2/2)
cdfvar_Normal(x) = ifelse(x < nsigma^2, 0.0, 1.0)
cdfvar_Cauchy(x) = cdf(InvGamma_Cauchy, x)
cdfvar_TDist1(x) = cdf(InvGamma_TDist1, x)
cdfvar_TDist2(x) = cdf(InvGamma_TDist2, x)

figure(figsize=(12,4))

subplot2grid((1,20), (0,0), rowspan=1, colspan=10)
x = collect(linspace(eps(),10.0,1000))
plot(x, cdfvar_Normal.(x.^2), label="Normal: \$\\delta_{$nsigma_str^2}\$", lw=0.8, ls="-")
plot(x, cdfvar_Cauchy.(x.^2), label="Cauchy: inv\$\\Gamma(1/2,$csigma_str^2/2)\$", lw=1.2, ls=":")
plot(x, cdfvar_TDist1.(x.^2), label="TDist1: inv\$\\Gamma($t1nu_str/2, $t1nu_str/2)\$", lw=0.8, ls="--")
plot(x, cdfvar_TDist2.(x.^2), label="TDist2: inv\$\\Gamma($t2nu_str/2, $t2sigma_str^2\\cdot$t2nu_str/2)\$", lw=0.8, ls="-.")
xlabel("standard deviation of the normal distribution")
ylabel("cumulative probability")
grid("on", ls=":")
legend(fontsize=9)
title("The predictive distributions of standard deviations", fontsize=10)

subplot2grid((1,20), (0,10), rowspan=1, colspan=10)
x = collect(linspace(eps(),40,1000))
plot(x, cdfvar_Normal.(x.^2), label="Normal: \$\\delta_{$nsigma_str^2}\$", lw=0.8, ls="-")
plot(x, cdfvar_Cauchy.(x.^2), label="Cauchy: inv\$\\Gamma(1/2,$csigma_str^2/2)\$", lw=1.2, ls=":")
plot(x, cdfvar_TDist1.(x.^2), label="TDist1: inv\$\\Gamma($t1nu_str/2, $t1nu_str/2)\$", lw=0.8, ls="--")
plot(x, cdfvar_TDist2.(x.^2), label="TDist2: inv\$\\Gamma($t2nu_str/2, $t2sigma_str^2\\cdot$t2nu_str/2)\$", lw=0.8, ls="-.")
xlabel("standard deviation of the normal distribution")
ylabel("cumulative probability")
grid("on", ls=":")
legend(fontsize=9)
title("The predictive distributions of standard deviations", fontsize=10)

tight_layout();

InvGamma(nu, sigma) = InverseGamma(nu/2,nu*sigma^2/2)
cdfvar_Normal(x) = ifelse(x < nsigma^2, 0.0, 1.0)
cdfvar_Cauchy(x) = cdf(InvGamma_Cauchy, x)
cdfvar_TDist1(x) = cdf(InvGamma_TDist1, x)
cdfvar_TDist2(x) = cdf(InvGamma_TDist2, x)

figure(figsize=(12,4))

subplot2grid((1,20), (0,0), rowspan=1, colspan=10)
x = collect(linspace(eps(),2.0,1000))
plot(x, cdfvar_Normal.(x.^2), label="Normal: \$\\delta_{$nsigma_str^2}\$", lw=0.8, ls="-")
plot(x, cdfvar_Cauchy.(x.^2), label="Cauchy: inv\$\\Gamma(1/2,$csigma_str^2/2)\$", lw=1.2, ls=":")
plot(x, cdfvar_TDist1.(x.^2), label="TDist1: inv\$\\Gamma($t1nu_str/2, $t1nu_str/2)\$", lw=0.8, ls="--")
plot(x, cdfvar_TDist2.(x.^2), label="TDist2: inv\$\\Gamma($t2nu_str/2, $t2sigma_str^2\\cdot$t2nu_str/2)\$", lw=0.8, ls="-.")
xlabel("standard deviation of the normal distribution")
ylabel("cumulative probability")
grid("on", ls=":")
legend(fontsize=9)
title("The predictive distributions of standard deviations", fontsize=10)

subplot2grid((1,20), (0,10), rowspan=1, colspan=10)
x = collect(linspace(15,40,1000))
plot(x, cdfvar_Normal.(x.^2), label="Normal: \$\\delta_{$nsigma_str^2}\$", lw=0.8, ls="-")
plot(x, cdfvar_Cauchy.(x.^2), label="Cauchy: inv\$\\Gamma(1/2,$csigma_str^2/2)\$", lw=1.2, ls=":")
plot(x, cdfvar_TDist1.(x.^2), label="TDist1: inv\$\\Gamma($t1nu_str/2, $t1nu_str/2)\$", lw=0.8, ls="--")
plot(x, cdfvar_TDist2.(x.^2), label="TDist2: inv\$\\Gamma($t2nu_str/2, $t2sigma_str^2\\cdot$t2nu_str/2)\$", lw=0.8, ls="-.")
xlabel("standard deviation of the normal distribution")
ylabel("cumulative probability")
grid("on", ls=":")
legend(fontsize=9)
title("The predictive distributions of standard deviations", fontsize=10)

tight_layout();

@show InvGamma_Cauchy.invd.α, InvGamma_Cauchy.θ
@show InvGamma_TDist1.invd.α, InvGamma_TDist1.θ
@show InvGamma_TDist2.invd.α, InvGamma_TDist2.θ
;

function constpred_Normal(Y)
    f = fit_mle(Normal,Y)
    pred(x) = pdf(Normal(f.μ,f.σ),x)
    return pred
end

function maxloglikelihood_Normal(Y)
    pred = constpred_Normal(Y)
    n = length(Y)
    maxloglik = sum(log(pred(Y[i])) for i in 1:n)
    return maxloglik
end

aic_Normal(Y) = -2*maxloglikelihood_Normal(Y) + 2*2

function loocv_Normal(Y)
    loocv = 0.0
    n = length(Y)
    for i in 1:n
        pred = constpred_Normal(Y[[1:i-1;i+1:n]])
        loocv += log(pred(Y[i]))
    end
    return -2*loocv
end

@show -2*maxloglikelihood_Normal(Y)
@show aic_Normal(Y)
@show loocv_Normal(Y);

negloglik_Cauchy(Y,mu,sigma) = sigma > zero(sigma) ? -sum(logpdf(Cauchy(mu,sigma),Y[i]) for i in 1:length(Y)) : Inf
negloglik_Cauchy(Y,x) = negloglik_Cauchy(Y,x[1],x[2])
   
function constpred_Cauchy(Y; init=[0.0, 1.0])
    negloglik(x) = negloglik_Cauchy(Y,x)
    result = optimize(negloglik, init)
    a = result.minimizer
    pred(x) = pdf(Cauchy(a[1],a[2]),x)
    return pred
end

function maxloglikelihood_Cauchy(Y)
    pred = constpred_Cauchy(Y)
    n = length(Y)
    maxloglik = sum(log(pred(Y[i])) for i in 1:n)
    return maxloglik
end

aic_Cauchy(Y) = -2*maxloglikelihood_Cauchy(Y) + 2*2

function loocv_Cauchy(Y)
    loocv = 0.0
    n = length(Y)
    for i in 1:n
        pred = constpred_Cauchy(Y[[1:i-1;i+1:n]])
        loocv += log(pred(Y[i]))
    end
    return -2*loocv
end

@show -2*maxloglikelihood_Cauchy(Y)
@show aic_Cauchy(Y)
@show loocv_Cauchy(Y);

logpdf_TDist1(mu,nu,x) = nu > zero(nu) ? logpdf(TDist(nu), x-mu) : -Inf
negloglik_TDist1(Y,mu,nu) = -sum(logpdf_TDist1(mu,nu,Y[i]) for i in 1:length(Y))
negloglik_TDist1(Y,x) = negloglik_TDist1(Y,x[1],x[2])

function constpred_TDist1(Y; init=[0, 1.0])
    negloglik(x) = negloglik_TDist1(Y,x)
    result = optimize(negloglik, init)
    a = result.minimizer
    pred(x) = exp(logpdf_TDist1(a[1],a[2],x))
    return pred
end

function maxloglikelihood_TDist1(Y)
    pred = constpred_TDist1(Y)
    n = length(Y)
    maxloglik = sum(log(pred(Y[i])) for i in 1:n)
    return maxloglik
end

aic_TDist1(Y) = -2*maxloglikelihood_TDist1(Y) + 2*2

function loocv_TDist1(Y)
    loocv = 0.0
    n = length(Y)
    for i in 1:n
        pred = constpred_TDist1(Y[[1:i-1;i+1:n]])
        loocv += log(pred(Y[i]))
    end
    return -2*loocv
end

@show -2*maxloglikelihood_TDist1(Y)
@show aic_TDist1(Y)
@show loocv_TDist1(Y);

logpdf_TDist2(mu,nu,sigma,x) = (nu > zero(nu) && sigma > zero(nu)) ? logpdf(TDist(nu),(x-mu)/sigma)-log(sigma) : -Inf
negloglik_TDist2(Y,mu,nu,sigma) = -sum(logpdf_TDist2(mu,nu,sigma,Y[i]) for i in 1:length(Y))
negloglik_TDist2(Y,x) = negloglik_TDist2(Y,x[1],x[2],x[3])

function constpred_TDist2(Y; init=[0.0, 1.0, 1.0])
    negloglik(x) = negloglik_TDist2(Y,x)
    result = optimize(negloglik, init)
    a = result.minimizer
    pred(x) = exp(logpdf_TDist2(a[1],a[2],a[3],x))
    return pred
end

function maxloglikelihood_TDist2(Y)
    pred = constpred_TDist2(Y)
    n = length(Y)
    maxloglik = sum(log(pred(Y[i])) for i in 1:n)
    return maxloglik
end

aic_TDist2(Y) = -2*maxloglikelihood_TDist2(Y) + 2*3

function loocv_TDist2(Y)
    loocv = 0.0
    n = length(Y)
    for i in 1:n
        pred = constpred_TDist2(Y[[1:i-1;i+1:n]])
        loocv += log(pred(Y[i]))
    end
    return -2*loocv
end

@show -2*maxloglikelihood_TDist2(Y)
@show aic_TDist2(Y)
@show loocv_TDist2(Y);

println("-- Normal --")
@show aic_Normal(Y)
@show loocv_Normal(Y)
println("\n-- Cauchy --")
@show aic_Cauchy(Y)
@show loocv_Cauchy(Y)
println("\n-- TDist1 --")
@show aic_TDist1(Y)
@show loocv_TDist1(Y)
println("\n-- TDist2 --")
@show aic_TDist2(Y)
@show loocv_TDist2(Y)
;

p_Normal = [constpred_Normal(Y[[1:i-1;i+1:length(Y)]]) for i in 1:n]
loss_Normal = [-2*log(p_Normal[i](Y[i])) for i in 1:n]

p_Cauchy = [constpred_Cauchy(Y[[1:i-1;i+1:length(Y)]]) for i in 1:n]
loss_Cauchy = [-2*log(p_Cauchy[i](Y[i])) for i in 1:n]

p_TDist1 = [constpred_TDist1(Y[[1:i-1;i+1:length(Y)]]) for i in 1:n]
loss_TDist1 = [-2*log(p_TDist1[i](Y[i])) for i in 1:n]

p_TDist2 = [constpred_TDist2(Y[[1:i-1;i+1:length(Y)]]) for i in 1:n]
loss_TDist2 = [-2*log(p_TDist2[i](Y[i])) for i in 1:n]

function cyc_color(i)
    colorlist = [
        "#1f77b4",
        "#ff7f0e",
        "#2ca02c",
        "#d62728",
        "#9467bd",
        "#8c564b",
        "#e377c2",
        "#7f7f7f",
        "#bcbd22",
        "#17becf",
    ]
    n = length(colorlist)
    return colorlist[mod1(i,n)]
end

function cyc_style(i)
    stylelist = [
        "dashed",
        "dashdot",
        "dotted",
        (0, (4.5, 1.5, 1, 1.5, 1, 1.5, 1, 1.5)),
        (0, (3, 1.5, 1, 1.5, 3, 1.5, 1, 1.5)),
        (0, (2.5, 1.5, 2.5, 1.5, 2.5, 1.5, 1, 1.5)),
        (0, (2.5, 1, 2.5, 1, 2.5, 1, 2.5, 1)),
        (0, (5, 1.5, 2, 1.5, 2, 1.5)),
        (0, (5, 2, 2)),
        (0, (6, 3)),
    ]
    n = length(stylelist)
    return stylelist[mod1(i,n)]
end

@assert n == 7

figure(figsize=(12,5))
x = collect(linspace(-30,40,1000))

shape = (17,32)
loc(i) = (1+8*div(i,4),8*mod(i,4))
rs, cs = 8, 8

ax = subplot2grid(shape, loc(0), rowspan=rs, colspan=cs)
plot(x, -2*log.(pred_Normal.(x)), label="\$-2\\log p_\\ast\$", color="black", ls="-")
scatter(Y, -2*log.(pred_Normal.(Y)), color="black")
grid("on", linestyle=":")
legend()

for i in 1:n
    ax = subplot2grid(shape, loc(i), rowspan=rs, colspan=cs)
    plot(x, -2*log.(p_Normal[i].(x)), label="$i", color=cyc_color(i), ls=cyc_style(i))
    scatter([Y[i]], [-2*log(p_Normal[i](Y[i]))], color=cyc_color(i))
    grid("on", linestyle=":")
    legend()
end

suptitle("-2 log(PDFs of the leave-one-out predictive distributions for the Normal distribution model)")

tight_layout()

@assert n == 7

figure(figsize=(12,5))
x = collect(linspace(-30,40,1000))

shape = (17,32)
loc(i) = (1+8*div(i,4),8*mod(i,4))
rs, cs = 8, 8

ax = subplot2grid(shape, loc(0), rowspan=rs, colspan=cs)
plot(x, -2*log.(pred_Cauchy.(x)), label="\$-2\\log p_\\ast\$", color="black", ls="-")
scatter(Y, -2*log.(pred_Cauchy.(Y)), color="black")
grid("on", linestyle=":")
legend()

for i in 1:n
    ax = subplot2grid(shape, loc(i), rowspan=rs, colspan=cs)
    plot(x, -2*log.(p_Cauchy[i].(x)), label="$i", color=cyc_color(i), ls=cyc_style(i))
    scatter([Y[i]], [-2*log(p_Cauchy[i](Y[i]))], color=cyc_color(i))
    grid("on", linestyle=":")
    legend()
end

suptitle("-2 log(PDFs of the leave-one-out predictive distributions for the Cauchy distribution model)")

tight_layout()

@assert n == 7

figure(figsize=(12,5))
x = collect(linspace(-30,40,1000))

shape = (17,32)
loc(i) = (1+8*div(i,4),8*mod(i,4))
rs, cs = 8, 8

ax = subplot2grid(shape, loc(0), rowspan=rs, colspan=cs)
plot(x, -2*log.(pred_TDist1.(x)), label="\$-2\\log p_\\ast\$", color="black", ls="-")
scatter(Y, -2*log.(pred_TDist1.(Y)), color="black")
grid("on", linestyle=":")
legend()

for i in 1:n
    ax = subplot2grid(shape, loc(i), rowspan=rs, colspan=cs)
    plot(x, -2*log.(p_TDist1[i].(x)), label="$i", color=cyc_color(i), ls=cyc_style(i))
    scatter([Y[i]], [-2*log(p_TDist1[i](Y[i]))], color=cyc_color(i))
    grid("on", linestyle=":")
    legend()
end

suptitle("-2 log(PDFs of the leave-one-out predictive distributions for the t-distribution model 1)")

tight_layout()

@assert n == 7

figure(figsize=(12,5))
x = collect(linspace(-30,40,1000))

shape = (17,32)
loc(i) = (1+8*div(i,4),8*mod(i,4))
rs, cs = 8, 8

ax = subplot2grid(shape, loc(0), rowspan=rs, colspan=cs)
plot(x, -2*log.(pred_TDist2.(x)), label="\$-2\\log p_\\ast\$", color="black", ls="-")
scatter(Y, -2*log.(pred_TDist2.(Y)), color="black")
grid("on", linestyle=":")
legend()

for i in 1:n
    ax = subplot2grid(shape, loc(i), rowspan=rs, colspan=cs)
    plot(x, -2*log.(p_TDist2[i].(x)), label="$i", color=cyc_color(i), ls=cyc_style(i))
    scatter([Y[i]], [-2*log(p_TDist2[i](Y[i]))], color=cyc_color(i))
    grid("on", linestyle=":")
    legend()
end

suptitle("-2 log(PDFs of the leave-one-out predictive distributions for the t-distribution model 2)")

tight_layout()

function constparam_Normal(Y)
    d = fit_mle(Normal,Y)
    return d.μ, d.σ
end

function constparam_Cauchy(Y; init=[0.0, 1.0])
    negloglik(x) = negloglik_Cauchy(Y,x)
    result = optimize(negloglik, init)
    a = result.minimizer
    return a[1],a[2]
end

function constparam_TDist1(Y; init=[0.0, 1.0])
    negloglik(x) = negloglik_TDist1(Y,x)
    result = optimize(negloglik, init)
    a = result.minimizer
    return a[1],a[2]
end

function constparam_TDist2(Y; init=[0.0, 1.0, 1.0])
    negloglik(x) = negloglik_TDist2(Y,x)
    result = optimize(negloglik, init)
    a = result.minimizer
    return a[1],a[2],a[3]
end

Y

YY = [Y; [9.801, 9.701, 9.958]]

function comp_Normal(Y,YY;xmax=2)
    
@show aic_Normal(Y)/(2*length(Y))
@show loocv_Normal(Y)/(2*length(Y))
println()
@show aic_Normal(YY)/(2*length(YY))
@show loocv_Normal(YY)/(2*length(YY))
println()
@show a = constparam_Normal(Y)
@show nmu = a[1]
@show nsigma = a[2]
cdfvar_Normal(x) = x < nsigma^2 ? 0.0 : 1.0
println()
@show a = constparam_Normal(YY)
@show nnmu = a[1]
@show nnsigma = a[2]
cdfvar_NNormal(x) = x < nnsigma^2 ? 0.0 : 1.0
sleep(0.2)

figure(figsize=(6,3.6))
x = collect(linspace(0,xmax,1000))
plot(x, cdfvar_Normal.(x.^2),  label="Y")
plot(x, cdfvar_NNormal.(x.^2), label="YY", ls="--")
xlabel("std")
grid("on", ls=":")
title("PDFs")
legend()

end

comp_Normal(Y,YY,xmax=25)

function comp_Cauchy(Y,YY)

@show aic_Cauchy(Y)/(2*length(Y))
@show loocv_Cauchy(Y)/(2*length(Y))
println()
@show aic_Cauchy(YY)/(2*length(YY))
@show loocv_Cauchy(YY)/(2*length(YY))
println()
@show a = constparam_Cauchy(Y)
@show cmu = a[1]
@show csigma = a[2]
cdfvar_Cauchy(x) = cdf(InverseGamma(1/2,csigma^2/2), x)
pdfvar_Cauchy(x) = pdf(InverseGamma(1/2,csigma^2/2), x)
println()
@show a = constparam_Cauchy(YY)
@show ccmu = a[1]
@show ccsigma = a[2]
cdfvar_CCauchy(x) = cdf(InverseGamma(1/2,ccsigma^2/2), x)
pdfvar_CCauchy(x) = pdf(InverseGamma(1/2,ccsigma^2/2), x)
sleep(0.2)

figure(figsize=(9,3.6))
    
subplot(1,2,1)
x = collect(linspace(0,2,1000))
plot(x, pdfvar_Cauchy.(x), label="Y")
plot(x, pdfvar_CCauchy.(x),  label="YY", ls="--")
xlabel("variance")
grid("on", ls=":")
title("PDFs")
legend()

subplot(1,2,2)
x = collect(linspace(0,2,1000))
plot(x, cdfvar_Cauchy.(x.^2), label="Y")
plot(x, cdfvar_CCauchy.(x.^2),  label="YY", ls="--")
xlabel("std")
grid("on", ls=":")
title("CDFs")
legend()

tight_layout()    

end

comp_Cauchy(Y,YY)

function comp_TDist1(Y,YY)

@show aic_TDist1(Y)/(2*length(Y))
@show loocv_TDist1(Y)/(2*length(Y))
println()
@show aic_TDist1(YY)/(2*length(YY))
@show loocv_TDist1(YY)/(2*length(YY))
println()
@show a = constparam_TDist1(Y)
@show td1mu = a[1]
@show td1nu = a[2]
cdfvar_TDist1(x) = cdf(InverseGamma(td1nu/2,td1nu/2), x)
pdfvar_TDist1(x) = pdf(InverseGamma(td1nu/2,td1nu/2), x)
println()
@show a = constparam_TDist1(YY)
@show ttd1mu = a[1]
@show ttd1nu = a[2]
cdfvar_TTDist1(x) = cdf(InverseGamma(ttd1nu/2,ttd1nu/2), x)
pdfvar_TTDist1(x) = pdf(InverseGamma(ttd1nu/2,ttd1nu/2), x)
sleep(0.2)

figure(figsize=(9,3.6))
    
subplot(1,2,1)
x = collect(linspace(0,2,1000))
plot(x, pdfvar_TDist1.(x),  label="Y")
plot(x, pdfvar_TTDist1.(x), label="YY", ls="--")
xlabel("variance")
grid("on", ls=":")
title("PDFs")
legend()

subplot(1,2,2)
x = collect(linspace(0,2,1000))
plot(x, cdfvar_TDist1.(x.^2),  label="Y")
plot(x, cdfvar_TTDist1.(x.^2), label="YY", ls="--")
xlabel("std")
grid("on", ls=":")
title("CDFs")
legend()

tight_layout()
    
end

comp_TDist1(Y,YY)

function comp_TDist2(Y,YY)

@show aic_TDist2(Y)/(2*length(Y))
@show loocv_TDist2(Y)/(2*length(Y))
println()
@show aic_TDist2(YY)/(2*length(YY))
@show loocv_TDist2(YY)/(2*length(YY))
println()
@show a = constparam_TDist2(Y)
@show td2mu = a[1]
@show td2nu = a[2]
@show td2sigma = a[3]
cdfvar_TDist2(x) = cdf(InverseGamma(td2nu/2,td2sigma^2*td2nu/2), x)
pdfvar_TDist2(x) = pdf(InverseGamma(td2nu/2,td2sigma^2*td2nu/2), x)
println()
@show a = constparam_TDist2(YY)
@show ttd2mu = a[1]
@show ttd2nu = a[2]
@show ttd2sigma = a[3]
cdfvar_TTDist2(x) = cdf(InverseGamma(ttd2nu/2,ttd2sigma^2*ttd2nu/2), x)
pdfvar_TTDist2(x) = pdf(InverseGamma(ttd2nu/2,ttd2sigma^2*ttd2nu/2), x)
sleep(0.2)

figure(figsize=(9,3.6))
    
subplot(1,2,1)
x = collect(linspace(0,2,1000))
plot(x, pdfvar_TDist2.(x),  label="Y",)
plot(x, pdfvar_TTDist2.(x), label="YY", ls="--")
xlabel("variance")
grid("on", ls=":")
title("PDFs")
legend()

subplot(1,2,2)
x = collect(linspace(0,2,1000))
plot(x, cdfvar_TDist2.(x.^2),  label="Y")
plot(x, cdfvar_TTDist2.(x.^2), label="YY", ls="--")
xlabel("std")
grid("on", ls=":")
title("CDFs")
legend()

tight_layout()
    
end

comp_TDist2(Y,YY)

YYY  = 0.5*randn(128)
YYYY = 0.5*randn(512)

@show mean(YYY), std(YYY)
@show mean(YYYY), std(YYYY);

comp_Normal(YYY,YYYY,xmax=1)

comp_Cauchy(YYY,YYYY)

comp_TDist1(YYY,YYYY)

comp_TDist2(YYY,YYYY)

println("-- Normal --")
@show aic_Normal(YYYY)
@show loocv_Normal(YYYY)
println("\n-- Cauchy --")
@show aic_Cauchy(YYYY)
@show loocv_Cauchy(YYYY)
println("\n-- TDist1 --")
@show aic_TDist1(YYYY)
@show loocv_TDist1(YYYY)
println("\n-- TDist2 --")
@show aic_TDist2(YYYY)
@show loocv_TDist2(YYYY)
;

Ymu = 0.0
Ynu = 1.0
Ysigma = 0.5
YT = Ymu .+ Ysigma .* rand(TDist(Ynu),256)
YN = Ymu .+ Ysigma .* randn(256)
@show Ymu, Ynu, Ysigma
@show mean(YT), std(YT)
@show mean(YN), std(YN);

@time comp_Normal(YT, YN, xmax=1)

@time comp_TDist2(YT,YN)