diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ec726e7..4c1a9dd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,9 +1,6 @@ name: CI on: pull_request: - branches: - - main - - dev push: branches: - main @@ -17,6 +14,7 @@ jobs: fail-fast: false matrix: version: + - '1.6' - '1.4' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'. - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. - 'nightly' diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index 40ebdca..0000000 --- a/.travis.yml +++ /dev/null @@ -1,39 +0,0 @@ -language: julia -os: - - linux - # - osx -# arch: arm64-graviton2 - -julia: - # - 1.2 - - 1.4 - - nightly - -notifications: - email: false - -after_success: - # - julia -e 'import Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder()); Coveralls.submit(process_folder())' - -branches: - only: - - main - - /^v\d+\.\d+(\.\d+)?(-\S*)?$/ - - dev - -jobs: - allow_failures: - - julia: nightly - fast_finish: true - include: - - stage: "Documentation" - julia: 1.4 - os: linux - script: - # - julia --project=perf/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); - # Pkg.instantiate()' - # - julia --project=perf/ perf/samplers.jl - - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); - Pkg.instantiate()' - - julia --project=docs/ docs/make.jl - after_success: skip \ No newline at end of file diff --git a/Project.toml b/Project.toml index 3ac2ab4..38a4ecb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LRMoE" uuid = "b6847444-3cfa-4637-857f-9c8c615f8351" authors = ["Spark Tseung "] -version = "0.3.1" +version = "0.3.2" [deps] Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5" diff --git a/appveyor.yml b/appveyor.yml deleted file mode 100644 index 2877a93..0000000 --- a/appveyor.yml +++ /dev/null @@ -1,42 +0,0 @@ -environment: - matrix: - - julia_version: 1 - - julia_version: nightly - -platform: - - x86 - - x64 - -# # Uncomment the following lines to allow failures on nightly julia -# # (tests will run but not make your overall status red) -matrix: - allow_failures: - - julia_version: nightly - -branches: - only: - - main - - /release-.*/ - -notifications: - - provider: Email - on_build_success: false - on_build_failure: false - on_build_status_changed: false - -install: - - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1")) - -build_script: - - echo "%JL_BUILD_SCRIPT%" - - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%" - -test_script: - - echo "%JL_TEST_SCRIPT%" - - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%" - -# # Uncomment to support code coverage upload. Should only be enabled for packages -# # which would have coverage gaps without running on Windows -# on_success: -# - echo "%JL_CODECOV_SCRIPT%" -# - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%" \ No newline at end of file diff --git a/docs/src/dummytest.md b/docs/src/dummytest.md deleted file mode 100644 index d7c77d2..0000000 --- a/docs/src/dummytest.md +++ /dev/null @@ -1,3 +0,0 @@ -# Dummy - -dummy \ No newline at end of file diff --git a/src/LRMoE.jl b/src/LRMoE.jl index f058e3a..7ea091a 100644 --- a/src/LRMoE.jl +++ b/src/LRMoE.jl @@ -79,15 +79,10 @@ export # loglikelihood functions pdf, logpdf, cdf, logcdf, - rowlogsumexp, expert_ll_pos, - expert_tn_pos, - expert_tn_bar_pos, + rowlogsumexp, expert_ll, expert_tn, expert_tn_bar, - expert_ll_pos_list, - expert_tn_pos_list, - expert_tn_bar_pos_list, expert_ll_list, expert_tn_list, expert_tn_bar_list, loglik_aggre_dim, diff --git a/src/experts/continuous/burr.jl b/src/experts/continuous/burr.jl index 0ade617..053fc51 100644 --- a/src/experts/continuous/burr.jl +++ b/src/experts/continuous/burr.jl @@ -134,9 +134,14 @@ mean(d::BurrExpert) = mean(LRMoE.Burr(d.k, d.c, d.λ)) var(d::BurrExpert) = var(LRMoE.Burr(d.k, d.c, d.λ)) quantile(d::BurrExpert, p) = quantile(LRMoE.Burr(d.k, d.c, d.λ), p) function lev(d::BurrExpert, u) - uu = 1 / (1 + (u / d.λ)^d.c) - return d.λ * gamma(float(1 + 1 / d.c)) * gamma(float(d.k - 1 / d.c)) / - gamma(float(d.k)) * beta_inc(1 + 1 / d.c, d.k - 1 / d.c, 1 - uu)[1] + u * uu^d.k + if isinf(u) + return mean(d) + else + uu = 1 / (1 + (u / d.λ)^d.c) + return d.λ * gamma(float(1 + 1 / d.c)) * gamma(float(d.k - 1 / d.c)) / + gamma(float(d.k)) * beta_inc(1 + 1 / d.c, d.k - 1 / d.c, 1 - uu)[1] + + u * uu^d.k + end end excess(d::BurrExpert, u) = mean(d) - lev(d, u) diff --git a/src/experts/continuous/gamma.jl b/src/experts/continuous/gamma.jl index 2427358..ee9410e 100644 --- a/src/experts/continuous/gamma.jl +++ b/src/experts/continuous/gamma.jl @@ -25,7 +25,7 @@ end ## Outer constructors GammaExpert(k::Real, θ::Real) = GammaExpert(promote(k, θ)...) GammaExpert(k::Integer, θ::Integer) = GammaExpert(float(k), float(θ)) -GammaExpert() = GammaExpert(1.0, 1.0) +GammaExpert() = GammaExpert(2.0, 1.0) ## Conversion function convert(::Type{GammaExpert{T}}, k::S, θ::S) where {T<:Real,S<:Real} @@ -135,8 +135,12 @@ mean(d::GammaExpert) = mean(Distributions.Gamma(d.k, d.θ)) var(d::GammaExpert) = var(Distributions.Gamma(d.k, d.θ)) quantile(d::GammaExpert, p) = quantile(Distributions.Gamma(d.k, d.θ), p) function lev(d::GammaExpert, u) - return d.θ * d.k * gamma_inc(float(d.k + 1), u / d.θ, 0)[1] + - u * (1 - gamma_inc(float(d.k), u / d.θ, 0)[1]) + if isinf(u) + return mean(d) + else + return d.θ * d.k * gamma_inc(float(d.k + 1), u / d.θ, 0)[1] + + u * (1 - gamma_inc(float(d.k), u / d.θ, 0)[1]) + end end excess(d::GammaExpert, u) = mean(d) - lev(d, u) diff --git a/src/experts/continuous/inversegaussian.jl b/src/experts/continuous/inversegaussian.jl index cdbaa31..2679dfd 100644 --- a/src/experts/continuous/inversegaussian.jl +++ b/src/experts/continuous/inversegaussian.jl @@ -141,10 +141,14 @@ mean(d::InverseGaussianExpert) = mean(Distributions.InverseGaussian(d.μ, d.λ)) var(d::InverseGaussianExpert) = var(Distributions.InverseGaussian(d.μ, d.λ)) quantile(d::InverseGaussianExpert, p) = quantile(Distributions.InverseGaussian(d.μ, d.λ), p) function lev(d::InverseGaussianExpert, u) - z = (u - d.μ) / d.μ - y = (u + d.μ) / d.μ - return u - d.μ * z * cdf.(Normal(0, 1), z * (d.λ / u)^0.5) - - d.μ * y * exp(2 * d.λ / d.μ) * cdf.(Normal(0, 1), -y * (d.λ / u)^0.5) + if isinf(u) + return mean(d) + else + z = (u - d.μ) / d.μ + y = (u + d.μ) / d.μ + return u - d.μ * z * cdf.(Normal(0, 1), z * (d.λ / u)^0.5) - + d.μ * y * exp(2 * d.λ / d.μ) * cdf.(Normal(0, 1), -y * (d.λ / u)^0.5) + end end excess(d::InverseGaussianExpert, u) = mean(d) - lev(d, u) diff --git a/src/experts/continuous/lognormal.jl b/src/experts/continuous/lognormal.jl index f00e345..caba230 100644 --- a/src/experts/continuous/lognormal.jl +++ b/src/experts/continuous/lognormal.jl @@ -111,8 +111,12 @@ mean(d::LogNormalExpert) = mean(Distributions.LogNormal(d.μ, d.σ)) var(d::LogNormalExpert) = var(Distributions.LogNormal(d.μ, d.σ)) quantile(d::LogNormalExpert, p) = quantile(Distributions.LogNormal(d.μ, d.σ), p) function lev(d::LogNormalExpert, u) - return exp(d.μ + 0.5 * d.σ^2) * cdf.(Normal(d.μ + d.σ^2, d.σ), log(u)) + - u * (1 - cdf.(Normal(d.μ, d.σ), log(u))) + if isinf(u) + return mean(d) + else + return exp(d.μ + 0.5 * d.σ^2) * cdf.(Normal(d.μ + d.σ^2, d.σ), log(u)) + + u * (1 - cdf.(Normal(d.μ, d.σ), log(u))) + end end excess(d::LogNormalExpert, u) = mean(d) - lev(d, u) @@ -131,7 +135,7 @@ function _int_obs_logY(d::LogNormalExpert, yl, yu, expert_ll) return log(yl) else return exp(-expert_ll) * ( - d.σ * invsqrt2π * 0.5 * _diff_dens_series(d, yl, yu) + + d.σ * invsqrt2π * _diff_dens_series(d, yl, yu) + d.μ * _diff_dist_series(d, yl, yu) ) end @@ -140,7 +144,7 @@ end function _int_lat_logY(d::LogNormalExpert, tl, tu, expert_tn_bar) return exp(-expert_tn_bar) * ( d.μ - ( - d.σ * invsqrt2π * 0.5 * _diff_dens_series(d, tl, tu) + + d.σ * invsqrt2π * _diff_dens_series(d, tl, tu) + d.μ * _diff_dist_series(d, tl, tu) ) ) diff --git a/src/experts/continuous/weibull.jl b/src/experts/continuous/weibull.jl index 7f50130..178d712 100644 --- a/src/experts/continuous/weibull.jl +++ b/src/experts/continuous/weibull.jl @@ -129,8 +129,12 @@ mean(d::WeibullExpert) = mean(Distributions.Weibull(d.k, d.θ)) var(d::WeibullExpert) = var(Distributions.Weibull(d.k, d.θ)) quantile(d::WeibullExpert, p) = quantile(Distributions.Weibull(d.k, d.θ), p) function lev(d::WeibullExpert, u) - return d.θ * gamma_inc(float(1 / d.k + 1), (u / d.θ)^d.k, 0)[1] * - gamma(float(1 / d.k + 1)) + u * exp(-(u / d.θ)^d.k) + if isinf(u) + return mean(d) + else + return d.θ * gamma_inc(float(1 / d.k + 1), (u / d.θ)^d.k, 0)[1] * + gamma(float(1 / d.k + 1)) + u * exp(-(u / d.θ)^d.k) + end end excess(d::WeibullExpert, u) = mean(d) - lev(d, u) diff --git a/test/fit_functions.jl b/test/archive/fit_functions.jl similarity index 100% rename from test/fit_functions.jl rename to test/archive/fit_functions.jl diff --git a/test/fit_functions_new.jl b/test/archive/fit_functions_new.jl similarity index 100% rename from test/fit_functions_new.jl rename to test/archive/fit_functions_new.jl diff --git a/test/paramsinit.jl b/test/archive/paramsinit.jl similarity index 100% rename from test/paramsinit.jl rename to test/archive/paramsinit.jl diff --git a/test/predict.jl b/test/archive/predict.jl similarity index 100% rename from test/predict.jl rename to test/archive/predict.jl diff --git a/test/real_init.jl b/test/archive/real_init.jl similarity index 100% rename from test/real_init.jl rename to test/archive/real_init.jl diff --git a/test/dummytest.jl b/test/dummytest.jl deleted file mode 100644 index 6edeaf4..0000000 --- a/test/dummytest.jl +++ /dev/null @@ -1,10 +0,0 @@ -using Test - -@testset "dummytest" begin - - @test 1.0+1.0 == 2.0 - @test 2.0+2.0 == 4.0 - @test 3.0+3.0 == 6.0 - # @test 4.0+4.0 == 8.0 - -end \ No newline at end of file diff --git a/test/expert_loglik.jl b/test/expert_loglik.jl deleted file mode 100644 index 3335ae5..0000000 --- a/test/expert_loglik.jl +++ /dev/null @@ -1,827 +0,0 @@ -using Test -using Distributions -using StatsFuns - -using LRMoE - -@testset "count with exposures" begin - - # mean(LRMoE.GammaCountExpert(1.0, 1.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1.0)) - # mean(LRMoE.GammaCountExpert(3.0, 1.0)) - # mean(LRMoE.GammaCountExpert(4.0, 1.0)) - # mean(LRMoE.GammaCountExpert(5.0, 1.0)) - - # mean(LRMoE.GammaCountExpert(1.0, 2.0)) - # mean(LRMoE.GammaCountExpert(2.0, 2.0)) - # mean(LRMoE.GammaCountExpert(3.0, 2.0)) - # mean(LRMoE.GammaCountExpert(4.0, 2.0)) - # mean(LRMoE.GammaCountExpert(5.0, 2.0)) - - # mean(LRMoE.GammaCountExpert(1.0, 3.0)) - # mean(LRMoE.GammaCountExpert(2.0, 3.0)) - # mean(LRMoE.GammaCountExpert(3.0, 3.0)) - # mean(LRMoE.GammaCountExpert(4.0, 3.0)) - # mean(LRMoE.GammaCountExpert(5.0, 3.0)) - - # mean(LRMoE.GammaCountExpert(1.0, 1/1.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/2.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/3.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/4.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/5.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/6.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/7.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/8.0)) - # mean(LRMoE.GammaCountExpert(1.0, 1/9.0)) - - # var(LRMoE.GammaCountExpert(1.0, 1/1.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/2.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/3.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/4.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/5.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/6.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/7.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/8.0)) - # var(LRMoE.GammaCountExpert(1.0, 1/9.0)) - - # mean(LRMoE.GammaCountExpert(2.0, 1/1.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/2.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/3.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/4.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/5.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/6.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/7.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/8.0)) - # mean(LRMoE.GammaCountExpert(2.0, 1/9.0)) - - # var(LRMoE.GammaCountExpert(2.0, 1/1.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/2.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/3.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/4.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/5.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/6.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/7.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/8.0)) - # var(LRMoE.GammaCountExpert(2.0, 1/9.0)) - - # mean(LRMoE.NegativeBinomialExpert(1.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(2.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(3.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(4.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(5.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(6.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(7.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(8.0, 0.20)) - # mean(LRMoE.NegativeBinomialExpert(9.0, 0.20)) - - # var(LRMoE.NegativeBinomialExpert(1.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(2.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(3.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(4.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(5.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(6.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(7.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(8.0, 0.20)) - # var(LRMoE.NegativeBinomialExpert(9.0, 0.20)) - -end - - - - -@testset "expert_ll: continuous, NonZI" begin - - # param1 = 3.0 - # param2 = 2.0 - # param3 = 5.0 - - # l = Distributions.LogNormal(param1, param2) - # x = rand(l, 20) - - # tmp = LRMoE.BurrExpert(param1, param2, param3) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.GammaExpert(param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.InverseGaussianExpert(param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.LogNormalExpert(param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.WeibullExpert(param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), -Inf) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), 0.0) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - -end - - - -@testset "expert_ll: continuous, ZI" begin - - # param0 = rand(Distributions.Uniform(0, 1), 1)[1] - # param1 = 3.0 - # param2 = 2.0 - # param3 = 5.0 - - # l = Distributions.LogNormal(param1, param2) - # x = rand(l, 20) .+ 0.05 - - # tmp = LRMoE.ZIBurrExpert(param0, param1, param2, param3) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0)) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6) - - - # tmp = LRMoE.ZIGammaExpert(param0, param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0)) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6) - - - # tmp = LRMoE.ZIInverseGaussianExpert(param0, param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0)) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6) - - - # tmp = LRMoE.ZILogNormalExpert(param0, param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0)) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6) - - - # tmp = LRMoE.ZIWeibullExpert(param0, param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll_exact.(tmp, 0), log(param0)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.((1-param0).*exp.(LRMoE.logpdf.(tmp, x)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log.(param0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.((1-param0).*exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)))), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), log.(1), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.((1-param0).*exp.(logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, 0.75 .*x) .- logcdf.(tmp, 1.25 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.((1-param0).*exp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x)))), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log(1-param0)) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(param0 .+ (1-param0).*exp.(log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, 0.25 .*x) .- logcdf.(tmp, 1.5 .*x))))), atol = 1e-6) - -end - - - -@testset "expert_ll: discrete, NonZI" begin - - # param1 = 20 - # param2 = rand(Distributions.Uniform(0, 1), 1)[1] - # param3 = 5.0 - - # l = Distributions.Binomial(param1, param2) - # x = rand(l, 20) - - # tmp = LRMoE.BinomialExpert(param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # tmp = LRMoE.GammaCountExpert(param1, param3) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1] - # # tmp_expo = LRMoE.GammaCountExpert(param1, param3/expo) - - # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), LRMoE.logpdf.(tmp_expo, x)) - - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), LRMoE.logpdf.(tmp_expo, x)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0.0)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0)) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6) - - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x))) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(LRMoE.logpdf.(tmp_expo, 0)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.NegativeBinomialExpert(param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1] - # # tmp_expo = LRMoE.NegativeBinomialExpert(param1*expo, param2) - - # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), LRMoE.logpdf.(tmp_expo, x)) - - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), LRMoE.logpdf.(tmp_expo, x)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0.0)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0)) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6) - - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x))) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(LRMoE.logpdf.(tmp_expo, 0)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.PoissonExpert(param1) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0.0)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), LRMoE.logpdf.(tmp, 0)) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(LRMoE.logpdf.(tmp, 0)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1] - # # tmp_expo = LRMoE.PoissonExpert(param1*expo) - - # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), LRMoE.logpdf.(tmp_expo, x)) - - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), LRMoE.logpdf.(tmp_expo, x)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0.0)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), LRMoE.logpdf.(tmp_expo, 0)) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6) - - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x))) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(LRMoE.logpdf.(tmp_expo, 0)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6) - -end - - - -@testset "expert_ll: discrete, ZI" begin - - # param0 = rand(Distributions.Uniform(0, 1), 1)[1] - # param1 = 20 - # param2 = rand(Distributions.Uniform(0, 1), 1)[1] - # param3 = 5.0 - - # l = Distributions.Binomial(param1, param2) - # x = rand(l, 20) .+ 1 - - # tmp = LRMoE.ZIBinomialExpert(param0, param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # tmp = LRMoE.ZIGammaCountExpert(param0, param1, param3) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1] - # # tmp_expo = LRMoE.ZIGammaCountExpert(param0, param1, param3/expo) - - # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x)) - - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6) - - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x))) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))), atol = 1e-06) - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.ZINegativeBinomialExpert(param0, param1, param2) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1] - # # tmp_expo = LRMoE.ZINegativeBinomialExpert(param0, param1*expo, param2) - - # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x)) - - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6) - - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x))) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))), atol = 1e-06) - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6) - - - # tmp = LRMoE.ZIPoissonExpert(param0, param1) - - # @test isapprox(LRMoE.expert_ll_exact.(tmp, x), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), log.(1-param0) .+ LRMoE.logpdf.(tmp, x)) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), log.(1-param0) .+ logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf), logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, -1) .- logcdf.(tmp, Inf)), atol = 1e-6) - # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf), log.(1-param0) .+ logcdf.(tmp, 1.25 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp, 1.25 .*x)), atol = 1e-6) - - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x)), atol = 1e-06) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))) - # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log.(1-param0) .+ logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) - # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp, 0.0)))), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x), log1mexp.(logcdf.(tmp, 1.5 .*x) .+ log1mexp.(logcdf.(tmp, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp, 1.5 .*x))), atol = 1e-6) - - # # expo = rand(Distributions.Uniform(0.5, 5), 1)[1] - # # tmp_expo = LRMoE.ZIPoissonExpert(param0, param1*expo) - - # # @test isapprox(LRMoE.expert_ll_exact.(tmp, x, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x)) - - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf, exposure = expo), log.(1-param0) .+ LRMoE.logpdf.(tmp_expo, x)) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, ceil.(x) .- 1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0, Inf, Inf, exposure = expo), logcdf.(tmp_expo, Inf) .+ log1mexp.(logcdf.(tmp_expo, -1) .- logcdf.(tmp_expo, Inf)), atol = 1e-6) - # # @test isapprox(LRMoE.expert_ll.(tmp, 0.0, 0.75 .*x, 1.25 .*x, Inf, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.25 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.75 .*x) .- 1) .- logcdf.(tmp_expo, 1.25 .*x)), atol = 1e-6) - - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf, exposure = expo), fill(0.0, length(x)), atol = 1e-06) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))) - # # @test isapprox(LRMoE.expert_tn.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log.(1-param0) .+ logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x)), atol = 1e-6) - - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf, exposure = expo), fill(-Inf, length(x))) - # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, 0.0, 0.0, 0.0, exposure = expo), log1mexp(log(param0 + (1-param0) .* exp.(LRMoE.logpdf.(tmp_expo, 0.0)))), atol = 1e-06) - # # # @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.25 .*x, 0.75 .*x, 1.25 .*x, 1.5 .*x, exposure = expo), log1mexp.(logcdf.(tmp_expo, 1.5 .*x) .+ log1mexp.(logcdf.(tmp_expo, ceil.(0.25 .*x) .- 1) .- logcdf.(tmp_expo, 1.5 .*x))), atol = 1e-6) - -end - - - -@testset "exposurize_expert" begin - - # param0 = rand(Distributions.Uniform(0, 1), 1)[1] - # param1 = rand(Distributions.Uniform(10, 20), 1)[1] - # param2 = rand(Distributions.Uniform(0, 1), 1)[1] - # param3 = rand(Distributions.Uniform(5, 10), 1)[1] - - # expo = rand(Distributions.Uniform(0.5, 20), 1)[1] - - # tmp = LRMoE.BurrExpert(param1, param2, param3) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - # tmp = LRMoE.ZIBurrExpert(param0, param1, param2, param3) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - - # tmp = LRMoE.GammaExpert(param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - # tmp = LRMoE.ZIGammaExpert(param0, param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - - # tmp = LRMoE.InverseGaussianExpert(param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - # tmp = LRMoE.ZIInverseGaussianExpert(param0, param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - - # tmp = LRMoE.LogNormalExpert(param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - # tmp = LRMoE.ZILogNormalExpert(param0, param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - - # tmp = LRMoE.WeibullExpert(param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - # tmp = LRMoE.ZIWeibullExpert(param0, param1, param2) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - - - # tmp = LRMoE.BinomialExpert(20, param0) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - # tmp = LRMoE.ZIBinomialExpert(0.20, 20, param0) - # @test tmp == LRMoE.exposurize_expert(tmp, exposure = expo) - - - # tmp = LRMoE.GammaCountExpert(param1, param2) - # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo) - # @test isapprox([LRMoE.params(tmp)...] .* [1, 1/expo], [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...], atol = 1e-6) - # tmp = LRMoE.ZIGammaCountExpert(param0, param1, param2) - # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo) - # @test isapprox([LRMoE.params(tmp)...] .* [1, 1, 1/expo], [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...], atol = 1e-6) - - # tmp = LRMoE.NegativeBinomialExpert(param1, param0) - # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo) - # @test [LRMoE.params(tmp)...] .* [expo, 1] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...] - # tmp = LRMoE.ZINegativeBinomialExpert(param0, param1, param0) - # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo) - # @test [LRMoE.params(tmp)...] .* [1, expo, 1] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...] - - # tmp = LRMoE.PoissonExpert(param1) - # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo) - # @test [LRMoE.params(tmp)...] .* [expo] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...] - # tmp = LRMoE.ZIPoissonExpert(param0, param1) - # @test tmp != LRMoE.exposurize_expert(tmp, exposure = expo) - # @test [LRMoE.params(tmp)...] .* [1, expo] == [LRMoE.params(LRMoE.exposurize_expert(tmp, exposure = expo))...] - - - - - # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3); - # BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)] - - # expo = rand(Distributions.Uniform(0.5, 20), 15) - - # @test size(LRMoE.exposurize_model(model, exposure = expo)) == (2, 4, 15) - - - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,1,:] == fill(mean(model[1,1]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,2,:] == fill(mean(model[1,2]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,3,:] == fill(mean(model[1,3]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,4,:] == fill(mean(model[1,4]), 15) - - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,1,:] == fill(mean(model[2,1]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,2,:] == mean(model[2,2]) .* expo - - - # model_expod = LRMoE.exposurize_model(model, exposure = expo) - # Y_tmp = hcat(rand(Distributions.LogNormal(2, 1), 15), rand(Distributions.Poisson(5), 15)) - - # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[1,:], model_expod[:,:,1])) == (2, 4) - # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[5,:], model_expod[:,:,5])) == (2, 4) - # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[12,:], model_expod[:,:,12])) == (2, 4) - # @test size(LRMoE.expert_ll_list_exact(Y_tmp, model_expod)) == (2, 4, 15) - - # idx = rand(1:15, 1)[1] - # mat_tmp = LRMoE.expert_ll_ind_mat_exact(Y_tmp[idx,:], model_expod[:,:,idx]) - # @test mat_tmp[1,1] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1]) - # @test mat_tmp[1,2] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1]) - # @test mat_tmp[1,3] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1]) - # @test mat_tmp[1,4] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1]) - # @test mat_tmp[2,1] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2]) - # @test mat_tmp[2,2] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2]) - # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2]) - # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2]) - - # idx = rand(1:15, 1)[1] - # cube_tmp = LRMoE.expert_ll_list_exact(Y_tmp, model_expod) - # @test cube_tmp[1,1,idx] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1]) - # @test cube_tmp[1,2,idx] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1]) - # @test cube_tmp[1,3,idx] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1]) - # @test cube_tmp[1,4,idx] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1]) - # @test cube_tmp[2,1,idx] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2]) - # @test cube_tmp[2,2,idx] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2]) - # @test cube_tmp[2,3,idx] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2]) - # @test cube_tmp[2,4,idx] == LRMoE.expert_ll_exact(model_expod[2,4,idx], Y_tmp[idx,2]) - - -end - - - -@testset "exposurize_model" begin - # param0 = rand(Distributions.Uniform(0, 1), 1)[1] - # param1 = rand(Distributions.Uniform(10, 20), 1)[1] - # param2 = rand(Distributions.Uniform(0, 1), 1)[1] - # param3 = rand(Distributions.Uniform(5, 10), 1)[1] - - # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3); - # BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)] - - # expo = rand(Distributions.Uniform(0.5, 20), 15) - - # @test size(LRMoE.exposurize_model(model, exposure = expo)) == (2, 4, 15) - - - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,1,:] == fill(mean(model[1,1]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,2,:] == fill(mean(model[1,2]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,3,:] == fill(mean(model[1,3]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[1,4,:] == fill(mean(model[1,4]), 15) - - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,1,:] == fill(mean(model[2,1]), 15) - # @test mean.(LRMoE.exposurize_model(model, exposure = expo))[2,2,:] == mean(model[2,2]) .* expo - - - # model_expod = LRMoE.exposurize_model(model, exposure = expo) - # Y_tmp = hcat(rand(Distributions.LogNormal(2, 1), 15), rand(Distributions.Poisson(5), 15)) - - # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[1,:], model_expod[:,:,1])) == (2, 4) - # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[5,:], model_expod[:,:,5])) == (2, 4) - # @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[12,:], model_expod[:,:,12])) == (2, 4) - # @test size(LRMoE.expert_ll_list_exact(Y_tmp, model_expod)) == (2, 4, 15) - - # idx = rand(1:15, 1)[1] - # mat_tmp = LRMoE.expert_ll_ind_mat_exact(Y_tmp[idx,:], model_expod[:,:,idx]) - # @test mat_tmp[1,1] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1]) - # @test mat_tmp[1,2] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1]) - # @test mat_tmp[1,3] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1]) - # @test mat_tmp[1,4] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1]) - # @test mat_tmp[2,1] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2]) - # @test mat_tmp[2,2] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2]) - # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2]) - # @test mat_tmp[2,3] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2]) - - # idx = rand(1:15, 1)[1] - # cube_tmp = LRMoE.expert_ll_list_exact(Y_tmp, model_expod) - # @test cube_tmp[1,1,idx] == LRMoE.expert_ll_exact(model_expod[1,1,idx], Y_tmp[idx,1]) - # @test cube_tmp[1,2,idx] == LRMoE.expert_ll_exact(model_expod[1,2,idx], Y_tmp[idx,1]) - # @test cube_tmp[1,3,idx] == LRMoE.expert_ll_exact(model_expod[1,3,idx], Y_tmp[idx,1]) - # @test cube_tmp[1,4,idx] == LRMoE.expert_ll_exact(model_expod[1,4,idx], Y_tmp[idx,1]) - # @test cube_tmp[2,1,idx] == LRMoE.expert_ll_exact(model_expod[2,1,idx], Y_tmp[idx,2]) - # @test cube_tmp[2,2,idx] == LRMoE.expert_ll_exact(model_expod[2,2,idx], Y_tmp[idx,2]) - # @test cube_tmp[2,3,idx] == LRMoE.expert_ll_exact(model_expod[2,3,idx], Y_tmp[idx,2]) - # @test cube_tmp[2,4,idx] == LRMoE.expert_ll_exact(model_expod[2,4,idx], Y_tmp[idx,2]) - - - # Y_tmp = hcat(0.2 .* Y_tmp[:,1], 0.8 .* Y_tmp[:,1], 1.2 .* Y_tmp[:,1], 5 .* Y_tmp[:,1], fill(0, length(Y_tmp[:,2])), Y_tmp[:,2], Y_tmp[:,2], fill(Inf, length(Y_tmp[:,2]))) - - # @test size(LRMoE.expert_ll_ind_mat(Y_tmp[1,:], model_expod[:,:,1])) == (2, 4) - # @test size(LRMoE.expert_ll_ind_mat(Y_tmp[5,:], model_expod[:,:,5])) == (2, 4) - # @test size(LRMoE.expert_ll_ind_mat(Y_tmp[12,:], model_expod[:,:,12])) == (2, 4) - # @test size(LRMoE.expert_ll_list(Y_tmp, model_expod)) == (2, 4, 15) - - # idx = rand(1:15, 1)[1] - # mat_tmp = LRMoE.expert_ll_ind_mat(Y_tmp[idx,:], model_expod[:,:,idx]) - # @test mat_tmp[1,1] == LRMoE.expert_ll(model_expod[1,1,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test mat_tmp[1,2] == LRMoE.expert_ll(model_expod[1,2,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test mat_tmp[1,3] == LRMoE.expert_ll(model_expod[1,3,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test mat_tmp[1,4] == LRMoE.expert_ll(model_expod[1,4,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test mat_tmp[2,1] == LRMoE.expert_ll(model_expod[2,1,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - # @test mat_tmp[2,2] == LRMoE.expert_ll(model_expod[2,2,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - # @test mat_tmp[2,3] == LRMoE.expert_ll(model_expod[2,3,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - # @test mat_tmp[2,3] == LRMoE.expert_ll(model_expod[2,3,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - - # idx = rand(1:15, 1)[1] - # cube_tmp = LRMoE.expert_ll_list(Y_tmp, model_expod) - # @test cube_tmp[1,1,idx] == LRMoE.expert_ll(model_expod[1,1,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test cube_tmp[1,2,idx] == LRMoE.expert_ll(model_expod[1,2,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test cube_tmp[1,3,idx] == LRMoE.expert_ll(model_expod[1,3,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test cube_tmp[1,4,idx] == LRMoE.expert_ll(model_expod[1,4,idx], Y_tmp[idx,1], Y_tmp[idx,2], Y_tmp[idx,3], Y_tmp[idx,4]) - # @test cube_tmp[2,1,idx] == LRMoE.expert_ll(model_expod[2,1,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - # @test cube_tmp[2,2,idx] == LRMoE.expert_ll(model_expod[2,2,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - # @test cube_tmp[2,3,idx] == LRMoE.expert_ll(model_expod[2,3,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - # @test cube_tmp[2,4,idx] == LRMoE.expert_ll(model_expod[2,4,idx], Y_tmp[idx,5], Y_tmp[idx,6], Y_tmp[idx,7], Y_tmp[idx,8]) - - # # ll related - # gate_tmp = rand(Distributions.Uniform(-1, 1), 15, 4) - - # dim_agg = LRMoE.loglik_aggre_dim(cube_tmp) - # @test size(dim_agg) == (15, 4) - # @test size(LRMoE.loglik_aggre_gate_dim(gate_tmp, dim_agg)) == (15, 4) - - - # # Simulation related - # # X_tmp = rand(Distributions.Uniform(-1, 1), 15, 7) - # # α_tmp = rand(Distributions.Uniform(-1, 1), 4, 7) - - # # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3); - # # BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)] - - # # expo = rand(Distributions.Uniform(0.5, 20), 15) - - # # model_expod = LRMoE.exposurize_model(model, exposure = expo) - - # # sim_tmp = LRMoE.sim_dataset(α_tmp, X_tmp, model, exposure = expo) - # # mean(sim_tmp, dims = 1) - - # # sim_tmp = LRMoE.sim_dataset(α_tmp, X_tmp, model) - # # mean(sim_tmp, dims = 1) - - -end \ No newline at end of file diff --git a/test/experts/continuous/burr.jl b/test/experts/continuous/burr.jl new file mode 100644 index 0000000..4e2895a --- /dev/null +++ b/test/experts/continuous/burr.jl @@ -0,0 +1,96 @@ +@testset "burr pdf, cdf, etc." begin + kk = [1.5, 2.0, 5.0, 10.0] + cc = [1.5, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for k in kk, c in cc, λ in λλ + l = LRMoE.Burr(k, c, λ) + r = LRMoE.BurrExpert(k, c, λ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "burr expert_ll" begin + kk = [1.5, 2.0, 5.0, 10.0] + cc = [1.5, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for k in kk, c in cc, λ in λλ + tmp = LRMoE.BurrExpert(k, c, λ) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0.0] .= -Inf + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0.0] .= -Inf + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0.0] .= 0.0 + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "burr exposurize" begin + d = BurrExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "burr lev/excess" begin + kk = [1.5, 2.0, 5.0, 10.0] + cc = [1.5, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for k in kk, c in cc, λ in λλ + tmp = BurrExpert(k, c, λ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == mean(Burr(k, c, λ)) + @test LRMoE.lev(tmp, Inf) == mean(Burr(k, c, λ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end \ No newline at end of file diff --git a/test/experts/continuous/gamma.jl b/test/experts/continuous/gamma.jl new file mode 100644 index 0000000..db9babf --- /dev/null +++ b/test/experts/continuous/gamma.jl @@ -0,0 +1,94 @@ +# Won't test: k <= 1.0 due to spurious likelihood at zero +@testset "gamma pdf, cdf, etc." begin + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for k in kk, θ in θθ + l = Distributions.Gamma(k, θ) + r = LRMoE.GammaExpert(k, θ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "gamma expert_ll" begin + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for k in kk, θ in θθ + tmp = LRMoE.GammaExpert(k, θ) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0.0] .= -Inf + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0.0] .= -Inf + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0.0] .= 0.0 + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "gamma exposurize" begin + d = GammaExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "gamma lev/excess" begin + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for k in kk, θ in θθ + tmp = GammaExpert(k, θ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == mean(Gamma(k, θ)) + @test LRMoE.lev(tmp, Inf) == mean(Gamma(k, θ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end \ No newline at end of file diff --git a/test/experts/continuous/inversegaussian.jl b/test/experts/continuous/inversegaussian.jl new file mode 100644 index 0000000..9c1e56b --- /dev/null +++ b/test/experts/continuous/inversegaussian.jl @@ -0,0 +1,93 @@ +@testset "inversegaussian pdf, cdf, etc." begin + μμ = [1.0, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for μ in μμ, λ in λλ + l = Distributions.InverseGaussian(μ, λ) + r = LRMoE.InverseGaussianExpert(μ, λ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "inversegaussian expert_ll" begin + μμ = [1.0, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for μ in μμ, λ in λλ + tmp = LRMoE.InverseGaussianExpert(μ, λ) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0.0] .= -Inf # x == 0.0 + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0.0] .= -Inf # x == 0.0 + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0.0] .= 0.0 # x == 0.0 + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "inversegaussian exposurize" begin + d = InverseGaussianExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "inversegaussian lev/excess" begin + μμ = [1.0, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for μ in μμ, λ in λλ + tmp = InverseGaussianExpert(μ, λ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == mean(InverseGaussian(μ, λ)) + @test LRMoE.lev(tmp, Inf) == mean(InverseGaussian(μ, λ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end \ No newline at end of file diff --git a/test/experts/continuous/lognormal.jl b/test/experts/continuous/lognormal.jl new file mode 100644 index 0000000..ee02092 --- /dev/null +++ b/test/experts/continuous/lognormal.jl @@ -0,0 +1,124 @@ +@testset "lognormal pdf, cdf, etc." begin + μμ = [1.0, 2.0, 5.0, 10.0] + σσ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for μ in μμ, σ in σσ + l = Distributions.LogNormal(μ, σ) + r = LRMoE.LogNormalExpert(μ, σ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "lognormal expert_ll" begin + μμ = [1.0, 2.0, 5.0, 10.0] + σσ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for μ in μμ, σ in σσ + tmp = LRMoE.LogNormalExpert(μ, σ) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0.0] .= -Inf # x == 0.0 + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0.0] .= -Inf # x == 0.0 + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0.0] .= 0.0 # x == 0.0 + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "lognormal exposurize" begin + d = LogNormalExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "lognormal lev/excess" begin + μμ = [1.0, 2.0, 5.0, 10.0] + σσ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for μ in μμ, σ in σσ + tmp = LogNormalExpert(μ, σ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == mean(LogNormal(μ, σ)) + @test LRMoE.lev(tmp, Inf) == mean(LogNormal(μ, σ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end + +# TODO: this is currently not working +# @testset "lognormal numerical integration" begin +# μμ = [1.0, 2.0, 5.0, 10.0] +# σσ = [1.0, 0.5, 0.1, 0.01] +# yl = [0.5, 1.0] # , 2.0, 5.0, 10.0, 20.0, 50.0, 100.0] +# yu_multiplier = [1.5] # [1.0, 1.5, 2.0, 5.0] + +# for μ in μμ, σ in σσ, m in yu_multiplier +# tmp = LogNormalExpert(μ, σ) +# d = LogNormal(μ, σ) +# yu = m .* yl +# expert_ll = +# LRMoE.expert_ll.(tmp, fill(0.0, length(yl)), yl, yu, fill(Inf, length(yl))) + +# fn_integrate_logY(d, y) = log(y) * pdf(d, y) + +# @test isapprox(LRMoE._int_obs_logY.(tmp, yl, yu, expert_ll), +# [ +# if yl_i == yu_i +# log(yl_i) +# else +# quadgk.(x -> fn_integrate_logY(d, x), yl_i, yu_i; rtol=1e-8)[1] * +# exp(-1.0 * e_ll) +# end +# for +# (yl_i, yu_i, e_ll) in zip(yl, yu, expert_ll) +# ], +# atol=1e-6) +# end +# end \ No newline at end of file diff --git a/test/experts/continuous/weibull.jl b/test/experts/continuous/weibull.jl new file mode 100644 index 0000000..00b0c3b --- /dev/null +++ b/test/experts/continuous/weibull.jl @@ -0,0 +1,94 @@ +# Won't test: k <= 1.0 due to spurious likelihood at zero +@testset "weibull pdf, cdf, etc." begin + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for k in kk, θ in θθ + l = Distributions.Weibull(k, θ) + r = LRMoE.WeibullExpert(k, θ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "weibull expert_ll" begin + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for k in kk, θ in θθ + tmp = LRMoE.WeibullExpert(k, θ) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0.0] .= -Inf + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0.0] .= -Inf + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0.0] .= 0.0 + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "weibull exposurize" begin + d = WeibullExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "weibull lev/excess" begin + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for k in kk, θ in θθ + tmp = WeibullExpert(k, θ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == mean(Weibull(k, θ)) + @test LRMoE.lev(tmp, Inf) == mean(Weibull(k, θ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end \ No newline at end of file diff --git a/test/experts/continuous/ziburr.jl b/test/experts/continuous/ziburr.jl new file mode 100644 index 0000000..b93d4a1 --- /dev/null +++ b/test/experts/continuous/ziburr.jl @@ -0,0 +1,130 @@ +@testset "zi_burr pdf, cdf, etc." begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + cc = [1.5, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, k in kk, c in cc, λ in λλ + l = LRMoE.Burr(k, c, λ) + r = LRMoE.ZIBurrExpert(p, k, c, λ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zi_burr expert_ll" begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + cc = [1.5, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, k in kk, c in cc, λ in λλ + tmp = LRMoE.ZIBurrExpert(p, k, c, λ) + + # exact observations + true_val = log(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0.0] .= log(p) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val, atol=1e-6) + + # right censoring + true_val = + log.( + (1 - p) .* + exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf))) + ) + true_val[x .== 0.0] .= + log.( + p .+ + (1 - p) .* + exp.( + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0.0) .- logcdf.(tmp, Inf)) + ) + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + true_val = fill(0.0, length(x)) + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + true_val = fill(-Inf, length(x)) + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + p .+ + (1 - p) .* + exp.( + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + ) + true_val[x .== 0.0] .= log(1 - p) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "zi_burr exposurize" begin + d = ZIBurrExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "zi_burr lev/excess" begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + cc = [1.5, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for p in pp, k in kk, c in cc, λ in λλ + tmp = ZIBurrExpert(p, k, c, λ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == (1 - p) * mean(Burr(k, c, λ)) + @test LRMoE.lev(tmp, Inf) == (1 - p) * mean(Burr(k, c, λ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end diff --git a/test/experts/continuous/zigamma.jl b/test/experts/continuous/zigamma.jl new file mode 100644 index 0000000..e0c8647 --- /dev/null +++ b/test/experts/continuous/zigamma.jl @@ -0,0 +1,127 @@ +@testset "zigamma pdf, cdf, etc." begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, k in kk, θ in θθ + l = Distributions.Gamma(k, θ) + r = LRMoE.ZIGammaExpert(p, k, θ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zigamma expert_ll" begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, k in kk, θ in θθ + tmp = LRMoE.ZIGammaExpert(p, k, θ) + + # exact observations + true_val = log(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0.0] .= log(p) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val, atol=1e-6) + + # right censoring + true_val = + log.( + (1 - p) .* + exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf))) + ) + true_val[x .== 0.0] .= + log.( + p .+ + (1 - p) .* + exp.( + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0.0) .- logcdf.(tmp, Inf)) + ) + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + true_val = fill(0.0, length(x)) + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + true_val = fill(-Inf, length(x)) + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + p .+ + (1 - p) .* + exp.( + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + ) + true_val[x .== 0.0] .= log(1 - p) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "zigamma exposurize" begin + d = ZIGammaExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "zigamma lev/excess" begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for p in pp, k in kk, θ in θθ + tmp = ZIGammaExpert(p, k, θ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == (1 - p) * mean(Gamma(k, θ)) + @test LRMoE.lev(tmp, Inf) == (1 - p) * mean(Gamma(k, θ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end diff --git a/test/experts/continuous/ziinversegaussian.jl b/test/experts/continuous/ziinversegaussian.jl new file mode 100644 index 0000000..3833de0 --- /dev/null +++ b/test/experts/continuous/ziinversegaussian.jl @@ -0,0 +1,127 @@ +@testset "ziinversegaussian pdf, cdf, etc." begin + pp = [0.05, 0.10, 0.50, 0.90] + μμ = [1.0, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, μ in μμ, λ in λλ + l = Distributions.InverseGaussian(μ, λ) + r = LRMoE.ZIInverseGaussianExpert(p, μ, λ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "ziinversegaussian expert_ll" begin + pp = [0.05, 0.10, 0.50, 0.90] + μμ = [1.0, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, μ in μμ, λ in λλ + tmp = LRMoE.ZIInverseGaussianExpert(p, μ, λ) + + # exact observations + true_val = log(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0.0] .= log(p) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val, atol=1e-6) + + # right censoring + true_val = + log.( + (1 - p) .* + exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf))) + ) + true_val[x .== 0.0] .= + log.( + p .+ + (1 - p) .* + exp.( + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0.0) .- logcdf.(tmp, Inf)) + ) + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + true_val = fill(0.0, length(x)) + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + true_val = fill(-Inf, length(x)) + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + p .+ + (1 - p) .* + exp.( + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + ) + true_val[x .== 0.0] .= log(1 - p) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "ziinversegaussian exposurize" begin + d = ZIInverseGaussianExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "ziinversegaussian lev/excess" begin + pp = [0.05, 0.10, 0.50, 0.90] + μμ = [1.0, 2.0, 5.0, 10.0] + λλ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for p in pp, μ in μμ, λ in λλ + tmp = ZIInverseGaussianExpert(p, μ, λ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == (1 - p) * mean(InverseGaussian(μ, λ)) + @test LRMoE.lev(tmp, Inf) == (1 - p) * mean(InverseGaussian(μ, λ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end diff --git a/test/experts/continuous/zilognormal.jl b/test/experts/continuous/zilognormal.jl new file mode 100644 index 0000000..e5d1d93 --- /dev/null +++ b/test/experts/continuous/zilognormal.jl @@ -0,0 +1,127 @@ +@testset "zilognormal pdf, cdf, etc." begin + pp = [0.05, 0.10, 0.50, 0.90] + μμ = [1.0, 2.0, 5.0, 10.0] + σσ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, μ in μμ, σ in σσ + l = Distributions.LogNormal(μ, σ) + r = LRMoE.ZILogNormalExpert(p, μ, σ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zilognormal expert_ll" begin + pp = [0.05, 0.10, 0.50, 0.90] + μμ = [1.0, 2.0, 5.0, 10.0] + σσ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, μ in μμ, σ in σσ + tmp = LRMoE.ZILogNormalExpert(p, μ, σ) + + # exact observations + true_val = log(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0.0] .= log(p) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val, atol=1e-6) + + # right censoring + true_val = + log.( + (1 - p) .* + exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf))) + ) + true_val[x .== 0.0] .= + log.( + p .+ + (1 - p) .* + exp.( + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0.0) .- logcdf.(tmp, Inf)) + ) + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + true_val = fill(0.0, length(x)) + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + true_val = fill(-Inf, length(x)) + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + p .+ + (1 - p) .* + exp.( + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + ) + true_val[x .== 0.0] .= log(1 - p) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "zilognormal exposurize" begin + d = ZILogNormalExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "zilognormal lev/excess" begin + pp = [0.05, 0.10, 0.50, 0.90] + μμ = [1.0, 2.0, 5.0, 10.0] + σσ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for p in pp, μ in μμ, σ in σσ + tmp = ZILogNormalExpert(p, μ, σ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == (1 - p) * mean(LogNormal(μ, σ)) + @test LRMoE.lev(tmp, Inf) == (1 - p) * mean(LogNormal(μ, σ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end diff --git a/test/experts/continuous/ziweibull.jl b/test/experts/continuous/ziweibull.jl new file mode 100644 index 0000000..8135580 --- /dev/null +++ b/test/experts/continuous/ziweibull.jl @@ -0,0 +1,127 @@ +@testset "ziweibullpdf, cdf, etc." begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, k in kk, θ in θθ + l = Distributions.Weibull(k, θ) + r = LRMoE.ZIWeibullExpert(p, k, θ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "ziweibullexpert_ll" begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:0.5:100.0) + + for p in pp, k in kk, θ in θθ + tmp = LRMoE.ZIWeibullExpert(p, k, θ) + + # exact observations + true_val = log(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0.0] .= log(p) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val, atol=1e-6) + + # right censoring + true_val = + log.( + (1 - p) .* + exp.(logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, x) .- logcdf.(tmp, Inf))) + ) + true_val[x .== 0.0] .= + log.( + p .+ + (1 - p) .* + exp.( + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, 0.0) .- logcdf.(tmp, Inf)) + ) + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, 0.75 .* x) .- logcdf.(tmp, 1.25 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + true_val = fill(0.0, length(x)) + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + (1 - p) .* + exp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + true_val[x .== 0.0] .= log(p) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + true_val = fill(-Inf, length(x)) + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), true_val, atol=1e-6) + + # censoring + truncation + true_val = + log.( + p .+ + (1 - p) .* + exp.( + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, 0.25 .* x) .- logcdf.(tmp, 1.5 .* x)) + ) + ) + ) + true_val[x .== 0.0] .= log(1 - p) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "ziweibullexposurize" begin + d = ZIWeibullExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end + +@testset "ziweibulllev/excess" begin + pp = [0.05, 0.10, 0.50, 0.90] + kk = [1.5, 2.0, 5.0, 10.0] + θθ = [1.0, 0.5, 0.1, 0.01] + x = [0.0, Inf] + + for p in pp, k in kk, θ in θθ + tmp = ZIWeibullExpert(p, k, θ) + @test LRMoE.lev(tmp, 0.0) == 0.0 + @test LRMoE.excess(tmp, 0.0) == (1 - p) * mean(Weibull(k, θ)) + @test LRMoE.lev(tmp, Inf) == (1 - p) * mean(Weibull(k, θ)) + @test LRMoE.excess(tmp, Inf) == 0.0 + end +end diff --git a/test/experts/discrete/binomial.jl b/test/experts/discrete/binomial.jl new file mode 100644 index 0000000..e1b0757 --- /dev/null +++ b/test/experts/discrete/binomial.jl @@ -0,0 +1,77 @@ +@testset "binomial pdf, cdf, etc." begin + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for n in nn, p in pp + l = Distributions.Binomial(n, p) + r = LRMoE.BinomialExpert(n, p) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "binomial expert_ll" begin + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for n in nn, p in pp + tmp = LRMoE.BinomialExpert(n, p) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "binomial exposurize" begin + d = BinomialExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end diff --git a/test/experts/discrete/gammacount.jl b/test/experts/discrete/gammacount.jl new file mode 100644 index 0000000..d6662d5 --- /dev/null +++ b/test/experts/discrete/gammacount.jl @@ -0,0 +1,77 @@ +@testset "gammacount pdf, cdf, etc." begin + mm = [1.5, 2, 5, 10, 25] + ss = [1.0, 2.0, 5.0, 10.0, 25.0] + x = collect(0.0:1:50.0) + + for m in mm, s in ss + l = LRMoE.GammaCount(m, s) + r = LRMoE.GammaCountExpert(m, s) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "gammacount expert_ll" begin + mm = [1.5, 2, 5, 10, 25] + ss = [1.0, 2.0, 5.0, 10.0, 25.0] + x = collect(0.0:1:50.0) + + for m in mm, s in ss + tmp = LRMoE.GammaCountExpert(m, s) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +# @testset "gammacount exposurize" begin +# d = GammaCountExpert() +# exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + +# for e in exposure_vec +# @test mean(LRMoE.exposurize_expert(d; exposure=e)) == mean(d) * e +# end +# end diff --git a/test/experts/discrete/negativebinomial.jl b/test/experts/discrete/negativebinomial.jl new file mode 100644 index 0000000..282715d --- /dev/null +++ b/test/experts/discrete/negativebinomial.jl @@ -0,0 +1,77 @@ +@testset "negativebinomial pdf, cdf, etc." begin + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for n in nn, p in pp + l = Distributions.NegativeBinomial(n, p) + r = LRMoE.NegativeBinomialExpert(n, p) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "negativebinomial expert_ll" begin + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for n in nn, p in pp + tmp = LRMoE.NegativeBinomialExpert(n, p) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "negativebinomial exposurize" begin + d = NegativeBinomialExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test mean(LRMoE.exposurize_expert(d; exposure=e)) == mean(d) * e + end +end diff --git a/test/experts/discrete/poisson.jl b/test/experts/discrete/poisson.jl new file mode 100644 index 0000000..944c4c5 --- /dev/null +++ b/test/experts/discrete/poisson.jl @@ -0,0 +1,75 @@ +@testset "poisson pdf, cdf, etc." begin + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:1:50.0) + + for λ in λλ + l = Distributions.Poisson(λ) + r = LRMoE.PoissonExpert(λ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "poisson expert_ll" begin + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:1:50.0) + + for λ in λλ + tmp = LRMoE.PoissonExpert(λ) + + # exact observations + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), LRMoE.logpdf.(tmp, x)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), LRMoE.logpdf.(tmp, x)) + + # right censoring + true_val = + logcdf.(tmp, Inf) .+ log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "poisson exposurize" begin + d = PoissonExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test mean(LRMoE.exposurize_expert(d; exposure=e)) == mean(d) * e + end +end diff --git a/test/experts/discrete/zibinomial.jl b/test/experts/discrete/zibinomial.jl new file mode 100644 index 0000000..aec1060 --- /dev/null +++ b/test/experts/discrete/zibinomial.jl @@ -0,0 +1,116 @@ +@testset "zibinomial pdf, cdf, etc." begin + p0p0 = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for p0 in p0p0, n in nn, p in pp + l = Distributions.Binomial(n, p) + r = LRMoE.ZIBinomialExpert(p0, n, p) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zibinomial expert_ll" begin + p0p0 = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for p0 in p0p0, n in nn, p in pp + tmp = LRMoE.ZIBinomialExpert(p0, n, p) + + # exact observations + true_val = log.(1 - p0) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val) + + true_val = log.(1 - p0) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), true_val) + + # right censoring + true_val = + log.(1 - p0) .+ logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(0.0) .- 1) .- logcdf.(tmp, Inf)) + ), + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.(1 - p0) .+ logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, 1.25 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.75 .* 0.0) .- 1) .- logcdf.(tmp, 1.25 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + log.(1 - p0) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, 1.5 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.25 .* 0.0) .- 1) .- logcdf.(tmp, 1.5 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + log.(1 - p0) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0] .= log1mexp( + log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "zibinomial exposurize" begin + d = ZIBinomialExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test LRMoE.exposurize_expert(d; exposure=e) == d + end +end diff --git a/test/experts/discrete/zigammacount.jl b/test/experts/discrete/zigammacount.jl new file mode 100644 index 0000000..6de6ae2 --- /dev/null +++ b/test/experts/discrete/zigammacount.jl @@ -0,0 +1,117 @@ +# Slightly smaller test sets due to numerical issues (-Inf vs small values) +@testset "zigammacount pdf, cdf, etc." begin + p0p0 = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95] + mm = [1.5, 2, 5, 10] + ss = [1.0, 2.0, 5.0, 10.0] + x = collect(0.0:1:20.0) + + for p0 in p0p0, m in mm, s in ss + l = LRMoE.GammaCount(m, s) + r = LRMoE.ZIGammaCountExpert(p0, m, s) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zigammacount expert_ll" begin + p0p0 = [0.01, 0.05, 0.10, 0.25, 0.90, 0.95] + mm = [1.5, 2, 5, 10] + ss = [1.0, 2.0, 5.0, 10.0] + x = collect(0.0:1:20.0) + + for p0 in p0p0, m in mm, s in ss + tmp = LRMoE.ZIGammaCountExpert(p0, m, s) + + # exact observations + true_val = log.(1 - p0) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val) + + true_val = log.(1 - p0) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), true_val) + + # right censoring + true_val = + log.(1 - p0) .+ logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(0.0) .- 1) .- logcdf.(tmp, Inf)) + ), + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.(1 - p0) .+ logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, 1.25 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.75 .* 0.0) .- 1) .- logcdf.(tmp, 1.25 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + log.(1 - p0) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, 1.5 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.25 .* 0.0) .- 1) .- logcdf.(tmp, 1.5 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + log.(1 - p0) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0] .= log1mexp( + log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +# @testset "zigammacount exposurize" begin +# d = ZINegativeBinomialExpert() +# exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + +# for e in exposure_vec +# @test mean(LRMoE.exposurize_expert(d; exposure=e)) == mean(d) * e +# end +# end diff --git a/test/experts/discrete/zinegativebinomial.jl b/test/experts/discrete/zinegativebinomial.jl new file mode 100644 index 0000000..14e58de --- /dev/null +++ b/test/experts/discrete/zinegativebinomial.jl @@ -0,0 +1,116 @@ +@testset "zinegativebinomial pdf, cdf, etc." begin + p0p0 = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for p0 in p0p0, n in nn, p in pp + l = Distributions.NegativeBinomial(n, p) + r = LRMoE.ZINegativeBinomialExpert(p0, n, p) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zinegativebinomial expert_ll" begin + p0p0 = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + nn = [1, 2, 5, 10, 25, 50, 100] + pp = [0.01, 0.05, 0.10, 0.50, 0.90, 0.95, 0.99] + x = collect(0.0:1:50.0) + + for p0 in p0p0, n in nn, p in pp + tmp = LRMoE.ZINegativeBinomialExpert(p0, n, p) + + # exact observations + true_val = log.(1 - p0) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val) + + true_val = log.(1 - p0) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), true_val) + + # right censoring + true_val = + log.(1 - p0) .+ logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(0.0) .- 1) .- logcdf.(tmp, Inf)) + ), + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.(1 - p0) .+ logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, 1.25 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.75 .* 0.0) .- 1) .- logcdf.(tmp, 1.25 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + log.(1 - p0) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0] .= log( + p0 + + (1 - p0) .* + exp.( + logcdf.(tmp, 1.5 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.25 .* 0.0) .- 1) .- logcdf.(tmp, 1.5 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + log.(1 - p0) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0] .= log1mexp( + log(p0 + (1 - p0) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "zinegativebinomial exposurize" begin + d = ZINegativeBinomialExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test mean(LRMoE.exposurize_expert(d; exposure=e)) == mean(d) * e + end +end diff --git a/test/experts/discrete/zipoisson.jl b/test/experts/discrete/zipoisson.jl new file mode 100644 index 0000000..c914150 --- /dev/null +++ b/test/experts/discrete/zipoisson.jl @@ -0,0 +1,114 @@ +@testset "zipoisson pdf, cdf, etc." begin + pp = [0.05, 0.10, 0.50, 0.90] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:1:50.0) + + for p in pp, λ in λλ + l = Distributions.Poisson(λ) + r = LRMoE.ZIPoissonExpert(p, λ) + @test LRMoE.logpdf.(r, x) ≈ Distributions.logpdf.(l, x) + @test LRMoE.logcdf.(r, x) ≈ Distributions.logcdf.(l, x) + @test LRMoE.pdf.(r, x) ≈ Distributions.pdf.(l, x) + @test LRMoE.cdf.(r, x) ≈ Distributions.cdf.(l, x) + end +end + +@testset "zipoisson expert_ll" begin + pp = [0.05, 0.10, 0.50, 0.90] + λλ = [1.0, 0.5, 0.1, 0.01] + x = collect(0.0:1:50.0) + + for p in pp, λ in λλ + tmp = LRMoE.ZIPoissonExpert(p, λ) + + # exact observations + true_val = log.(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p + (1 - p) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll_exact.(tmp, x), true_val) + + true_val = log.(1 - p) .+ LRMoE.logpdf.(tmp, x) + true_val[x .== 0] .= log(p + (1 - p) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, x, Inf), true_val) + + # right censoring + true_val = + log.(1 - p) .+ logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(x) .- 1) .- logcdf.(tmp, Inf)) + true_val[x .== 0] .= log( + p + + (1 - p) .* + exp.( + logcdf.(tmp, Inf) .+ + log1mexp.(logcdf.(tmp, ceil.(0.0) .- 1) .- logcdf.(tmp, Inf)) + ), + ) + @test isapprox(LRMoE.expert_ll.(tmp, 0.0, x, Inf, Inf), true_val, atol=1e-6) + + true_val = + log.(1 - p) .+ logcdf.(tmp, 1.25 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.75 .* x) .- 1) .- logcdf.(tmp, 1.25 .* x)) + true_val[x .== 0] .= log( + p + + (1 - p) .* + exp.( + logcdf.(tmp, 1.25 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.75 .* 0.0) .- 1) .- logcdf.(tmp, 1.25 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_ll.(tmp, 0.0, 0.75 .* x, 1.25 .* x, Inf), true_val, atol=1e-6 + ) + + # exact observations + @test isapprox(LRMoE.expert_tn.(tmp, 0.0, x, x, Inf), fill(0.0, length(x))) + + # censoring + truncation + true_val = + log.(1 - p) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + true_val[x .== 0] .= log( + p + + (1 - p) .* + exp.( + logcdf.(tmp, 1.5 .* 0.0) .+ + log1mexp.( + logcdf.(tmp, ceil.(0.25 .* 0.0) .- 1) .- logcdf.(tmp, 1.5 .* 0.0) + ) + ), + ) + @test isapprox( + LRMoE.expert_tn.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + + # exact observations + @test isapprox(LRMoE.expert_tn_bar.(tmp, 0.0, x, x, Inf), fill(-Inf, length(x))) + + # censoring + truncation + true_val = + log1mexp.( + log.(1 - p) .+ logcdf.(tmp, 1.5 .* x) .+ + log1mexp.(logcdf.(tmp, ceil.(0.25 .* x) .- 1) .- logcdf.(tmp, 1.5 .* x)) + ) + true_val[x .== 0] .= log1mexp( + log(p + (1 - p) .* exp.(LRMoE.logpdf.(tmp, 0.0))) + ) + @test isapprox( + LRMoE.expert_tn_bar.(tmp, 0.25 .* x, 0.75 .* x, 1.25 .* x, 1.5 .* x), + true_val, + atol=1e-6, + ) + end +end + +@testset "zipoisson exposurize" begin + d = ZIPoissonExpert() + exposure_vec = [0.5, 1.0, 2.0, 5.0, 10.0] + + for e in exposure_vec + @test mean(LRMoE.exposurize_expert(d; exposure=e)) == mean(d) * e + end +end diff --git a/test/experts/ll/exposurize.jl b/test/experts/ll/exposurize.jl new file mode 100644 index 0000000..aaec471 --- /dev/null +++ b/test/experts/ll/exposurize.jl @@ -0,0 +1,413 @@ +@testset "exposurize experts" begin + param0_list = [0.1, 0.5, 0.9] + param1_list = [10.0, 15.0, 20.0] + param2_list = [0.1, 0.5, 0.9] + param3_list = [5.0, 7.0, 10.0] + + sample_size = 21 + + expo_list = [0.5, 1.0, 2.0, 5.0, 10.0] + + for param0 in param0_list, + param1 in param1_list, + param2 in param2_list, + param3 in param3_list, + expo in expo_list + + tmp = LRMoE.BinomialExpert(20, param0) + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + tmp = LRMoE.ZIBinomialExpert(0.20, 20, param0) + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + + tmp = LRMoE.GammaCountExpert(param1, param2) + if expo != 1.0 + @test tmp != LRMoE.exposurize_expert(tmp; exposure=expo) + else + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + end + @test isapprox( + [LRMoE.params(tmp)...] .* [1, 1 / expo], + [LRMoE.params(LRMoE.exposurize_expert(tmp; exposure=expo))...], + atol=1e-6, + ) + + tmp = LRMoE.ZIGammaCountExpert(param0, param1, param2) + if expo != 1.0 + @test tmp != LRMoE.exposurize_expert(tmp; exposure=expo) + else + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + end + @test isapprox( + [LRMoE.params(tmp)...] .* [1, 1, 1 / expo], + [LRMoE.params(LRMoE.exposurize_expert(tmp; exposure=expo))...], + atol=1e-6, + ) + + tmp = LRMoE.NegativeBinomialExpert(param1, param0) + if expo != 1.0 + @test tmp != LRMoE.exposurize_expert(tmp; exposure=expo) + else + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + end + @test [LRMoE.params(tmp)...] .* [expo, 1] == + [LRMoE.params(LRMoE.exposurize_expert(tmp; exposure=expo))...] + + tmp = LRMoE.ZINegativeBinomialExpert(param0, param1, param0) + if expo != 1.0 + @test tmp != LRMoE.exposurize_expert(tmp; exposure=expo) + else + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + end + @test [LRMoE.params(tmp)...] .* [1, expo, 1] == + [LRMoE.params(LRMoE.exposurize_expert(tmp; exposure=expo))...] + + tmp = LRMoE.PoissonExpert(param1) + if expo != 1.0 + @test tmp != LRMoE.exposurize_expert(tmp; exposure=expo) + else + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + end + @test [LRMoE.params(tmp)...] .* [expo] == + [LRMoE.params(LRMoE.exposurize_expert(tmp; exposure=expo))...] + if expo != 1.0 + @test tmp != LRMoE.exposurize_expert(tmp; exposure=expo) + else + @test tmp == LRMoE.exposurize_expert(tmp; exposure=expo) + end + @test [LRMoE.params(tmp)...] .* [expo] == + [LRMoE.params(LRMoE.exposurize_expert(tmp; exposure=expo))...] + + model = [ + LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3) + BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0) + ] + + expos = fill(expo, sample_size) + + @test size(LRMoE.exposurize_model(model; exposure=expos)) == (2, 4, sample_size) + + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 1, :] == + fill(mean(model[1, 1]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 2, :] == + fill(mean(model[1, 2]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 3, :] == + fill(mean(model[1, 3]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 4, :] == + fill(mean(model[1, 4]), sample_size) + + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[2, 1, :] == + fill(mean(model[2, 1]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[2, 2, :] == + fill(mean(model[2, 2]) .* expo, sample_size) + + model_expod = LRMoE.exposurize_model(model; exposure=expos) + Y_tmp = hcat( + collect(0.0:0.5:10), collect(0:1:20) + ) + + @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[1, :], model_expod[:, :, 1])) == + (2, 4) + @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[5, :], model_expod[:, :, 5])) == + (2, 4) + @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[12, :], model_expod[:, :, 12])) == + (2, 4) + @test size(LRMoE.expert_ll_list_exact(Y_tmp, model_expod)) == (2, 4, sample_size) + + for idx in [1, 5, 10, 20, 21] + mat_tmp = LRMoE.expert_ll_ind_mat_exact(Y_tmp[idx, :], model_expod[:, :, idx]) + @test mat_tmp[1, 1] == + LRMoE.expert_ll_exact(model_expod[1, 1, idx], Y_tmp[idx, 1]) + @test mat_tmp[1, 2] == + LRMoE.expert_ll_exact(model_expod[1, 2, idx], Y_tmp[idx, 1]) + @test mat_tmp[1, 3] == + LRMoE.expert_ll_exact(model_expod[1, 3, idx], Y_tmp[idx, 1]) + @test mat_tmp[1, 4] == + LRMoE.expert_ll_exact(model_expod[1, 4, idx], Y_tmp[idx, 1]) + @test mat_tmp[2, 1] == + LRMoE.expert_ll_exact(model_expod[2, 1, idx], Y_tmp[idx, 2]) + @test mat_tmp[2, 2] == + LRMoE.expert_ll_exact(model_expod[2, 2, idx], Y_tmp[idx, 2]) + @test mat_tmp[2, 3] == + LRMoE.expert_ll_exact(model_expod[2, 3, idx], Y_tmp[idx, 2]) + @test mat_tmp[2, 3] == + LRMoE.expert_ll_exact(model_expod[2, 3, idx], Y_tmp[idx, 2]) + end + + for idx in [1, 5, 10, 20, 21] + cube_tmp = LRMoE.expert_ll_list_exact(Y_tmp, model_expod) + @test cube_tmp[1, 1, idx] == + LRMoE.expert_ll_exact(model_expod[1, 1, idx], Y_tmp[idx, 1]) + @test cube_tmp[1, 2, idx] == + LRMoE.expert_ll_exact(model_expod[1, 2, idx], Y_tmp[idx, 1]) + @test cube_tmp[1, 3, idx] == + LRMoE.expert_ll_exact(model_expod[1, 3, idx], Y_tmp[idx, 1]) + @test cube_tmp[1, 4, idx] == + LRMoE.expert_ll_exact(model_expod[1, 4, idx], Y_tmp[idx, 1]) + @test cube_tmp[2, 1, idx] == + LRMoE.expert_ll_exact(model_expod[2, 1, idx], Y_tmp[idx, 2]) + @test cube_tmp[2, 2, idx] == + LRMoE.expert_ll_exact(model_expod[2, 2, idx], Y_tmp[idx, 2]) + @test cube_tmp[2, 3, idx] == + LRMoE.expert_ll_exact(model_expod[2, 3, idx], Y_tmp[idx, 2]) + @test cube_tmp[2, 4, idx] == + LRMoE.expert_ll_exact(model_expod[2, 4, idx], Y_tmp[idx, 2]) + end + end +end + +@testset "exposurize model" begin + param0_list = [0.1, 0.5, 0.9] + param1_list = [10.0, 15.0, 20.0] + param2_list = [0.1, 0.5, 0.9] + param3_list = [5.0, 7.0, 10.0] + + sample_size = 21 + + expo_list = [0.5, 1.0, 2.0, 5.0, 10.0] + + for param0 in param0_list, + param1 in param1_list, + param2 in param2_list, + param3 in param3_list, + expo in expo_list + + model = [ + LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3) + BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0) + ] + + expos = fill(expo, sample_size) + + @test size(LRMoE.exposurize_model(model; exposure=expos)) == (2, 4, sample_size) + + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 1, :] == + fill(mean(model[1, 1]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 2, :] == + fill(mean(model[1, 2]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 3, :] == + fill(mean(model[1, 3]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[1, 4, :] == + fill(mean(model[1, 4]), sample_size) + + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[2, 1, :] == + fill(mean(model[2, 1]), sample_size) + @test mean.(LRMoE.exposurize_model(model; exposure=expos))[2, 2, :] == + fill(mean(model[2, 2]) .* expo, sample_size) + + model_expod = LRMoE.exposurize_model(model; exposure=expos) + Y_tmp = hcat( + collect(0.0:0.5:10), collect(0:1:20) + ) + + @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[1, :], model_expod[:, :, 1])) == + (2, 4) + @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[5, :], model_expod[:, :, 5])) == + (2, 4) + @test size(LRMoE.expert_ll_ind_mat_exact(Y_tmp[12, :], model_expod[:, :, 12])) == + (2, 4) + @test size(LRMoE.expert_ll_list_exact(Y_tmp, model_expod)) == (2, 4, sample_size) + + for idx in [1, 5, 10, 20, 21] + mat_tmp = LRMoE.expert_ll_ind_mat_exact(Y_tmp[idx, :], model_expod[:, :, idx]) + @test mat_tmp[1, 1] == + LRMoE.expert_ll_exact(model_expod[1, 1, idx], Y_tmp[idx, 1]) + @test mat_tmp[1, 2] == + LRMoE.expert_ll_exact(model_expod[1, 2, idx], Y_tmp[idx, 1]) + @test mat_tmp[1, 3] == + LRMoE.expert_ll_exact(model_expod[1, 3, idx], Y_tmp[idx, 1]) + @test mat_tmp[1, 4] == + LRMoE.expert_ll_exact(model_expod[1, 4, idx], Y_tmp[idx, 1]) + @test mat_tmp[2, 1] == + LRMoE.expert_ll_exact(model_expod[2, 1, idx], Y_tmp[idx, 2]) + @test mat_tmp[2, 2] == + LRMoE.expert_ll_exact(model_expod[2, 2, idx], Y_tmp[idx, 2]) + @test mat_tmp[2, 3] == + LRMoE.expert_ll_exact(model_expod[2, 3, idx], Y_tmp[idx, 2]) + @test mat_tmp[2, 3] == + LRMoE.expert_ll_exact(model_expod[2, 3, idx], Y_tmp[idx, 2]) + end + + for idx in [1, 5, 10, 20, 21] + cube_tmp = LRMoE.expert_ll_list_exact(Y_tmp, model_expod) + @test cube_tmp[1, 1, idx] == + LRMoE.expert_ll_exact(model_expod[1, 1, idx], Y_tmp[idx, 1]) + @test cube_tmp[1, 2, idx] == + LRMoE.expert_ll_exact(model_expod[1, 2, idx], Y_tmp[idx, 1]) + @test cube_tmp[1, 3, idx] == + LRMoE.expert_ll_exact(model_expod[1, 3, idx], Y_tmp[idx, 1]) + @test cube_tmp[1, 4, idx] == + LRMoE.expert_ll_exact(model_expod[1, 4, idx], Y_tmp[idx, 1]) + @test cube_tmp[2, 1, idx] == + LRMoE.expert_ll_exact(model_expod[2, 1, idx], Y_tmp[idx, 2]) + @test cube_tmp[2, 2, idx] == + LRMoE.expert_ll_exact(model_expod[2, 2, idx], Y_tmp[idx, 2]) + @test cube_tmp[2, 3, idx] == + LRMoE.expert_ll_exact(model_expod[2, 3, idx], Y_tmp[idx, 2]) + @test cube_tmp[2, 4, idx] == + LRMoE.expert_ll_exact(model_expod[2, 4, idx], Y_tmp[idx, 2]) + end + + Y_tmp = hcat( + 0.2 .* Y_tmp[:, 1], + 0.8 .* Y_tmp[:, 1], + 1.2 .* Y_tmp[:, 1], + 5 .* Y_tmp[:, 1], + fill(0, length(Y_tmp[:, 2])), + Y_tmp[:, 2], + Y_tmp[:, 2], + fill(Inf, length(Y_tmp[:, 2])), + ) + + @test size(LRMoE.expert_ll_ind_mat(Y_tmp[1, :], model_expod[:, :, 1])) == (2, 4) + @test size(LRMoE.expert_ll_ind_mat(Y_tmp[5, :], model_expod[:, :, 5])) == (2, 4) + @test size(LRMoE.expert_ll_ind_mat(Y_tmp[12, :], model_expod[:, :, 12])) == (2, 4) + @test size(LRMoE.expert_ll_list(Y_tmp, model_expod)) == (2, 4, sample_size) + + for idx in [1, 5, 10, 20, 21] + mat_tmp = LRMoE.expert_ll_ind_mat(Y_tmp[idx, :], model_expod[:, :, idx]) + @test mat_tmp[1, 1] == LRMoE.expert_ll( + model_expod[1, 1, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test mat_tmp[1, 2] == LRMoE.expert_ll( + model_expod[1, 2, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test mat_tmp[1, 3] == LRMoE.expert_ll( + model_expod[1, 3, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test mat_tmp[1, 4] == LRMoE.expert_ll( + model_expod[1, 4, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test mat_tmp[2, 1] == LRMoE.expert_ll( + model_expod[2, 1, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + @test mat_tmp[2, 2] == LRMoE.expert_ll( + model_expod[2, 2, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + @test mat_tmp[2, 3] == LRMoE.expert_ll( + model_expod[2, 3, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + @test mat_tmp[2, 3] == LRMoE.expert_ll( + model_expod[2, 3, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + end + + cube_tmp = LRMoE.expert_ll_list(Y_tmp, model_expod) + + for idx in [1, 5, 10, 20, 21] + @test cube_tmp[1, 1, idx] == LRMoE.expert_ll( + model_expod[1, 1, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test cube_tmp[1, 2, idx] == LRMoE.expert_ll( + model_expod[1, 2, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test cube_tmp[1, 3, idx] == LRMoE.expert_ll( + model_expod[1, 3, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test cube_tmp[1, 4, idx] == LRMoE.expert_ll( + model_expod[1, 4, idx], + Y_tmp[idx, 1], + Y_tmp[idx, 2], + Y_tmp[idx, 3], + Y_tmp[idx, 4], + ) + @test cube_tmp[2, 1, idx] == LRMoE.expert_ll( + model_expod[2, 1, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + @test cube_tmp[2, 2, idx] == LRMoE.expert_ll( + model_expod[2, 2, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + @test cube_tmp[2, 3, idx] == LRMoE.expert_ll( + model_expod[2, 3, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + @test cube_tmp[2, 4, idx] == LRMoE.expert_ll( + model_expod[2, 4, idx], + Y_tmp[idx, 5], + Y_tmp[idx, 6], + Y_tmp[idx, 7], + Y_tmp[idx, 8], + ) + end + + # ll related + Random.seed!(1234) + + gate_tmp = rand(Distributions.Uniform(-1, 1), sample_size, 4) + + dim_agg = LRMoE.loglik_aggre_dim(cube_tmp) + @test size(dim_agg) == (sample_size, 4) + @test size(LRMoE.loglik_aggre_gate_dim(gate_tmp, dim_agg)) == (sample_size, 4) + + # Simulation related + # X_tmp = rand(Distributions.Uniform(-1, 1), 15, 7) + # α_tmp = rand(Distributions.Uniform(-1, 1), 4, 7) + + # model = [LogNormalExpert(param1, param2) ZIGammaExpert(param0, param1, param2) WeibullExpert(param1, param2) BurrExpert(param1, param2, param3); + # BinomialExpert(20, param0) PoissonExpert(param1) ZIGammaCountExpert(param0, param1, param2) ZINegativeBinomialExpert(param0, param1, param0)] + + # expo = rand(Distributions.Uniform(0.5, 20), 15) + + # model_expod = LRMoE.exposurize_model(model, exposure = expo) + + # sim_tmp = LRMoE.sim_dataset(α_tmp, X_tmp, model, exposure = expo) + # mean(sim_tmp, dims = 1) + + # sim_tmp = LRMoE.sim_dataset(α_tmp, X_tmp, model) + # mean(sim_tmp, dims = 1) + end +end \ No newline at end of file diff --git a/test/gating.jl b/test/gating.jl index 4fce092..b21e9c5 100644 --- a/test/gating.jl +++ b/test/gating.jl @@ -1,18 +1,15 @@ -using Test -using Distributions - @testset "logit gating" begin + Random.seed!(1234) d = Distributions.Uniform(-2.0, 2.0) - for i in 1:20 + for i in 1:5 α = rand(d, 5, 20) - x = rand(d, 100, 20) + x = rand(d, 100, 20) gating = LogitGating(α, x) ax = x * α' - result = ax .- log.(sum(exp.(ax), dims = 2)) + result = ax .- log.(sum(exp.(ax); dims=2)) @test gating ≈ result - @test sum(exp.(gating), dims = 2) ≈ fill(1, size(gating)[1]) + @test sum(exp.(gating); dims=2) ≈ fill(1, size(gating)[1]) end - end \ No newline at end of file diff --git a/test/loglik.jl b/test/loglik.jl index a75a37b..95df255 100644 --- a/test/loglik.jl +++ b/test/loglik.jl @@ -3,216 +3,92 @@ using Distributions using StatsFuns @testset "loglik list (individual)" begin + μμ = [1.0, 2.0] + σσ = [1.0, 0.5] - d = Distributions.Gamma(1.0, 2.0) - - μμ = rand(d, 5) - σσ = rand(d, 5) - - for (μ, σ) in zip(μμ, σσ) + for μ in μμ, σ in σσ l = Distributions.LogNormal(μ, σ) - y = rand(l, 25) + y = collect(0:1:24) + sample_size = length(y) + expos = fill(1.0, sample_size) # LogNormal r = LRMoE.LogNormalExpert(μ, σ) - - # Y = hcat(y, y, y, y) - # model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ)] - - Y = hcat(fill(0, length(y)), y, y, fill(Inf, length(y)), fill(0, length(y)), 0.80.*y, 1.25.*y, fill(Inf, length(y))) - model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ); - LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.80, 1.5*μ, 2.0*σ)] - + + Y = hcat( + fill(0, sample_size), + y, + y, + fill(Inf, sample_size), + fill(0, sample_size), + 0.80 .* y, + 1.25 .* y, + fill(Inf, sample_size), + ) + + model_raw = [ + LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5 * μ, 0.6 * σ) LRMoE.LogNormalExpert(1.5 * μ, 2.0 * σ) + LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5 * μ, 0.6 * σ) LRMoE.ZILogNormalExpert(0.80, 1.5 * μ, 2.0 * σ) + ] + + model = LRMoE.exposurize_model(model_raw; exposure=expos) + # Test size - @test size(expert_ll_pos_list(Y, model))[1] == 2 - @test size(expert_ll_pos_list(Y, model)[1])[1] == 25 - @test size(expert_ll_pos_list(Y, model)[1])[2] == 3 - - @test size(expert_tn_pos_list(Y, model))[1] == 2 - @test size(expert_tn_pos_list(Y, model)[1])[1] == 25 - @test size(expert_tn_pos_list(Y, model)[1])[2] == 3 - - @test size(expert_tn_bar_pos_list(Y, model))[1] == 2 - @test size(expert_tn_bar_pos_list(Y, model)[1])[1] == 25 - @test size(expert_tn_bar_pos_list(Y, model)[1])[2] == 3 - - @test size(expert_ll_list(Y, model))[1] == 2 - @test size(expert_ll_list(Y, model)[1])[1] == 25 - @test size(expert_ll_list(Y, model)[1])[2] == 3 - - @test size(expert_tn_list(Y, model))[1] == 2 - @test size(expert_tn_list(Y, model)[1])[1] == 25 - @test size(expert_tn_list(Y, model)[1])[2] == 3 - - @test size(expert_tn_bar_list(Y, model))[1] == 2 - @test size(expert_tn_bar_list(Y, model)[1])[1] == 25 - @test size(expert_tn_bar_list(Y, model)[1])[2] == 3 - - # 1*1 model - y = rand(l, 1) - Y = [0 y y Inf] - model = reshape([LRMoE.LogNormalExpert(μ, σ)], (1, 1)) - @test expert_ll_pos_list(Y, model)[1] ≈ [Distributions.logpdf.(l, Y[2])] - @test expert_tn_pos_list(Y, model)[1] ≈ [0.0] - @test expert_tn_bar_pos_list(Y, model)[1] ≈ [-Inf] - @test expert_ll_list(Y, model)[1] ≈ [Distributions.logpdf.(l, Y[2])] - @test expert_tn_list(Y, model)[1] ≈ [0.0] - @test expert_tn_bar_list(Y, model)[1] ≈ [-Inf] - - y = rand(l, 1) - Y = [0 0.80*y 1.25*y Inf] - model = reshape([LRMoE.LogNormalExpert(μ, σ)], (1, 1)) - @test expert_ll_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_tn_pos_list(Y, model)[1] ≈ [0.0] - @test expert_tn_bar_pos_list(Y, model)[1] ≈ [-Inf] - @test expert_ll_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_tn_list(Y, model)[1] ≈ [0.0] - @test expert_tn_bar_list(Y, model)[1] ≈ [-Inf] - - y = rand(l, 1) - Y = [0.50*y 0.80*y 1.25*y 2.00*y] - model = reshape([LRMoE.LogNormalExpert(μ, σ)], (1, 1)) - @test expert_ll_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp.(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_tn_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))] - @test expert_tn_bar_pos_list(Y, model)[1] ≈ log1mexp.([Distributions.logcdf.(l, Y[4]) + log1mexp(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))]) - @test expert_ll_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp.(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_tn_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))] - @test expert_tn_bar_list(Y, model)[1] ≈ log1mexp.([Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))]) - - # 2*1 model - p = rand(Distributions.Uniform(0.0, 1.0), 1)[1] - y = rand(l, 2) .+ 0.5 - Y = [0.0 0.80*y[1] 1.25*y[1] Inf 0.0 0.80*y[2] 1.25*y[2] Inf] - model = reshape([LRMoE.LogNormalExpert(μ, σ); - LRMoE.ZILogNormalExpert(p, μ, σ)], (2, 1)) - @test expert_ll_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp.(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_ll_pos_list(Y, model)[2] ≈ [Distributions.logcdf.(l, Y[7]) + log1mexp.(Distributions.logcdf.(l, Y[6]) - Distributions.logcdf.(l, Y[7]))] - - @test expert_tn_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))] - @test expert_tn_pos_list(Y, model)[2] ≈ [Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))] - - @test expert_tn_bar_pos_list(Y, model)[1] ≈ log1mexp.([Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))]) - @test expert_tn_bar_pos_list(Y, model)[2] ≈ log1mexp.([Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))]) - - @test expert_ll_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp.(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_ll_list(Y, model)[2] ≈ [log(1-p) .+ Distributions.logcdf.(l, Y[7]) + log1mexp.(Distributions.logcdf.(l, Y[6]) - Distributions.logcdf.(l, Y[7]))] - - @test expert_tn_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))] - @test expert_tn_list(Y, model)[2] ≈ [Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))] - - @test expert_tn_bar_list(Y, model)[1] ≈ log1mexp.([Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))]) - @test expert_tn_bar_list(Y, model)[2] ≈ log1mexp.([Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))]) - - - p = rand(Distributions.Uniform(0.0, 1.0), 1)[1] - y = rand(l, 2) .+ 0.5 - Y = [0.0 0.80*y[1] 1.25*y[1] Inf 0.50*y[2] 0.80*y[2] 1.25*y[2] 2.50*y[2]] - model = reshape([LRMoE.LogNormalExpert(μ, σ); - LRMoE.ZILogNormalExpert(p, μ, σ)], (2, 1)) - @test expert_ll_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp.(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_ll_pos_list(Y, model)[2] ≈ [Distributions.logcdf.(l, Y[7]) + log1mexp.(Distributions.logcdf.(l, Y[6]) - Distributions.logcdf.(l, Y[7]))] - - @test expert_tn_pos_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))] - @test expert_tn_pos_list(Y, model)[2] ≈ [Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))] - - @test expert_tn_bar_pos_list(Y, model)[1] ≈ log1mexp.([Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))]) - @test expert_tn_bar_pos_list(Y, model)[2] ≈ log1mexp.([Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))]) - - @test expert_ll_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[3]) + log1mexp.(Distributions.logcdf.(l, Y[2]) - Distributions.logcdf.(l, Y[3]))] - @test expert_ll_list(Y, model)[2] ≈ [log(1-p) .+ Distributions.logcdf.(l, Y[7]) + log1mexp.(Distributions.logcdf.(l, Y[6]) - Distributions.logcdf.(l, Y[7]))] - - @test expert_tn_list(Y, model)[1] ≈ [Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))] - @test expert_tn_list(Y, model)[2] ≈ [log(1-p) .+ Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8]))] - - @test expert_tn_bar_list(Y, model)[1] ≈ log1mexp.([Distributions.logcdf.(l, Y[4]) + log1mexp.(Distributions.logcdf.(l, Y[1]) - Distributions.logcdf.(l, Y[4]))]) - @test expert_tn_bar_list(Y, model)[2] ≈ [log.(p .+ (1-p) .* exp.(log1mexp.(Distributions.logcdf.(l, Y[8]) + log1mexp.(Distributions.logcdf.(l, Y[5]) - Distributions.logcdf.(l, Y[8])))))] - - # 2*3 model - p1 = rand(Distributions.Uniform(0.0, 1.0), 1)[1] - p2 = rand(Distributions.Uniform(0.0, 1.0), 1)[1] - y = rand(l, 100, 2) .+ 0.5 - Y = [fill(0.0, length(y[:,1])) 0.80.*y[:,1] 1.25.*y[:,1] fill(Inf, length(y[:,1])) 0.50.*y[:,2] 0.80.*y[:,2] 1.25.*y[:,2] 2.50.*y[:,2]] - model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.ZILogNormalExpert(p1, 1.2*μ, 0.8*σ); - LRMoE.ZILogNormalExpert(p2, 2*μ, 0.5*σ) LRMoE.LogNormalExpert(0.9*μ, 0.75*σ)] - l1 = Distributions.LogNormal(μ, σ) - l2 = Distributions.LogNormal(1.2*μ, 0.8*σ) - l3 = Distributions.LogNormal(2*μ, 0.5*σ) - l4 = Distributions.LogNormal(0.9*μ, 0.75*σ) - - @test expert_ll_pos_list(Y, model)[1] ≈ hcat(Distributions.logcdf.(l1, Y[:,3]) + log1mexp.(Distributions.logcdf.(l1, Y[:,2]) - Distributions.logcdf.(l1, Y[:,3])), - Distributions.logcdf.(l2, Y[:,3]) + log1mexp.(Distributions.logcdf.(l2, Y[:,2]) - Distributions.logcdf.(l2, Y[:,3]))) - @test expert_ll_pos_list(Y, model)[2] ≈ hcat(Distributions.logcdf.(l3, Y[:,7]) + log1mexp.(Distributions.logcdf.(l3, Y[:,6]) - Distributions.logcdf.(l3, Y[:,7])), - Distributions.logcdf.(l4, Y[:,7]) + log1mexp.(Distributions.logcdf.(l4, Y[:,6]) - Distributions.logcdf.(l4, Y[:,7]))) - - @test expert_tn_pos_list(Y, model)[1] ≈ hcat(fill(0.0, length(Y[:,2])), - Distributions.logcdf.(l2, Y[:,4]) + log1mexp.(Distributions.logcdf.(l2, Y[:,1]) - Distributions.logcdf.(l2, Y[:,4]))) - @test expert_tn_pos_list(Y, model)[2] ≈ hcat(Distributions.logcdf.(l3, Y[:,8]) + log1mexp.(Distributions.logcdf.(l3, Y[:,5]) - Distributions.logcdf.(l3, Y[:,8])), - Distributions.logcdf.(l4, Y[:,8]) + log1mexp.(Distributions.logcdf.(l4, Y[:,5]) - Distributions.logcdf.(l4, Y[:,8]))) - - # @test expert_tn_bar_pos_list(Y, model)[1] ≈ hcat(fill(-Inf, length(Y[:,2])), - # log1mexp.(Distributions.logcdf.(l2, Y[:,4]) + log1mexp.(Distributions.logcdf.(l2, Y[:,1]) - Distributions.logcdf.(l2, Y[:,4])))) - # @test expert_tn_bar_pos_list(Y, model)[2] ≈ hcat(log1mexp.(Distributions.logcdf.(l3, Y[:,8]) + log1mexp.(Distributions.logcdf.(l3, Y[:,5]) - Distributions.logcdf.(l3, Y[:,8]))), - # log1mexp.(Distributions.logcdf.(l4, Y[:,8]) + log1mexp.(Distributions.logcdf.(l4, Y[:,5]) - Distributions.logcdf.(l4, Y[:,8])))) - - # @test expert_ll_list(Y, model)[1] ≈ hcat(Distributions.logcdf.(l1, Y[:,3]) + log1mexp.(Distributions.logcdf.(l1, Y[:,2]) - Distributions.logcdf.(l1, Y[:,3])), - # log(1-p1) .+ Distributions.logcdf.(l2, Y[:,3]) + log1mexp.(Distributions.logcdf.(l2, Y[:,2]) - Distributions.logcdf.(l2, Y[:,3]))) - # @test expert_ll_list(Y, model)[2] ≈ hcat(log(1-p2) .+ Distributions.logcdf.(l3, Y[:,7]) + log1mexp.(Distributions.logcdf.(l3, Y[:,6]) - Distributions.logcdf.(l3, Y[:,7])), - # Distributions.logcdf.(l4, Y[:,7]) + log1mexp.(Distributions.logcdf.(l4, Y[:,6]) - Distributions.logcdf.(l4, Y[:,7]))) - - @test expert_tn_list(Y, model)[1] ≈ hcat(fill(0.0, length(Y[:,2])), - fill(0.0, length(Y[:,2]))) - - - # @test expert_tn_list(Y, model)[2] ≈ hcat(log(1-p2) .+ Distributions.logcdf.(l3, Y[:,8]) .+ log1mexp.(Distributions.logcdf.(l3, Y[:,5]) - Distributions.logcdf.(l3, Y[:,8])), - # Distributions.logcdf.(l4, Y[:,8]) + log1mexp.(Distributions.logcdf.(l4, Y[:,5]) - Distributions.logcdf.(l4, Y[:,8]))) - - # @test isapprox(expert_tn_list(Y, model)[2], hcat(log(1-p2) .+ Distributions.logcdf.(l3, Y[:,8]) .+ log1mexp.(Distributions.logcdf.(l3, Y[:,5]) - Distributions.logcdf.(l3, Y[:,8])), - # Distributions.logcdf.(l4, Y[:,8]) + log1mexp.(Distributions.logcdf.(l4, Y[:,5]) - Distributions.logcdf.(l4, Y[:,8]))), atol = 1e-03) - - @test expert_tn_bar_list(Y, model)[1] ≈ hcat(fill(-Inf, length(Y[:,2])), - fill(-Inf, length(Y[:,2]))) - - # @test expert_tn_bar_list(Y, model)[2] ≈ hcat(log.(p2.+ (1-p2) .* (1 .- exp.(Distributions.logcdf.(l3, Y[:,8]) + log1mexp.(Distributions.logcdf.(l3, Y[:,5]) - Distributions.logcdf.(l3, Y[:,8]))))), - # log1mexp.(Distributions.logcdf.(l4, Y[:,8]) + log1mexp.(Distributions.logcdf.(l4, Y[:,5]) - Distributions.logcdf.(l4, Y[:,8])))) - - # @test isapprox(expert_tn_bar_list(Y, model)[2], hcat(log.(p2.+ (1-p2) .* (1 .- exp.(Distributions.logcdf.(l3, Y[:,8]) + log1mexp.(Distributions.logcdf.(l3, Y[:,5]) - Distributions.logcdf.(l3, Y[:,8]))))), - # log1mexp.(Distributions.logcdf.(l4, Y[:,8]) + log1mexp.(Distributions.logcdf.(l4, Y[:,5]) - Distributions.logcdf.(l4, Y[:,8])))), atol = 1e-03) - - end + @test size(expert_ll_list(Y, model)[:, :, 1]) == (2, 3) + @test size(expert_ll_list(Y, model))[3] == sample_size + @test size(expert_tn_list(Y, model)[:, :, 1]) == (2, 3) + @test size(expert_tn_list(Y, model))[3] == sample_size + @test size(expert_tn_bar_list(Y, model)[:, :, 1]) == (2, 3) + @test size(expert_tn_bar_list(Y, model))[3] == sample_size + end end @testset "loglik list (aggredated) " begin - d = Distributions.Gamma(1.0, 2.0) - - μμ = rand(d, 5) - σσ = rand(d, 5) + μμ = [1.0, 2.0] + σσ = [1.0, 0.5] - for (μ, σ) in zip(μμ, σσ) + for μ in μμ, σ in σσ l = Distributions.LogNormal(μ, σ) - y = rand(l, 25) + y = collect(0:1:24) + sample_size = length(y) + expos = fill(1.0, sample_size) # LogNormal r = LRMoE.LogNormalExpert(μ, σ) - - X = rand(Uniform(-1, 1), 25, 10) + + Random.seed!(1234) + + X = rand(Uniform(-1, 1), sample_size, 10) α = rand(Uniform(-1, 1), 3, 10) - α[3,:] .= 0.0 + α[3, :] .= 0.0 + + Y = hcat( + fill(0, sample_size), + y, + y, + fill(Inf, sample_size), + fill(0, sample_size), + 0.80 .* y, + 1.25 .* y, + fill(Inf, sample_size), + ) + model_raw = [ + LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5 * μ, 0.6 * σ) LRMoE.LogNormalExpert(1.5 * μ, 2.0 * σ) + LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5 * μ, 0.6 * σ) LRMoE.ZILogNormalExpert(0.80, 1.5 * μ, 2.0 * σ) + ] + model = LRMoE.exposurize_model(model_raw; exposure=expos) - Y = hcat(fill(0, length(y)), y, y, fill(Inf, length(y)), fill(0, length(y)), 0.80.*y, 1.25.*y, fill(Inf, length(y))) - model = [LRMoE.LogNormalExpert(μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.LogNormalExpert(1.5*μ, 2.0*σ); - LRMoE.ZILogNormalExpert(0.4, μ, σ) LRMoE.LogNormalExpert(0.5*μ, 0.6*σ) LRMoE.ZILogNormalExpert(0.80, 1.5*μ, 2.0*σ)] - # Test size - @test size(LogitGating(α, X))[1] == 25 + @test size(LogitGating(α, X))[1] == sample_size @test size(LogitGating(α, X))[2] == 3 gate = LogitGating(α, X) ll_ls = loglik_np(Y, gate, model) # Test number @test isa(ll_ls.ll, Number) - @test exp.(ll_ls.gate_expert_tn) + exp.(ll_ls.gate_expert_tn_bar) ≈ fill(1.0, length(ll_ls.gate_expert_tn)) - + @test exp.(ll_ls.gate_expert_tn) + exp.(ll_ls.gate_expert_tn_bar) ≈ + fill(1.0, length(ll_ls.gate_expert_tn)) end end \ No newline at end of file diff --git a/test/pdfcdf.jl b/test/pdfcdf.jl deleted file mode 100644 index b9bdc7e..0000000 --- a/test/pdfcdf.jl +++ /dev/null @@ -1,21 +0,0 @@ -using Test -using Distributions - -@testset "override pdf, cdf, etc." begin - - d = Distributions.Gamma(1.0, 2.0) - - μμ = rand(d, 50) - σσ = rand(d, 50) - - for (μ, σ) in zip(μμ, σσ) - l = Distributions.LogNormal(μ, σ) - r = LRMoE.LogNormalExpert(μ, σ) - x = rand(l, 100) - @test LRMoE.logpdf(r, x) ≈ Distributions.logpdf.(l, x) - @test LRMoE.logcdf(r, x) ≈ Distributions.logcdf.(l, x) - @test LRMoE.pdf(r, x) ≈ Distributions.pdf.(l, x) - @test LRMoE.cdf(r, x) ≈ Distributions.cdf.(l, x) - end - -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index a10a951..5836c9d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,33 +1,42 @@ -# using Distributions -# using PDMats # test dependencies -# using Distributed -# using StatsBase -# using LinearAlgebra -# using HypothesisTests - -# import JSON -# import ForwardDiff - using LRMoE using Distributions -# using Random +using QuadGK +using Random +using StatsFuns using Test const tests = [ - # "dummytest", - # "pdfcdf", - # "expert_loglik", - # "gating", - # "loglik", - # "fit_functions" + # gating + "gating", + # exposurize + "experts/ll/exposurize", + # continuous experts + "experts/continuous/burr", + "experts/continuous/gamma", + "experts/continuous/inversegaussian", + "experts/continuous/lognormal", + "experts/continuous/weibull", + # zi continuous experts + "experts/continuous/ziburr", + "experts/continuous/zigamma", + "experts/continuous/ziinversegaussian", + "experts/continuous/zilognormal", + "experts/continuous/ziweibull", + # discrete experts + "experts/discrete/binomial", + "experts/discrete/gammacount", + "experts/discrete/negativebinomial", + "experts/discrete/poisson", + # zi discrete experts + "experts/discrete/zibinomial", + "experts/discrete/zigammacount", + "experts/discrete/zinegativebinomial", + "experts/discrete/zipoisson", + # loglik + "loglik", ] -printstyled("Running tests:\n", color=:blue) - -# Random.seed!(345679) - -# to reduce redundancy, we might break this file down into seperate `$t * "_utils.jl"` files -# include("testutils.jl") +printstyled("Running tests:\n"; color=:blue) for t in tests @testset "Test $t" begin