From ba63fda96208b115781fd784e3048dfe70968afe Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 10 Jul 2023 16:11:14 +0100 Subject: [PATCH 01/10] Fix #2033 --- src/modes/OptimInterface.jl | 62 ++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index 8c174460b6..e8a30f6e67 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -11,16 +11,16 @@ import Printf """ ModeResult{ - V<:NamedArrays.NamedArray, - M<:NamedArrays.NamedArray, - O<:Optim.MultivariateOptimizationResults, + V<:NamedArrays.NamedArray, + M<:NamedArrays.NamedArray, + O<:Optim.MultivariateOptimizationResults, S<:NamedArrays.NamedArray } A wrapper struct to store various results from a MAP or MLE estimation. """ struct ModeResult{ - V<:NamedArrays.NamedArray, + V<:NamedArrays.NamedArray, O<:Optim.MultivariateOptimizationResults, M<:OptimLogDensity } <: StatsBase.StatisticalModel @@ -52,7 +52,7 @@ end function StatsBase.coeftable(m::ModeResult) # Get columns for coeftable. - terms = StatsBase.coefnames(m) + terms = String.(StatsBase.coefnames(m)) estimates = m.values.array[:,1] stderrors = StatsBase.stderror(m) tstats = estimates ./ stderrors @@ -113,7 +113,7 @@ mle = optimize(model, MLE()) mle = optimize(model, MLE(), NelderMead()) ``` """ -function Optim.optimize(model::Model, ::MLE, options::Optim.Options=Optim.Options(); kwargs...) +function Optim.optimize(model::Model, ::MLE, options::Optim.Options=Optim.Options(); kwargs...) return _mle_optimize(model, options; kwargs...) end function Optim.optimize(model::Model, ::MLE, init_vals::AbstractArray, options::Optim.Options=Optim.Options(); kwargs...) @@ -123,11 +123,11 @@ function Optim.optimize(model::Model, ::MLE, optimizer::Optim.AbstractOptimizer, return _mle_optimize(model, optimizer, options; kwargs...) end function Optim.optimize( - model::Model, - ::MLE, - init_vals::AbstractArray, - optimizer::Optim.AbstractOptimizer, - options::Optim.Options=Optim.Options(); + model::Model, + ::MLE, + init_vals::AbstractArray, + optimizer::Optim.AbstractOptimizer, + options::Optim.Options=Optim.Options(); kwargs... ) return _mle_optimize(model, init_vals, optimizer, options; kwargs...) @@ -159,7 +159,7 @@ map_est = optimize(model, MAP(), NelderMead()) ``` """ -function Optim.optimize(model::Model, ::MAP, options::Optim.Options=Optim.Options(); kwargs...) +function Optim.optimize(model::Model, ::MAP, options::Optim.Options=Optim.Options(); kwargs...) return _map_optimize(model, options; kwargs...) end function Optim.optimize(model::Model, ::MAP, init_vals::AbstractArray, options::Optim.Options=Optim.Options(); kwargs...) @@ -169,11 +169,11 @@ function Optim.optimize(model::Model, ::MAP, optimizer::Optim.AbstractOptimizer, return _map_optimize(model, optimizer, options; kwargs...) end function Optim.optimize( - model::Model, - ::MAP, - init_vals::AbstractArray, - optimizer::Optim.AbstractOptimizer, - options::Optim.Options=Optim.Options(); + model::Model, + ::MAP, + init_vals::AbstractArray, + optimizer::Optim.AbstractOptimizer, + options::Optim.Options=Optim.Options(); kwargs... ) return _map_optimize(model, init_vals, optimizer, options; kwargs...) @@ -190,43 +190,43 @@ end Estimate a mode, i.e., compute a MLE or MAP estimate. """ function _optimize( - model::Model, - f::OptimLogDensity, + model::Model, + f::OptimLogDensity, optimizer::Optim.AbstractOptimizer = Optim.LBFGS(), - args...; + args...; kwargs... ) return _optimize(model, f, DynamicPPL.getparams(f), optimizer, args...; kwargs...) end function _optimize( - model::Model, - f::OptimLogDensity, + model::Model, + f::OptimLogDensity, options::Optim.Options = Optim.Options(), - args...; + args...; kwargs... ) return _optimize(model, f, DynamicPPL.getparams(f), Optim.LBFGS(), args...; kwargs...) end function _optimize( - model::Model, - f::OptimLogDensity, + model::Model, + f::OptimLogDensity, init_vals::AbstractArray = DynamicPPL.getparams(f), options::Optim.Options = Optim.Options(), - args...; + args...; kwargs... ) return _optimize(model, f, init_vals, Optim.LBFGS(), options, args...; kwargs...) end function _optimize( - model::Model, - f::OptimLogDensity, - init_vals::AbstractArray = DynamicPPL.getparams(f), + model::Model, + f::OptimLogDensity, + init_vals::AbstractArray = DynamicPPL.getparams(f), optimizer::Optim.AbstractOptimizer = Optim.LBFGS(), options::Optim.Options = Optim.Options(), - args...; + args...; kwargs... ) # Convert the initial values, since it is assumed that users provide them @@ -243,7 +243,7 @@ function _optimize( @warn "Optimization did not converge! You may need to correct your model or adjust the Optim parameters." end - # Get the VarInfo at the MLE/MAP point, and run the model to ensure + # Get the VarInfo at the MLE/MAP point, and run the model to ensure # correct dimensionality. @set! f.varinfo = DynamicPPL.unflatten(f.varinfo, M.minimizer) @set! f.varinfo = invlink!!(f.varinfo, model) From f4c6876888f9b82c6371f2245847ca54a02beb24 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 10 Jul 2023 16:14:41 +0100 Subject: [PATCH 02/10] Propose additional columns --- src/modes/OptimInterface.jl | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index e8a30f6e67..f2556c8ab2 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -25,13 +25,13 @@ struct ModeResult{ M<:OptimLogDensity } <: StatsBase.StatisticalModel "A vector with the resulting point estimates." - values :: V + values::V "The stored Optim.jl results." - optim_result :: O + optim_result::O "The final log likelihood or log joint, depending on whether `MAP` or `MLE` was run." - lp :: Float64 + lp::Float64 "The evaluation function used to calculate the output." - f :: M + f::M end ############################# # Various StatsBase methods # @@ -53,11 +53,16 @@ end function StatsBase.coeftable(m::ModeResult) # Get columns for coeftable. terms = String.(StatsBase.coefnames(m)) - estimates = m.values.array[:,1] + estimates = m.values.array[:, 1] stderrors = StatsBase.stderror(m) - tstats = estimates ./ stderrors + zscore = estimates ./ stderrors + # p = 2 * (1 .- cdf.(Normal(0, 1), abs.(zscore))) - StatsBase.CoefTable([estimates, stderrors, tstats], ["estimate", "stderror", "tstat"], terms) + # # 95% confidence interval (CI) + # ci_low = estimates .+ quantile(Normal(0, 1), 0.025) .* stderrors + # ci_high = estimates .+ quantile(Normal(0, 1), 0.975) .* stderrors + + StatsBase.CoefTable([estimates, stderrors, zscore], ["estimate", "stderror", "tstat"], terms) end function StatsBase.informationmatrix(m::ModeResult; hessian_function=ForwardDiff.hessian, kwargs...) @@ -192,7 +197,7 @@ Estimate a mode, i.e., compute a MLE or MAP estimate. function _optimize( model::Model, f::OptimLogDensity, - optimizer::Optim.AbstractOptimizer = Optim.LBFGS(), + optimizer::Optim.AbstractOptimizer=Optim.LBFGS(), args...; kwargs... ) @@ -202,7 +207,7 @@ end function _optimize( model::Model, f::OptimLogDensity, - options::Optim.Options = Optim.Options(), + options::Optim.Options=Optim.Options(), args...; kwargs... ) @@ -212,8 +217,8 @@ end function _optimize( model::Model, f::OptimLogDensity, - init_vals::AbstractArray = DynamicPPL.getparams(f), - options::Optim.Options = Optim.Options(), + init_vals::AbstractArray=DynamicPPL.getparams(f), + options::Optim.Options=Optim.Options(), args...; kwargs... ) @@ -223,9 +228,9 @@ end function _optimize( model::Model, f::OptimLogDensity, - init_vals::AbstractArray = DynamicPPL.getparams(f), - optimizer::Optim.AbstractOptimizer = Optim.LBFGS(), - options::Optim.Options = Optim.Options(), + init_vals::AbstractArray=DynamicPPL.getparams(f), + optimizer::Optim.AbstractOptimizer=Optim.LBFGS(), + options::Optim.Options=Optim.Options(), args...; kwargs... ) From ba7d7f34caccdba8969df63d7393cb37a0133fb3 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 10 Jul 2023 17:09:28 +0100 Subject: [PATCH 03/10] Add other suggestions --- src/modes/OptimInterface.jl | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index f2556c8ab2..fb8d05687c 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -50,19 +50,23 @@ function Base.show(io::IO, m::ModeResult) show(io, m.values.array) end -function StatsBase.coeftable(m::ModeResult) +function StatsBase.coeftable(m::ModeResult; level::Real=0.95) # Get columns for coeftable. terms = String.(StatsBase.coefnames(m)) estimates = m.values.array[:, 1] stderrors = StatsBase.stderror(m) zscore = estimates ./ stderrors - # p = 2 * (1 .- cdf.(Normal(0, 1), abs.(zscore))) + p = 2 * (1 .- cdf.(Normal(0, 1), abs.(zscore))) - # # 95% confidence interval (CI) - # ci_low = estimates .+ quantile(Normal(0, 1), 0.025) .* stderrors - # ci_high = estimates .+ quantile(Normal(0, 1), 0.975) .* stderrors + # Confidence interval (CI) + q = quantile(Normal(0, 1), 1 - (1 - level) / 2) + ci_low = estimates .- q .* stderrors + ci_high = estimates .+ q .* stderrors - StatsBase.CoefTable([estimates, stderrors, zscore], ["estimate", "stderror", "tstat"], terms) + StatsBase.CoefTable( + [estimates, stderrors, zscore, p, ci_low, ci_high], + ["Coef.", "Std. Error", "z", "Pr(>|z|)", "Lower 95%", "Upper 95%"], + terms) end function StatsBase.informationmatrix(m::ModeResult; hessian_function=ForwardDiff.hessian, kwargs...) From f5dec2ee5ee38092dfd5dadd070ca65d113737e4 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 10 Jul 2023 18:22:22 +0100 Subject: [PATCH 04/10] Update src/modes/OptimInterface.jl Co-authored-by: David Widmann --- src/modes/OptimInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index fb8d05687c..9bdd33e261 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -56,7 +56,7 @@ function StatsBase.coeftable(m::ModeResult; level::Real=0.95) estimates = m.values.array[:, 1] stderrors = StatsBase.stderror(m) zscore = estimates ./ stderrors - p = 2 * (1 .- cdf.(Normal(0, 1), abs.(zscore))) + p = pvalue(Normal(), zscore; tail=:both) # Confidence interval (CI) q = quantile(Normal(0, 1), 1 - (1 - level) / 2) From 52f460510763aa5013993796e02bb6ea0cd13c90 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 10 Jul 2023 18:22:35 +0100 Subject: [PATCH 05/10] Update src/modes/OptimInterface.jl Co-authored-by: David Widmann --- src/modes/OptimInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index 9bdd33e261..db29422b4b 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -52,7 +52,7 @@ end function StatsBase.coeftable(m::ModeResult; level::Real=0.95) # Get columns for coeftable. - terms = String.(StatsBase.coefnames(m)) + terms = string.(StatsBase.coefnames(m)) estimates = m.values.array[:, 1] stderrors = StatsBase.stderror(m) zscore = estimates ./ stderrors From e5d7774cbe78bfdaa44a0c0801be3714f4eb41c4 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Mon, 10 Jul 2023 18:22:43 +0100 Subject: [PATCH 06/10] Update src/modes/OptimInterface.jl Co-authored-by: David Widmann --- src/modes/OptimInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index db29422b4b..a47458a24b 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -59,7 +59,7 @@ function StatsBase.coeftable(m::ModeResult; level::Real=0.95) p = pvalue(Normal(), zscore; tail=:both) # Confidence interval (CI) - q = quantile(Normal(0, 1), 1 - (1 - level) / 2) + q = quantile(Normal(), (1 + level) / 2) ci_low = estimates .- q .* stderrors ci_high = estimates .+ q .* stderrors From aae9d43c93c99996308de187799055b113ecb9b7 Mon Sep 17 00:00:00 2001 From: Hong Ge <3279477+yebai@users.noreply.github.com> Date: Mon, 10 Jul 2023 19:47:09 +0100 Subject: [PATCH 07/10] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 198739423d..5b271c5c06 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Turing" uuid = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" -version = "0.26.3" +version = "0.26.4" [deps] AbstractMCMC = "80f14c24-f653-4e6a-9b94-39d6b0f70001" From 41e133129a46fbd466d48cb8e740c14fda857ff8 Mon Sep 17 00:00:00 2001 From: Dominique Makowski Date: Tue, 11 Jul 2023 07:28:11 +0100 Subject: [PATCH 08/10] import SsStatsAPI --- src/modes/OptimInterface.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index a47458a24b..52b582ec61 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -7,6 +7,7 @@ import ..ForwardDiff import NamedArrays import StatsBase import Printf +import StatsAPI """ @@ -56,7 +57,7 @@ function StatsBase.coeftable(m::ModeResult; level::Real=0.95) estimates = m.values.array[:, 1] stderrors = StatsBase.stderror(m) zscore = estimates ./ stderrors - p = pvalue(Normal(), zscore; tail=:both) + p = StatsAPI.pvalue(Normal(), zscore; tail=:both) # Confidence interval (CI) q = quantile(Normal(), (1 + level) / 2) From db2f9b41251da52b6cc5edbd830c0611a5d59bad Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 11 Jul 2023 09:35:04 +0200 Subject: [PATCH 09/10] Update Project.toml --- Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 5b271c5c06..a4cc4686aa 100644 --- a/Project.toml +++ b/Project.toml @@ -32,6 +32,7 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" @@ -62,6 +63,7 @@ Requires = "0.5, 1.0" SciMLBase = "1.37.1" Setfield = "0.8, 1" SpecialFunctions = "0.7.2, 0.8, 0.9, 0.10, 1, 2" +StatsAPI = "1.6" StatsBase = "0.32, 0.33, 0.34" StatsFuns = "0.8, 0.9, 1" Tracker = "0.2.3" From b27b3940fc8b8453023ba3ae63af78a58759e1ce Mon Sep 17 00:00:00 2001 From: Hong Ge <3279477+yebai@users.noreply.github.com> Date: Tue, 11 Jul 2023 22:51:45 +0100 Subject: [PATCH 10/10] Update src/modes/OptimInterface.jl Co-authored-by: David Widmann --- src/modes/OptimInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modes/OptimInterface.jl b/src/modes/OptimInterface.jl index 52b582ec61..f8c288d03b 100644 --- a/src/modes/OptimInterface.jl +++ b/src/modes/OptimInterface.jl @@ -57,7 +57,7 @@ function StatsBase.coeftable(m::ModeResult; level::Real=0.95) estimates = m.values.array[:, 1] stderrors = StatsBase.stderror(m) zscore = estimates ./ stderrors - p = StatsAPI.pvalue(Normal(), zscore; tail=:both) + p = map(z -> StatsAPI.pvalue(Normal(), z; tail=:both), zscore) # Confidence interval (CI) q = quantile(Normal(), (1 + level) / 2)