see description

- Updated method argument types - Changed correlation types to singleton types resembling Base.RoundingMode - In place methods now return the mutated object - Irrational constants now provided by IrrationalConstants.jl - Root finding function now type stable and has better interface
SchisslerGroup · Feb 19, 2024 · 89609ac · 89609ac
1 parent 3e6fb2b
commit 89609ac
Show file tree

Hide file tree

Showing 14 changed files with 225 additions and 213 deletions.
diff --git a/Project.toml b/Project.toml
@@ -10,6 +10,7 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
@@ -26,20 +27,21 @@ Aqua = "0.8"
 Distributions = "0.22, 0.23, 0.24, 0.25"
 FastGaussQuadrature = "0.4, 0.5, 1"
 HypergeometricFunctions = "0.2, 0.3"
+IrrationalConstants = "0.2"
 IntervalArithmetic = "0.16, 0.17, 0.18, 0.19, 0.20"
 IntervalRootFinding = "0.5"
 IterTools = "1.3, 1.4"
-LsqFit = "0.10, 0.12, 0.13, 0.14, 0.15"
 LinearAlgebra = "1"
+LsqFit = "0.10, 0.12, 0.13, 0.14, 0.15"
 Polynomials = "1, 2, 3"
+PrecompileTools = "1"
 QuadGK = "2.4, 2.6"
+SharedArrays = "1"
 SpecialFunctions = "1, 2"
+Statistics = "1"
 StatsBase = "0.32, 0.33, 0.34"
-julia = "1.6, 1.7, 1.8, 1.9, 1.10"
 Test = "1"
-PrecompileTools = "1"
-Statistics = "1"
-SharedArrays = "1"
+julia = "1.6, 1.7, 1.8, 1.9, 1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

diff --git a/src/Bigsimr.jl b/src/Bigsimr.jl
@@ -4,6 +4,7 @@ using Base.Threads: @threads
 using Distributions
 using FastGaussQuadrature: gausshermite
 using HypergeometricFunctions: _₂F₁
+using IrrationalConstants
 using IntervalArithmetic: interval, mid
 using IntervalRootFinding: roots, Krawczyk
 using IterTools: subsets
@@ -14,54 +15,53 @@ using QuadGK: quadgk
 using SharedArrays
 using SpecialFunctions: erfc, erfcinv
 using StatsBase: corspearman, corkendall
+using Statistics: clampcor
 
 import Distributions: mean, std, quantile, cdf, pdf, var, params
-import LinearAlgebra: diag, inv, logdet
-import Statistics: cor, clampcor
+import Statistics
 
 
 struct ValidCorrelationError <: Exception end
 
 
-abstract type Correlation end
+"""
+    CorType
+
+A type used for specifiying the type of correlation. Supported correlations are:
+
+    - [`Pearson`](@ref)
+    - [`Spearman`](@ref)
+    - [`Kendall`](@ref)
+"""
+struct CorType{T} end
 
 """
-    Pearson <: Correlation
+    Pearson
 
 Pearson's ``r`` product-moment correlation
 """
-struct Pearson <: Correlation end
+const Pearson = CorType{:Pearson}()
 
 """
-    Spearman <: Correlation
+    Spearman
 
 Spearman's ``ρ`` rank correlation
 """
-struct Spearman <: Correlation end
+const Spearman = CorType{:Spearman}()
 
 """
-    Kendall <: Correlation
+    Kendall
 
 Kendall's ``τ`` rank correlation
 """
-struct Kendall <: Correlation end
+const Kendall = CorType{:Kendall}()
+
 
 
 const UD  = UnivariateDistribution
 const CUD = ContinuousUnivariateDistribution
 const DUD = DiscreteUnivariateDistribution
 
-const sqrt2_f64      ::Float64 = sqrt(2)
-const invsqrt2_f64   ::Float64 = inv(sqrt(2))
-const invsqrtpi_f64  ::Float64 = inv(sqrt(π))
-const invsqrt2pi_f64 ::Float64 = inv(sqrt(2π))
-
-const sqrt2_f32      ::Float32 = sqrt2_f64
-const invsqrt2_f32   ::Float32 = invsqrt2_f64
-const invsqrtpi_f32  ::Float32 = invsqrtpi_f64
-const invsqrt2pi_f32 ::Float32 = invsqrt2pi_f64
-
-
 
 include("utils.jl")
 include("random_vector.jl")
@@ -80,35 +80,35 @@ include("GSDist.jl")
 
 
 
-export 
+export
     # correlation types
-    Correlation, 
-    Pearson, 
-    Spearman, 
+    CorType,
+    Pearson,
+    Spearman,
     Kendall,
     # random vector generation
-    rvec, 
+    rvec,
     rmvn,
     # correlation calculation
     cor,
     cor_fast,
-    cor_convert, 
-    cor_bounds, 
+    cor_convert,
+    cor_bounds,
     # nearest correlation
-    cor_nearPD, 
-    cor_fastPD, 
+    cor_nearPD,
+    cor_fastPD,
     cor_fastPD!,
     # random correlation generation
-    cor_randPD, 
+    cor_randPD,
     cor_randPSD,
     # correlation Utils
     iscorrelation,
-    cor_constrain, 
+    cor_constrain,
     cor_constrain!,
     cov2cor,
     cov2cor!,
     # pearson methods
-    pearson_match, 
+    pearson_match,
     pearson_bounds
 
 

diff --git a/src/cor_bounds.jl b/src/cor_bounds.jl
@@ -1,5 +1,5 @@
 """
-    cor_bounds(d1::UnivariateDistribution, d2::UnivariateDistribution, C::Type{<:Correlation}=Pearson; n_samples::Real=100_000)
+    cor_bounds(d1::UnivariateDistribution, d2::UnivariateDistribution, cortype::CorType=Pearson; n_samples=100_000)
 
 Compute the stochastic lower and upper correlation bounds between two marginal
 distributions.
@@ -35,29 +35,29 @@ julia> cor_bounds(A, B, Spearman)
 (lower = -1.0, upper = 1.0)
 ```
 """
