Rel 0.0.7 - Rely in StatsFuns' logsumexp

Author: goedman

Project.toml

@@ -9,12 +9,14 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 StanSample = "c1514b29-d3a0-5178-b312-660c88baa699"
 StatisticalRethinking = "2d09df54-9d0f-5258-8220-54c2a3d4fbee"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 JSON = "0.21"
 StanSample = "3.0"
 StatisticalRethinking = "3.2"
+StatsFuns = "0.9"
 StatsPlots = "0.14"
 julia = "1"

README.md

@@ -50,9 +50,6 @@ The cars WAIC example requires RDatasets.jl to be installed and functioning.
 `gpinv()` -
     Inverse Generalised Pareto distribution function.
 
-`logsumexp()` -
-    Sum of a vector where numbers are represented by their logarithms.
-
 `waic()` -
     Compute WAIC for a loglikelihood matrix.
 
@@ -65,9 +62,6 @@ The cars WAIC example requires RDatasets.jl to be installed and functioning.
 `var2()` -
     Uncorrected variance.
 
-`log_sum_exp()` -
-    Compute logarithmic sum of a vector.
-
 ### Corresponding R code
 
 The corresponding R code can be found in [R package called

examples/cars_waic/cars.jl

@@ -37,14 +37,13 @@ rc = stan_sample(cars_stan_model; data)
 if success(rc)
     cars_df = read_samples(cars_stan_model; output_format=:dataframe)
     precis(cars_df[:, [:a, :b, :sigma]])
-
     nt_cars = read_samples(cars_stan_model);
 end
 
 log_lik = nt_cars.log_lik'
-ns, n = size(log_lik)
-lppd = [log_sum_exp(log_lik[:, i] .- log(ns)) for i in 1:n]
-pwaic = [var(log_lik[:, i]) for i in 1:n]
+n_sam, n_obs = size(log_lik)
+lppd = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs)
+pwaic = [var(log_lik[:, i]) for i in 1:n_obs]
 -2(sum(lppd) - sum(pwaic)) |> display
 
 waic(log_lik) |> display

notebooks/cars.jl

@@ -75,12 +75,12 @@ end
 if success(rc)
 	nt_cars = read_samples(cars_stan_model);
 	log_lik = nt_cars.log_lik'
-	ns, n = size(log_lik)
 end
 
 # ╔═╡ 20ed768a-6008-11eb-13f4-458ca1a29592
 begin
-	lppd = [StatsFuns.logsumexp(log_lik[:, i] .- log(ns)) for i in 1:n]
+    ns, n = size(log_lik)
+	lppd = reshape(logsumexp(log_lik .- log(ns); dims=1), n)
 	pwaic = [var(log_lik[:, i]) for i in 1:n]
 	-2(sum(lppd) - sum(pwaic))
 end

src/PSIS.jl

@@ -32,13 +32,15 @@ arXiv preprint arXiv:1507.02646.
 
 module PSIS
 
+using StatsFuns
+
 psis_path = @__DIR__
 
 include("psisloo.jl")
 include("psislw.jl")
 include("gpdfitnew.jl")
 include("gpinv.jl")
-include("logsumexp.jl")
+#include("logsumexp.jl")
 include("waic.jl")
 include("pk_utilities.jl")
 

src/psisloo.jl

@@ -38,7 +38,7 @@ function psisloo(log_lik::AbstractArray, wcpp::Int64=20, wtrunc::Float64=3/4)
     lwp, ks = psislw(-lw, wcpp, wtrunc)
 
     lwp += lw
-    loos = logsumexp(lwp, 1)
+    loos = reshape(logsumexp(lwp, dims=1), size(lw, 2))
     loo = sum(loos)
 
     return loo, loos, ks

src/waic.jl

@@ -4,11 +4,11 @@ var2(x) = mean(x.^2) .- mean(x)^2
 function waic( ll::AbstractArray; pointwise=false , log_lik="log_lik" , kwargs... )
 
     n_samples, n_obs = size(ll)
-    lpd = zeros(n_obs)
+    #lpd = zeros(n_obs)
     pD = zeros(n_obs)
 
+    lpd = reshape(logsumexp(ll .- log(n_samples); dims=1), n_obs)
     for i in 1:n_obs 
-        lpd[i] = logsumexp(ll[:,i]) .- log(n_samples)
         pD[i] = var2(ll[:,i])
     end