Test for approx equality with R

Author: Closed-Limelike-Curves

Project.toml

@@ -4,6 +4,7 @@ authors = ["Carlos Parada <[email protected]>"]
 version = "0.1.1"
 
 [deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"

src/GPD.jl

@@ -55,7 +55,7 @@ function gpdfit(
     # build pointwise estimates of ξ and θ at each grid point
     θ_hats = similar(sample, m)
     ξ_hats = similar(sample, m)
-    @turbo @. θ_hats = 1 / sample[len] + (1 - sqrt((m+1) / $(1:m))) / prior / quartile
+    @turbo @. θ_hats = 1 / sample[len] + (1 - sqrt(m / ($(1:m)-.5))) / prior / quartile
     @tullio threads=false ξ_hats[x] := log1p(-θ_hats[x] * sample[y]) |> _ / len
     log_like = similar(ξ_hats)
     # Calculate profile log-likelihood at each estimate:

src/ImportanceSampling.jl

@@ -64,7 +64,7 @@ function psis(
     weights::AbstractArray{F} = similar(log_ratios)
     # Shift ratios by maximum to prevent overflow
     # then shift by log of posterior sample size to avoid loss of precision for small values
-    @tturbo @. weights = exp(log_ratios + $log(post_sample_size) - $maximum(log_ratios; dims=2))
+    @tturbo @. weights = exp(log_ratios - $maximum(log_ratios; dims=2))
 
     rel_eff = generate_rel_eff(weights, dims, rel_eff, source)
     check_input_validity_psis(reshape(log_ratios, dims), rel_eff)
@@ -85,10 +85,6 @@ function psis(
         @tturbo @. weights = log(weights)
     end
 
-    if dims[3] == 1
-        weights = dropdims(weights; dims=3)  # Reshape as array
-    end
-
     return Psis(
         weights,
         ξ,
@@ -138,13 +134,12 @@ Do PSIS on a single vector, smoothing its tail values.
 
 # Arguments
 
-  - `is_ratios::AbstractVector{AbstractFloat}`: A vector of (not necessarily
-  normalized) importance sampling ratios.
+- `is_ratios::AbstractVector{AbstractFloat}`: A vector of importance sampling ratios, 
+scaled to have a maximum of 1.
 
 # Returns
 
-  - `T<:AbstractFloat`: ξ, the shape parameter for the GPD; larger numbers indicate
-  thicker tails.
+- `T<:AbstractFloat`: ξ, the shape parameter for the GPD; big numbers indicate thick tails.
 
 # Extended help
 
@@ -169,8 +164,8 @@ function do_psis_i!(
     cutoff = sorted_ratios[tail_start-1]
     ξ = psis_smooth_tail!(tail, cutoff)
 
-    # truncate at max of raw wts; because is_ratios ∝ len / maximum, max(raw weights) = len
-    clamp!(sorted_ratios, 0, len)
+    # truncate at max of raw wts (which is 1)
+    clamp!(sorted_ratios, 0, 1)
     # unsort the ratios to their original position:
     is_ratios .= @views sorted_ratios[invperm(ordering)]
 
@@ -202,7 +197,7 @@ function psis_smooth_tail!(tail::AbstractVector{T}, cutoff::T) where {T<:Abstrac
     # save time not sorting since tail is already sorted
     ξ, σ = GPD.gpdfit(tail)
     if ξ ≠ Inf
-        @turbo @. tail = GPD.gpd_quantile($(1:len) / (len + 1), ξ, σ) + cutoff
+        @turbo @. tail = GPD.gpd_quantile(($(1:len) - .5) / len, ξ, σ) + cutoff
     end
     return ξ
 end
@@ -229,7 +224,7 @@ observations).
 """
 struct Psis{
     F<:AbstractFloat,
-    AF<:AbstractArray{F},
+    AF<:AbstractArray{F,3},
     VF<:AbstractVector{F},
     I<:Integer,
     VI<:AbstractVector{I},

test/Manifest.toml

@@ -60,12 +60,6 @@ git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
 uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
 version = "0.1.1"
 
-[[BenchmarkTools]]
-deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
-git-tree-sha1 = "01ca3823217f474243cc2c8e6e1d1f45956fe872"
-uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.0.0"
-
 [[BufferedStreams]]
 deps = ["Compat", "Test"]
 git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f"
@@ -409,12 +403,6 @@ version = "1.1.0"
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
-[[PkgBenchmark]]
-deps = ["BenchmarkTools", "Dates", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Pkg", "Printf", "TerminalLoggers", "UUIDs"]
-git-tree-sha1 = "e4a10b7cdb7ec836850e43a4cee196f4e7b02756"
-uuid = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
-version = "0.2.12"
-
 [[PooledArrays]]
 deps = ["DataAPI", "Future"]
 git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"

test/Project.toml

@@ -1,6 +1,5 @@
 [deps]
 Libz = "2ec943e9-cfe8-584d-b93d-64dcb6d567b7"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
-PkgBenchmark = "32113eaa-f34f-5b0d-bd6c-c81e245fc73d"
 RData = "df47a6cb-8c03-5eed-afd8-b6050d6c41da"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

test/runtests.jl

@@ -24,6 +24,7 @@ logPsis = psis(logLikelihoodArray; lw=true)
     @test sqrt(mean((with_rel_eff.weights ./ rWeights .- 1).^2)) ≤ .001
     # RMSE less than .2% when using InferenceDiagnostics' ESS
     @test sqrt(mean((juliaPsis.weights ./ rWeights .- 1).^2)) ≤ .002
-    @test count(juliaPsis.weights .≉ matrixPsis.weights) ≤ 1
+    @test count(with_rel_eff.weights .≈ rWeights) ≤ 10
+    @test count(juliaPsis.weights .≉ matrixPsis.weights) ≤ 10
     @test sqrt(mean((logPsis.weights .- log.(rWeights)).^2)) ≤ .001
 end