Rel 0.0.4 - Added waic() and pk utilities

Author: goedman

Project.toml

@@ -1,20 +1,20 @@
 name = "PSIS"
 uuid = "fb740163-aa3c-59c1-9c12-c3f890714cde"
-authors = ["Rob J Goedman <goedman@icloud.com"]
+authors = ["@alvaro1101, Rob J Goedman <goedman@icloud.com"]
 version = "0.0.3"
 
 [deps]
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StanSample = "c1514b29-d3a0-5178-b312-660c88baa699"
-StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 StatisticalRethinking = "2d09df54-9d0f-5258-8220-54c2a3d4fbee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 JSON = "0.21"
-StatsPlots = "0.14"
 StanSample = "3.0"
 StatisticalRethinking = "3.2"
+StatsPlots = "0.14"
 julia = "1"

src/PSIS.jl

@@ -39,6 +39,8 @@ include("psislw.jl")
 include("gpdfitnew.jl")
 include("gpinv.jl")
 include("logsumexp.jl")
+include("waic.jl")
+include("pk_utilities.jl")
 
 export
     psis_path,

src/pk_utilities.jl

@@ -0,0 +1,18 @@
+function pk_qualify(pk::Vector{Float64})
+    pk_good = sum(pk .<= 0.5)
+    pk_ok = length(pk[pk .<= 0.7]) - pk_good
+    pk_bad = length(pk[pk .<= 1]) - pk_good - pk_ok
+    (good=pk_good, ok=pk_ok, bad=pk_bad, very_bad=sum(pk .> 1))
+end
+
+function pk_plot(pk::Vector{Float64}, title="PSIS diagnostic plot.")
+      scatter(pk1, xlab="Datapoint", ylab="Pareto shape k",
+    marker=2.5, lab="Pk points", leg=:topleft)
+  hline!([0.5], lab="pk = 0.5");hline!([0.7], lab="pk = 0.7")
+  hline!([1], lab="pk = 1.0")
+  title!(title)
+end
+
+export
+    pk_qualify,
+    pk_plot

src/waic.jl

@@ -0,0 +1,43 @@
+
+var2(x) = mean(x.^2) .- mean(x)^2
+
+function log_sum_exp(x) 
+    xmax = maximum(x)
+    xsum = sum(exp.(x .- xmax))
+    xmax + log(xsum)
+end
+
+function waic( ll::AbstractArray; pointwise=FALSE , log_lik="log_lik" , kwargs... )
+    
+    n_samples, n_obs = size(ll)
+    lpd <- zeros(n_obs)
+    pD <- zeros(n_obs)
+
+    for i in 1:n_obs 
+        lpd[i] = log_sum_exp(ll[:,i]) .- log(n_samples)
+        pD[i] = var2(ll[:,i])
+    end
+
+    waic_vec = (-2) .* ( lpd - pD )
+    if pointwise==FALSE
+        waic = sum(waic_vec)
+        lpd = sum(lpd)
+        pD = sum(pD)
+    else 
+        waic = waic_vec
+    end
+
+    try 
+        se = sqrt( n_obs*var2(waic_vec) )
+    catch e
+        prinrln(e)
+        se = nothing
+    end
+
+    (WAIC=waic, lppd=lpd, penalty=pD, std_err=se)
+end
+
+export
+    var2,
+    log_sum_exp,
+    waic