Bug fix

Author: Closed-Limelike-Curves

.lh/src/ESS.jl.json

@@ -3,7 +3,7 @@
     "activeCommit": 0,
     "commits": [
         {
-            "activePatchIndex": 11,
+            "activePatchIndex": 12,
             "patches": [
                 {
                     "date": 1626309971342,
@@ -52,6 +52,10 @@
                 {
                     "date": 1626454525488,
                     "content": "Index: \n===================================================================\n--- \n+++ \n@@ -2,9 +2,9 @@\n using MCMCDiagnosticTools\n using LoopVectorization\n using Tullio\n \n-export relative_eff, psis_ess\n+export relative_eff, psis_ess, psis_n_eff\n \n \"\"\"\n     relative_eff(sample::AbstractArray{AbstractFloat, 3}; method=FFTESSMethod())\n \n@@ -56,4 +56,8 @@\n     @warn \"PSIS ESS not adjusted based on MCMC ESS. MCSE and ESS estimates \" *\n           \"will be overoptimistic if samples are autocorrelated.\"\n     return psis_ess(weights, ones(size(weights)))\n end\n+\n+function psis_n_eff(args...; kwargs...)  # Alias for compatibility with R version\n+    return psis_ess(args...; kwargs...)\n+end\n\\ No newline at end of file\n"
+                },
+                {
+                    "date": 1626484291838,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -11,10 +11,10 @@\n Calculate the relative efficiency of an MCMC chain, i.e. the effective sample size divided\n by the nominal sample size.\n \"\"\"\n function relative_eff(\n-    sample::AbstractArray{T,3}; method=FFTESSMethod()\n-) where {T<:AbstractFloat}\n+    sample::AbstractArray{T, 3}; method=FFTESSMethod()\n+) where {T <: AbstractFloat}\n     dims = size(sample)\n     post_sample_size = dims[2] * dims[3]\n     # Only need ESS, not rhat\n     ess_sample = inv.(permutedims(sample, [2, 1, 3]))\n@@ -33,27 +33,27 @@\n Vehtari et al. 2019.\n \n # Arguments\n \n-- `weights`: A set of importance sampling weights derived from PSIS.\n-- `r_eff`: The relative efficiency of the MCMC chains from which PSIS samples were derived.\n-See `?relative_eff` to calculate `r_eff`.\n+  - `weights`: A set of importance sampling weights derived from PSIS.\n+  - `r_eff`: The relative efficiency of the MCMC chains from which PSIS samples were derived.\n+    See `?relative_eff` to calculate `r_eff`.\n \"\"\"\n function psis_ess(\n     weights::AbstractVector{T}, r_eff::AbstractVector{T}\n-) where {T<:AbstractFloat}\n+) where {T <: AbstractFloat}\n     @tullio sum_of_squares := weights[x]^2\n     return r_eff ./ sum_of_squares\n end\n \n function psis_ess(\n     weights::AbstractMatrix{T}, r_eff::AbstractVector{T}\n-) where {T<:AbstractFloat}\n+) where {T <: AbstractFloat}\n     @tullio sum_of_squares[x] := weights[x, y]^2\n     return @tturbo r_eff ./ sum_of_squares\n end\n \n-function psis_ess(weights::AbstractMatrix{T}) where {T<:AbstractFloat}\n+function psis_ess(weights::AbstractMatrix{T}) where {T <: AbstractFloat}\n     @warn \"PSIS ESS not adjusted based on MCMC ESS. MCSE and ESS estimates \" *\n           \"will be overoptimistic if samples are autocorrelated.\"\n     return psis_ess(weights, ones(size(weights)))\n end\n"
                 }
             ],
             "date": 1626309971342,

.lh/src/GPD.jl.json

@@ -3,7 +3,7 @@
     "activeCommit": 0,
     "commits": [
         {
-            "activePatchIndex": 7,
+            "activePatchIndex": 14,
             "patches": [
                 {
                     "date": 1626310029646,
@@ -36,6 +36,34 @@
                 {
                     "date": 1626362509519,
                     "content": "Index: \n===================================================================\n--- \n+++ \n@@ -44,9 +44,8 @@\n \n \n     prior = 3\n     grid_size = min_grid_pts + isqrt(len) # isqrt = floor sqrt\n-    grid_size = grid_size - (grid_size % 32) # multiples of 32 speed up calcs\n     n_0 = 10  # determines how strongly to nudge ξ towards .5\n     quartile::T = sample[(len+2)÷4]\n \n \n"
+                },
+                {
+                    "date": 1626484293993,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -33,9 +33,9 @@\n     sample::AbstractVector{T};\n     wip::Bool=true,\n     min_grid_pts::Integer=30,\n     sort_sample::Bool=false,\n-) where {T<:AbstractFloat}\n+) where {T <: AbstractFloat}\n \n     len = length(sample)\n     # sample must be sorted, but we can skip if sample is already sorted\n     if sort_sample\n@@ -52,18 +52,18 @@\n     # build pointwise estimates of ξ and θ at each grid point\n     θ_hats = similar(sample, grid_size)\n     ξ_hats = similar(sample, grid_size)\n     @turbo @. θ_hats =\n-        inv(sample[len]) + (1 - sqrt((grid_size+1) / $(1:grid_size))) / prior / quartile\n-    @tullio threads=false ξ_hats[x] := log1p(-θ_hats[x] * sample[y]) |> _ / len\n+        inv(sample[len]) + (1 - sqrt((grid_size + 1) / $(1:grid_size))) / prior / quartile\n+    @tullio threads = false ξ_hats[x] := log1p(-θ_hats[x] * sample[y]) |> _ / len\n     log_like = similar(ξ_hats)\n     # Calculate profile log-likelihood at each estimate:\n     @turbo @. log_like = len * (log(-θ_hats / ξ_hats) - ξ_hats - 1)\n     # Calculate weights from log-likelihood:\n     weights = ξ_hats  # Reuse preallocated array (which is no longer in use)\n-    @tullio threads=false weights[y] = exp(log_like[x] - log_like[y]) |> inv\n+    @tullio threads = false weights[y] = exp(log_like[x] - log_like[y]) |> inv\n     # Take weighted mean:\n-    @tullio threads=false θ_hat := weights[x] * θ_hats[x]\n+    @tullio threads = false θ_hat := weights[x] * θ_hats[x]\n \n     ξ::T = calc_ξ(sample, θ_hat)\n     σ::T = -ξ / θ_hat\n \n@@ -90,9 +90,9 @@\n # Returns\n \n A quantile of the Generalized Pareto Distribution.\n \"\"\"\n-function gpd_quantile(p, k::T, sigma::T) where {T<:AbstractFloat}\n+function gpd_quantile(p, k::T, sigma::T) where {T <: AbstractFloat}\n     return sigma * expm1(-k * log1p(-p)) / k\n end\n \n \n@@ -100,9 +100,9 @@\n     calc_ξ(sample, θ_hat)\n \n Calculate ξ, the parameter for the GPD.\n \"\"\"\n-function calc_ξ(sample::AbstractVector{T}, θ_hat::T) where {T<:AbstractFloat}\n+function calc_ξ(sample::AbstractVector{T}, θ_hat::T) where {T <: AbstractFloat}\n     ξ = zero(T)\n     @turbo for i in eachindex(sample)\n         ξ += log1p(-θ_hat * sample[i]) / length(sample)\n     end\n"
+                },
+                {
+                    "date": 1626484439615,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -52,18 +52,18 @@\n     # build pointwise estimates of ξ and θ at each grid point\n     θ_hats = similar(sample, grid_size)\n     ξ_hats = similar(sample, grid_size)\n     @turbo @. θ_hats =\n-        inv(sample[len]) + (1 - sqrt((grid_size + 1) / $(1:grid_size))) / prior / quartile\n-    @tullio threads = false ξ_hats[x] := log1p(-θ_hats[x] * sample[y]) |> _ / len\n+        inv(sample[len]) + (1 - sqrt((grid_size+1) / $(1:grid_size))) / prior / quartile\n+    @tullio threads=false ξ_hats[x] := log1p(-θ_hats[x] * sample[y]) |> _ / len\n     log_like = similar(ξ_hats)\n     # Calculate profile log-likelihood at each estimate:\n     @turbo @. log_like = len * (log(-θ_hats / ξ_hats) - ξ_hats - 1)\n     # Calculate weights from log-likelihood:\n     weights = ξ_hats  # Reuse preallocated array (which is no longer in use)\n-    @tullio threads = false weights[y] = exp(log_like[x] - log_like[y]) |> inv\n+    @tullio threads=false weights[y] = exp(log_like[x] - log_like[y]) |> inv\n     # Take weighted mean:\n-    @tullio threads = false θ_hat := weights[x] * θ_hats[x]\n+    @tullio threads=false θ_hat := weights[x] * θ_hats[x]\n \n     ξ::T = calc_ξ(sample, θ_hat)\n     σ::T = -ξ / θ_hat\n \n"
+                },
+                {
+                    "date": 1626484457349,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -101,11 +101,11 @@\n \n Calculate ξ, the parameter for the GPD.\n \"\"\"\n function calc_ξ(sample::AbstractVector{T}, θ_hat::T) where {T <: AbstractFloat}\n-    ξ = zero(T)\n+    ξ::T = zero(T)\n     @turbo for i in eachindex(sample)\n         ξ += log1p(-θ_hat * sample[i]) / length(sample)\n     end\n-    return ξ::T\n+    return ξ\n end\n \n"
+                },
+                {
+                    "date": 1626484527691,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -62,8 +62,9 @@\n     weights = ξ_hats  # Reuse preallocated array (which is no longer in use)\n     @tullio threads=false weights[y] = exp(log_like[x] - log_like[y]) |> inv\n     # Take weighted mean:\n     @tullio threads=false θ_hat := weights[x] * θ_hats[x]\n+    @tullio threads=false ξ := log1p(-θ_hat * sample[i]) |> _ / len\n \n     ξ::T = calc_ξ(sample, θ_hat)\n     σ::T = -ξ / θ_hat\n \n@@ -94,18 +95,4 @@\n function gpd_quantile(p, k::T, sigma::T) where {T <: AbstractFloat}\n     return sigma * expm1(-k * log1p(-p)) / k\n end\n \n-\n-\"\"\"\n-    calc_ξ(sample, θ_hat)\n-\n-Calculate ξ, the parameter for the GPD.\n-\"\"\"\n-function calc_ξ(sample::AbstractVector{T}, θ_hat::T) where {T <: AbstractFloat}\n-    ξ::T = zero(T)\n-    @turbo for i in eachindex(sample)\n-        ξ += log1p(-θ_hat * sample[i]) / length(sample)\n-    end\n-    return ξ\n-end\n-\n"
+                },
+                {
+                    "date": 1626486963522,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -9,9 +9,9 @@\n         sample::AbstractVector{T<:AbstractFloat}; \n         wip::Bool=true, \n         min_grid_pts::Integer=30, \n         sort_sample::Bool=false\n-        ) -> (ξ::T, σ::T)\n+    ) -> (ξ::T, σ::T)\n \n Return a named list of estimates for the parameters ξ (shape) and σ (scale) of the\n generalized Pareto distribution (GPD), assuming the location parameter is 0.\n \n"
+                },
+                {
+                    "date": 1626487070812,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -64,9 +64,8 @@\n     # Take weighted mean:\n     @tullio threads=false θ_hat := weights[x] * θ_hats[x]\n     @tullio threads=false ξ := log1p(-θ_hat * sample[i]) |> _ / len\n \n-    ξ::T = calc_ξ(sample, θ_hat)\n     σ::T = -ξ / θ_hat\n \n     # Drag towards .5 to reduce variance for small len\n     if wip\n"
+                },
+                {
+                    "date": 1626487134955,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -62,10 +62,10 @@\n     weights = ξ_hats  # Reuse preallocated array (which is no longer in use)\n     @tullio threads=false weights[y] = exp(log_like[x] - log_like[y]) |> inv\n     # Take weighted mean:\n     @tullio threads=false θ_hat := weights[x] * θ_hats[x]\n-    @tullio threads=false ξ := log1p(-θ_hat * sample[i]) |> _ / len\n-\n+    @tullio threads=false ξ := log1p(-θ_hat * sample[i])\n+    ξ /= len\n     σ::T = -ξ / θ_hat\n \n     # Drag towards .5 to reduce variance for small len\n     if wip\n"
                 }
             ],
             "date": 1626310029646,

.lh/src/ImportanceSampling.jl.json

@@ -3,7 +3,7 @@
     "activeCommit": 0,
     "commits": [
         {
-            "activePatchIndex": 36,
+            "activePatchIndex": 43,
             "patches": [
                 {
                     "date": 1626309654641,
@@ -152,6 +152,34 @@
                 {
                     "date": 1626454092733,
                     "content": "Index: \n===================================================================\n--- \n+++ \n@@ -71,9 +71,9 @@\n     end\n     \n     @tullio norm_const[i] := weights[i, j]\n     weights .= weights ./ norm_const\n-    ess = psis_n_eff(weights, r_eff)\n+    ess = psis_ess(weights, r_eff)\n \n     weights = reshape(weights, dims)\n \n     if log_weights\n"
+                },
+                {
+                    "date": 1626458228667,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -60,9 +60,9 @@\n     # Shift ratios by maximum to prevent overflow\n     @tturbo @. weights = exp(log_ratios - $maximum(log_ratios; dims=2))\n \n     r_eff = _generate_r_eff(weights, dims, r_eff, source)\n-    check_input_validity_psis(reshape(log_ratios, dims), r_eff)\n+    _check_input_validity_psis(reshape(log_ratios, dims), r_eff)\n \n     tail_length = similar(log_ratios, Int, data_size)\n     ξ = similar(log_ratios, data_size)\n     @inbounds Threads.@threads for i in eachindex(tail_length)\n@@ -134,9 +134,9 @@\n \n     # Define and check tail\n     tail_start = len - tail_length + 1  # index of smallest tail value\n     @views tail = sorted_ratios[tail_start:len]\n-    check_tail(tail)\n+    _check_tail(tail)\n \n     # Get value just before the tail starts:\n     cutoff = sorted_ratios[tail_start - 1]\n     ξ = psis_smooth_tail!(tail, cutoff)\n@@ -246,9 +246,9 @@\n \n \"\"\"\n Make sure all inputs to `psis` are valid.\n \"\"\"\n-function check_input_validity_psis(\n+function _check_input_validity_psis(\n     log_ratios::AbstractArray{T,3}, r_eff::AbstractVector{T}\n ) where {T<:AbstractFloat}\n     if any(isnan, log_ratios)\n         throw(DomainError(\"Invalid input for `log_ratios` (contains NaN values).\"))\n@@ -271,9 +271,9 @@\n \n \"\"\"\n Check the tail to make sure a GPD fit is possible.\n \"\"\"\n-function check_tail(tail::AbstractVector{T}) where {T<:AbstractFloat}\n+function _check_tail(tail::AbstractVector{T}) where {T<:AbstractFloat}\n     if maximum(tail) ≈ minimum(tail)\n         throw(\n             ArgumentError(\n                 \"Unable to fit generalized Pareto distribution: all tail values are the\" *\n"
+                },
+                {
+                    "date": 1626484622160,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -242,9 +242,9 @@\n         return r_eff\n     end\n end\n \n-\n+    \n \"\"\"\n Make sure all inputs to `psis` are valid.\n \"\"\"\n function _check_input_validity_psis(\n"
+                },
+                {
+                    "date": 1626484665188,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -4,11 +4,11 @@\n \n const LIKELY_ERROR_CAUSES = \"\"\"\n 1. Bugs in the program that generated the sample, or otherwise incorrect input variables. \n 2. Your chains failed to converge. Check your diagnostics. \n-3. You do not have enough posterior samples (Less than ~100 samples) -- try sampling more values.\n+3. You do not have enough posterior samples (ESS < ~100).\n \"\"\"\n-const MIN_TAIL_LEN = 16  # Minimum size of a tail for PSIS to give sensible answers\n+const MIN_TAIL_LEN = 10  # Minimum size of a tail for PSIS to give sensible answers\n const SAMPLE_SOURCES = [\"mcmc\", \"vi\", \"other\"]\n \n export Psis, psis\n \n@@ -251,9 +251,9 @@\n     log_ratios::AbstractArray{T,3}, r_eff::AbstractVector{T}\n ) where {T<:AbstractFloat}\n     if any(isnan, log_ratios)\n         throw(DomainError(\"Invalid input for `log_ratios` (contains NaN values).\"))\n-    elseif any(isinf, log_ratios)\n+    elseif any(isinf, log_rati    os)\n         throw(DomainError(\"Invalid input for `log_ratios` (contains infinite values).\"))\n     elseif isempty(log_ratios)\n         throw(ArgumentError(\"Invalid input for `log_ratios` (array is empty).\"))\n     elseif any(isnan, r_eff)\n"
+                },
+                {
+                    "date": 1626484737810,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -6,9 +6,9 @@\n 1. Bugs in the program that generated the sample, or otherwise incorrect input variables. \n 2. Your chains failed to converge. Check your diagnostics. \n 3. You do not have enough posterior samples (ESS < ~100).\n \"\"\"\n-const MIN_TAIL_LEN = 10  # Minimum size of a tail for PSIS to give sensible answers\n+const MIN_TAIL_LEN = 5  # Minimum size of a tail for PSIS to give sensible answers\n const SAMPLE_SOURCES = [\"mcmc\", \"vi\", \"other\"]\n \n export Psis, psis\n \n@@ -251,9 +251,9 @@\n     log_ratios::AbstractArray{T,3}, r_eff::AbstractVector{T}\n ) where {T<:AbstractFloat}\n     if any(isnan, log_ratios)\n         throw(DomainError(\"Invalid input for `log_ratios` (contains NaN values).\"))\n-    elseif any(isinf, log_rati    os)\n+    elseif any(isinf, log_ratios)\n         throw(DomainError(\"Invalid input for `log_ratios` (contains infinite values).\"))\n     elseif isempty(log_ratios)\n         throw(ArgumentError(\"Invalid input for `log_ratios` (array is empty).\"))\n     elseif any(isnan, r_eff)\n"
+                },
+                {
+                    "date": 1626484881897,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -292,9 +292,9 @@\n end\n \n \n \"\"\"\n-Assume that all objects belong to a single chain if chain index is missing. Warn user.\n+Assume that all objects belong to a single chain if chain index is missing. Inform user.\n \"\"\"\n function _assume_one_chain(log_ratios)\n     @info \"Chain information was not provided; \" *\n           \"all samples are assumed to be drawn from a single chain.\"\n"
+                },
+                {
+                    "date": 1626486225465,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -81,10 +81,10 @@\n     end\n \n     if any(ξ .≥ .7)\n         @warn \"Some Pareto k values are very high (>0.7), indicating that PSIS has \" * \n-        \"failed to approximate the true distribution for these points. Treat them \" *\n-        \"with caution.\"\n+        \"failed to approximate the true distribution for these points. Treat these \" *\n+        \"estimates with caution.\"\n     elseif any(ξ .≥ .5)\n         @info \"Some Pareto k values are slightly high (>0.5); convergence may be slow \" *\n         \"and MCSE estimates may be slight underestimates.\"\n     end\n"
+                },
+                {
+                    "date": 1626486967591,
+                    "content": "Index: \n===================================================================\n--- \n+++ \n@@ -17,9 +17,9 @@\n         log_ratios::AbstractArray{T:>AbstractFloat}, \n         r_eff; \n         source::String=\"mcmc\", \n         log_weights::Bool=false\n-        ) -> Psis\n+    ) -> Psis\n \n Implements Pareto-smoothed importance sampling (PSIS).\n \n # Arguments\n"
                 }
             ],
             "date": 1626309654640,