-function cor_bounds(d1::UD, d2::UD, C::Type{<:Correlation}=Pearson; n_samples=100_000)
-    return _cor_bounds(d1, d2, C, Int(n_samples))
+function cor_bounds(d1::UD, d2::UD, cortype::CorType=Pearson; n_samples=100_000)
+    return _cor_bounds(d1, d2, cortype, Int(n_samples))
 end
 
-function _cor_bounds(d1, d2, C, n)
+function _cor_bounds(d1::UD, d2::UD, cortype::CorType, n::Int)
     a = rand(d1, n)
     b = rand(d2, n)
 
     sort!(a)
     sort!(b)
-    
-    upper = cor(a, b, C)
+
+    upper = cor(a, b, cortype)
 
     reverse!(b)
-    lower = cor(a, b, C)
+    lower = cor(a, b, cortype)
 
     (lower = lower, upper = upper)
 end
 
 
 
 """
-    cor_bounds(margins::AbstractVector{<:UD}, C::Type{<:Correlation}=Pearson; n_samples::Real=100_000)
+    cor_bounds(margins::AbstractVector{<:UnivariateDistribution}, cortype::CorType=Pearson; n_samples=100_000)
 
 Compute the stochastic pairwise lower and upper correlation bounds between a set
 of marginal distributions.
@@ -67,25 +67,22 @@ The possible correlation types are:
   * [`Spearman`](@ref)
   * [`Kendall`](@ref)
 """
-function cor_bounds(margins::AbstractVector{<:UD}, C::Type{<:Correlation}=Pearson; n_samples=100_000)
-    return _cor_bounds(margins, C, Int(n_samples))
+function cor_bounds(margins::AbstractVector{<:UD}, cortype::CorType=Pearson; n_samples=100_000)
+    return _cor_bounds(margins, cortype, Int(n_samples))
 end
 
-function _cor_bounds(margins, C, n)
+function _cor_bounds(margins::AbstractVector{<:UD}, cortype::CorType, n::Int)
     d = length(margins)
     lower, upper = zeros(Float64, d, d), zeros(Float64, d, d)
 
     @threads for i in collect(subsets(1:d, Val{2}()))
-        l, u = cor_bounds(margins[i[1]], margins[i[2]], C; n_samples=n)
+        l, u = cor_bounds(margins[i[1]], margins[i[2]], cortype; n_samples=n)
         lower[i...] = l
         upper[i...] = u
     end
 
-    lower .= Symmetric(lower)
     cor_constrain!(lower)
-
-    upper .= Symmetric(upper)
     cor_constrain!(upper)
-    
+
     (lower = lower, upper = upper)
 end
diff --git a/src/cor_fastPD.jl b/src/cor_fastPD.jl
@@ -1,15 +1,12 @@
-function _fast_pca!(X, λ, P, n)
+function _fast_pca!(X::AbstractMatrix{T}, λ::AbstractVector{T}, P::AbstractMatrix{T}, n::Int) where {T}
     r = sum(λ .> 0)
-
     r == n && return X
-
-    T = eltype(X)
     s = n - r
 
     r == 0 && return fill!(X, zero(T))
 
-    if r == 1 
-        X .= (λ[1] * λ[1]) * (P[:,1] * P[:,1]')   
+    if r == 1
+        X .= (λ[1] * λ[1]) * (P[:,1] * P[:,1]')
         return X
     end
 
@@ -31,16 +28,16 @@ end
 
 
 """
-    cor_fastPD!(R, τ=1e-6)
+    cor_fastPD!(R::AbstractMatrix{<:Real}, tau=1e-6)
 
-Same as [`cor_fastPD`](@ref), but saves space by overwriting the input `R`
-instead of creating a copy.
+Same as [`cor_fastPD`](@ref), but saves space by overwriting the input `R` instead of
+creating a copy.
 
 See also: [`cor_fastPD`](@ref), [`cor_nearPD`](@ref)
 
 # Examples
 ```jldoctest
-julia> import LinearAlgebra: isposdef
+julia> using LinearAlgebra: isposdef
 
 julia> r = [1.00 0.82 0.56 0.44; 0.82 1.00 0.28 0.85; 0.56 0.28 1.00 0.22; 0.44 0.85 0.22 1.00]
 4×4 Matrix{Float64}:
@@ -53,17 +50,20 @@ julia> isposdef(r)
 false
 
 julia> cor_fastPD!(r)
-
+4×4 Matrix{Float64}:
+ 1.0       0.817095  0.559306  0.440514
+ 0.817095  1.0       0.280196  0.847352
+ 0.559306  0.280196  1.0       0.219582
+ 0.440514  0.847352  0.219582  1.0
 
 julia> isposdef(r)
 true
 ```
 """
-function cor_fastPD!(R::AbstractMatrix{<:Real}, tau=1e-6)
-    T = eltype(R)
+function cor_fastPD!(R::AbstractMatrix{T}, tau=1e-6) where {T<:Real}
     n  = size(R, 1)
     tau  = max(eps(T), tau)
-    
+
     R .= Symmetric(R, :U)
     R[diagind(R)] .= (one(T) - tau)
 
@@ -74,9 +74,7 @@ function cor_fastPD!(R::AbstractMatrix{<:Real}, tau=1e-6)
     _fast_pca!(R, λ, P, n)
 
     R[diagind(R)] .+= tau
-    R .= cov2cor(R)
-
-    return R
+    return cov2cor!(R)
 end