Add custom warning for k greater than 1 (#2)

* Add Printf as dependency

* Customize warning for Pareto k greater than 1

* Increment version number

Author: sethaxen

Project.toml

@@ -1,10 +1,11 @@
 name = "PSIS"
 uuid = "ce719bf2-d5d0-4fb9-925d-10a81b42ad04"
 authors = ["Seth Axen <[email protected]> and contributors"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
@@ -14,8 +15,8 @@ julia = "1"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]

src/PSIS.jl

@@ -2,6 +2,7 @@ module PSIS
 
 using Statistics: mean
 using LinearAlgebra: dot
+using Printf: @sprintf
 
 export psis, psis!
 
@@ -88,17 +89,15 @@ function psis!(logw, r_eff=1.0; sorted=issorted(logw), normalize=false)
         @inbounds logw_tail = @views logw[perm[tail_range]]
         if logw_max - first(logw_tail) < eps(eltype(logw_tail)) / 100
             @warn "Cannot fit the generalized Pareto distribution because all tail " *
-                  "values are the same"
+                "values are the same"
         else
             logw_tail .-= logw_max
             @inbounds logu = logw[perm[icut]] - logw_max
 
             _, k_hat = psis_tail!(logw_tail, logu, M)
             logw_tail .+= logw_max
 
-            k_hat ≥ 0.7 &&
-                @warn "Pareto k statistic exceeded 0.7. Resulting importance sampling estimates " *
-                      "are likely to be unstable."
+            check_pareto_k(k_hat)
         end
     end
 
@@ -109,6 +108,18 @@ function psis!(logw, r_eff=1.0; sorted=issorted(logw), normalize=false)
     return logw, k_hat
 end
 
+function check_pareto_k(k)
+    if k ≥ 1
+        @warn "Pareto k=$(@sprintf("%.2g", k)) ≥ 1. Resulting importance sampling " *
+            "estimates are likely to be unstable and are unlikely to converge with " *
+            "additional samples."
+    elseif k ≥ 0.7
+        @warn "Pareto k=$(@sprintf("%.2g", k)) ≥ 0.7. Resulting importance sampling " *
+            "estimates are likely to be unstable."
+    end
+    return nothing
+end
+
 tail_length(r_eff, S) = min(cld(S, 5), ceil(Int, 3 * sqrt(S / r_eff)))
 
 function psis_tail!(logw, logu, M=length(logw))

test/psis.jl

@@ -95,9 +95,35 @@ end
         @test k > 0.7
         msg = String(take!(io))
         @test occursin(
-            "Warning: Pareto k statistic exceeded 0.7. Resulting importance sampling estimates are likely to be unstable",
+            "Resulting importance sampling estimates are likely to be unstable", msg
+        )
+
+        io = IOBuffer()
+        with_logger(SimpleLogger(io)) do
+            PSIS.check_pareto_k(1.1)
+        end
+        msg = String(take!(io))
+        @test occursin(
+            "Warning: Pareto k=1.1 ≥ 1. Resulting importance sampling estimates are likely to be unstable and are unlikely to converge with additional samples.",
             msg,
         )
+
+        io = IOBuffer()
+        with_logger(SimpleLogger(io)) do
+            PSIS.check_pareto_k(0.8)
+        end
+        msg = String(take!(io))
+        @test occursin(
+            "Warning: Pareto k=0.8 ≥ 0.7. Resulting importance sampling estimates are likely to be unstable.",
+            msg,
+        )
+
+        io = IOBuffer()
+        with_logger(SimpleLogger(io)) do
+            PSIS.check_pareto_k(0.69)
+        end
+        msg = String(take!(io))
+        @test isempty(msg)
     end
 
     has_loo() && @testset "consistent with loo" begin