diff --git a/Project.toml b/Project.toml index 1352a11..05078da 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TransformVariables" uuid = "84d833dd-6860-57f9-a1a7-6da5db126cff" authors = ["Tamas K. Papp "] -version = "0.8.8" +version = "0.8.9" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" @@ -11,7 +11,6 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -22,7 +21,9 @@ ChangesOfVariables = "0.1" DocStringExtensions = "0.8, 0.9" ForwardDiff = "0.10" InverseFunctions = "0.1" +LinearAlgebra = "1.6" LogExpFunctions = "0.3" +Random = "1.6" SimpleUnPack = "1" StaticArrays = "1" julia = "1.6" diff --git a/src/special_arrays.jl b/src/special_arrays.jl index 360752e..dc19350 100644 --- a/src/special_arrays.jl +++ b/src/special_arrays.jl @@ -4,21 +4,31 @@ export UnitVector, UnitSimplex, CorrCholeskyFactor, corr_cholesky_factor #### building blocks #### +""" +$(SIGNATURES) + +`log(abs(…))` of the derivative of `tanh`, calculated accurately. +""" +function _tanh_logabsderiv(x) + d = 2*x + log(4) + d - 2 * log1pexp(d) +end + """ (y, r, ℓ) = $SIGNATURES Given ``x ∈ ℝ`` and ``0 ≤ r ≤ 1``, return `(y, r′)` such that -1. ``y² + r′² = r²``, +1. ``y² + (r′)² = r²``, 2. ``y: |y| ≤ r`` is mapped with a bijection from `x`. `ℓ` is the log Jacobian (whether it is evaluated depends on `flag`). """ @inline function l2_remainder_transform(flag::LogJacFlag, x, r) - z = 2*logistic(x) - 1 - (z * √r, r*(1 - abs2(z)), - flag isa NoLogJac ? flag : log(2) + logistic_logjac(x) + 0.5*log(r)) + # note that 1-tanh(x)^2 = sech(x)^2 + (tanh(x) * √r, r*sech(x)^2, + flag isa NoLogJac ? flag : _tanh_logabsderiv(x) + 0.5*log(r)) end """ @@ -26,7 +36,7 @@ end Inverse of [`l2_remainder_transform`](@ref) in `x` and `y`. """ -@inline l2_remainder_inverse(y, r) = logit((y/√r+1)/2), r-abs2(y) +@inline l2_remainder_inverse(y, r) = atanh(y/√r), r-y^2 #### #### UnitVector diff --git a/test/Project.toml b/test/Project.toml index cc42be0..a987591 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,6 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c" LogDensityProblemsAD = "996a588d-648d-4e1f-a8f0-a84b347e47b1" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/test/runtests.jl b/test/runtests.jl index 0824947..9c80888 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -659,3 +659,15 @@ end y.b[3] += 1 @test x[4] == z + 1 end + +@testset "near-singular Cholesky factor" begin + x = [8.348500225024523, -3.80486310849193, -15.115725300837742, 5.840812234057503, + 7.548980701857334, 0.6546495312434718, 2.863837638357627, 1.0081703617568052, + -38.543769810398466, -14.252165683848483, -22.75952203884357, 1.9543987098768612, + 5.415229912144962, -1.4360948924991273, 4.957606068283541, -5.443369115798325, + -2.536087079311158, -2.0710241403850635, -0.982209305513312, 6.821758239096414, + 5.925173901833287] + t = corr_cholesky_factor(7) + U = transform(t, x) + @test isfinite(logabsdet(U)[1]) +end