Merge pull request #79 from DanielVandH/refine

Refinement methods

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "ProfileLikelihood"
 uuid = "3fca794e-44fa-49a6-b6a0-8cd45572ba0a"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

docs/src/docstrings.md

@@ -24,7 +24,10 @@ ProfileLikelihood.ProfileLikelihoodSolution
 ProfileLikelihood.ConfidenceInterval
 profile 
 replace_profile!
+refine_profile!
 ProfileLikelihood.set_next_initial_estimate!(::Any, ::Any, ::Any, ::Any, ::Any)
+ProfileLikelihood.get_confidence_intervals!
+ProfileLikelihood.reach_min_steps!
 ```
 
 ## BivariateProfileLikelihoodSolution 

src/ProfileLikelihood.jl

@@ -54,6 +54,7 @@ export get_range
 export replace_profile!
 export construct_profile_grids
 export bivariate_profile
+export refine_profile!
 
 function __init__()
     @require CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" begin

src/bivariate.jl

@@ -55,10 +55,10 @@ function bivariate_profile(prob::LikelihoodProblem, sol::LikelihoodSolution, n::
     parallel=Val(false),
     next_initial_estimate_method=Val(:nearest),
     kwargs...) where {M,F}
-    parallel = Val(SciMLBase._unwrap_val(parallel))
-    normalise = Val(SciMLBase._unwrap_val(normalise))
-    confidence_region_method = Val(SciMLBase._unwrap_val(confidence_region_method))
-    next_initial_estimate_method = Val(SciMLBase._unwrap_val(next_initial_estimate_method))
+    parallel = _Val(parallel)
+    normalise = _Val(normalise)
+    confidence_region_method = _Val(confidence_region_method)
+    next_initial_estimate_method = _Val(next_initial_estimate_method)
 
     ## Extract the problem and solution 
     opt_prob, mles, ℓmax = extract_problem_and_solution(prob, sol)
@@ -173,7 +173,7 @@ end
 @inline function prepare_cache_vectors(n, mles::AbstractVector{T}, res, num_params, normalise, ℓmax; parallel=Val(false)) where {T}
     profile_vals = OffsetArray(zeros(T, 2res + 1, 2res + 1), -res:res, -res:res)
     other_mles = OffsetArray([zeros(T, num_params - 2) for _ in 1:(2res+1), _ in 1:(2res+1)], -res:res, -res:res)
-    normalise = SciMLBase._unwrap_val(normalise)
+    normalise = take_val(normalise)
     profile_vals[0, 0] = normalise ? zero(T) : ℓmax
     other_mles[0, 0] .= mles[Not(n[1], n[2])]
     if parallel == Val(false)
@@ -264,7 +264,7 @@ end
     fixed_prob = setup_layer_node_problem!(fixed_vals, grid, I, sub_cache, other_mle, layer, restricted_prob,
         next_initial_estimate_method, cache, n)
     soln = solve(fixed_prob, alg; kwargs...)
-    normalise = SciMLBase._unwrap_val(normalise)
+    normalise = take_val(normalise)
     @inbounds profile_vals[I] = -soln.objective - ℓmax * !normalise
     @inbounds other_mle[I] = soln.u
     return compare_to_threshold(any_above_threshold, profile_vals, threshold, I)

src/plotting.jl

@@ -10,6 +10,8 @@ const ALPHABET = join('a':'z')
         shade_ci=true, 
         fig_kwargs=nothing, 
         axis_kwargs=nothing,
+        show_points=false,
+        markersize=9,
         latex_names = Dict(vars .=> [LaTeXStrings.L"\theta_{i}" for i in SciMLBase.sym_to_index.(vars, Ref(prof))])) 
      
 Plot results from a profile likelihood solution `prof`.
@@ -27,6 +29,8 @@ Plot results from a profile likelihood solution `prof`.
 - `shade_ci=true`: Whether to shade the area under the profile between the confidence interval.
 - `fig_kwargs=nothing`: Extra keyword arguments for `Figure` (see the Makie docs).
 - `axis_kwargs=nothing`: Extra keyword arguments for `Axis` (see the Makie docs).
+- `show_points=false`: Whether to show the profile data. 
+- `markersize=9`: The marker size used for `show_points`.
 - `latex_names = Dict(vars .=> [LaTeXStrings.L"\theta_{i}" for i in SciMLBase.sym_to_index.(vars, Ref(prof))]))`: LaTeX names to use for the parameters. Defaults to `θᵢ`, where `i` is the index of the parameter. 
 
 # Output 
@@ -41,6 +45,8 @@ function plot_profiles(prof::ProfileLikelihoodSolution, vars=profiled_parameters
     shade_ci=true,
     fig_kwargs=nothing,
     axis_kwargs=nothing,
+    show_points=false,
+    markersize=9,
     latex_names=Dict(vars .=> [LaTeXStrings.L"\theta_{%$i}" for i in SciMLBase.sym_to_index.(vars, Ref(prof))]))
     num_plots = vars isa Symbol ? 1 : length(vars)
     _, _, plot_positions = choose_grid_layout(num_plots, ncol, nrow)
@@ -53,9 +59,9 @@ function plot_profiles(prof::ProfileLikelihoodSolution, vars=profiled_parameters
     for (ℓ, k) in itr
         i, j = plot_positions[ℓ]
         if axis_kwargs !== nothing
-            plot_profile!(prof[k], fig, ℓ, k, i, j, spline, true_vals[k], show_mles ? get_likelihood_solution(prof)[k] : nothing, shade_ci, latex_names[k]; axis_kwargs)
+            plot_profile!(prof[k], fig, ℓ, k, i, j, spline, true_vals[k], show_mles ? get_likelihood_solution(prof)[k] : nothing, shade_ci, latex_names[k], show_points, markersize; axis_kwargs)
         else
-            plot_profile!(prof[k], fig, ℓ, k, i, j, spline, true_vals[k], show_mles ? get_likelihood_solution(prof)[k] : nothing, shade_ci, latex_names[k])
+            plot_profile!(prof[k], fig, ℓ, k, i, j, spline, true_vals[k], show_mles ? get_likelihood_solution(prof)[k] : nothing, shade_ci, latex_names[k], show_points, markersize)
         end
     end
     return fig
@@ -150,7 +156,8 @@ function choose_grid_layout(num_plots, cols, rows)
 end
 
 function plot_profile!(prof::ProfileLikelihoodSolutionView, fig, ℓ, k, i, j,
-    spline, true_vals, mle_val=nothing, shade_ci=true, param_name=LaTeXStrings.L"\theta_{%$ℓ}"; axis_kwargs=nothing)
+    spline, true_vals, mle_val=nothing, shade_ci=true, param_name=LaTeXStrings.L"\theta_{%$ℓ}",
+    show_points=false, markersize=9; axis_kwargs=nothing)
     lower_ci, upper_ci = get_confidence_intervals(prof)
     θ_vals = get_parameter_values(prof)
     ℓ_vals = get_profile_values(prof)
@@ -196,6 +203,9 @@ function plot_profile!(prof::ProfileLikelihoodSolutionView, fig, ℓ, k, i, j,
     if !isnothing(mle_val)
         CairoMakie.vlines!(ax, [mle_val], color=:red, linetype=:dashed)
     end
+    if show_points
+        CairoMakie.scatter!(ax, θ_vals, ℓ_vals, color=:black, markersize=markersize)
+    end
     return nothing
 end
 function plot_profile!(prof::BivariateProfileLikelihoodSolutionView, fig, ℓ, (k, r), i, j,

src/prediction_intervals.jl

@@ -173,17 +173,25 @@ end
     return build_q_prototype(q, get_parent(prof), data)
 end
 
+@inline function spline_other_mles(parameter_values, other_mles)
+    spline = interpolate((parameter_values,), other_mles, Gridded(Linear()))
+    return spline
+end
+@inline function spline_other_mles(parameter_values_1, parameter_values_2, other_mles)
+    spline = interpolate((parameter_values_1, parameter_values_2), other_mles, Gridded(Linear()))
+    return spline
+end
 @inline function spline_other_mles(prof::ProfileLikelihoodSolutionView)
     data = get_parameter_values(prof)
     other_mles = get_other_mles(prof)
-    spline = interpolate((data,), other_mles, Gridded(Linear()))
+    spline = spline_other_mles(data, other_mles)
     return spline
 end
 @inline function spline_other_mles(prof::BivariateProfileLikelihoodSolutionView)
     grid_1 = get_parameter_values(prof, 1)
     grid_2 = get_parameter_values(prof, 2)
     other_mles = get_other_mles(prof)
-    spline = interpolate((grid_1, grid_2), other_mles, Gridded(Linear()))
+    spline = spline_other_mles(grid_1, grid_2, other_mles)
     return spline
 end
 @inline function spline_other_mles(prof::Union{ProfileLikelihoodSolution,BivariateProfileLikelihoodSolution})