From c44d8d81d9142a8feac346803f659e7d0d78d6fb Mon Sep 17 00:00:00 2001 From: Maximilian Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> Date: Tue, 15 Apr 2025 18:01:50 +0200 Subject: [PATCH 1/4] Documentation/fix docstrings (#260) Clean code, fix docstrings, fix docs --- Project.toml | 2 +- README.md | 2 +- docs/make.jl | 1 - docs/src/developer/implied.md | 14 ++-- docs/src/developer/loss.md | 2 +- docs/src/developer/observed.md | 4 +- docs/src/developer/optimizer.md | 4 +- docs/src/internals/files.md | 2 +- docs/src/internals/internals.md | 2 +- docs/src/performance/simulation.md | 4 -- docs/src/tutorials/backends/nlopt.md | 2 +- .../tutorials/construction/build_by_parts.md | 10 +-- .../construction/outer_constructor.md | 5 +- .../regularization/regularization.md | 3 +- src/additional_functions/helper.jl | 67 ++----------------- src/additional_functions/simulation.jl | 2 +- .../start_val/start_fabin3.jl | 21 +----- .../start_val/start_simple.jl | 1 - src/frontend/common.jl | 4 +- src/frontend/fit/fitmeasures/fit_measures.jl | 2 +- src/frontend/fit/fitmeasures/minus2ll.jl | 2 +- src/frontend/fit/summary.jl | 12 +--- src/frontend/pretty_printing.jl | 6 +- src/frontend/specification/RAMMatrices.jl | 2 +- src/frontend/specification/Sem.jl | 10 --- src/frontend/specification/StenoGraphs.jl | 1 - src/frontend/specification/documentation.jl | 27 ++++---- src/implied/RAM/generic.jl | 40 +++++------ src/implied/RAM/symbolic.jl | 30 ++++----- src/implied/abstract.jl | 1 - src/implied/empty.jl | 7 +- src/loss/ML/FIML.jl | 4 -- src/loss/ML/ML.jl | 7 +- src/loss/WLS/WLS.jl | 9 +-- src/loss/constant/constant.jl | 4 -- src/loss/regularization/ridge.jl | 4 -- src/objective_gradient_hessian.jl | 11 +-- src/observed/EM.jl | 54 +++++++-------- src/observed/abstract.jl | 2 +- src/observed/covariance.jl | 12 ++-- src/observed/data.jl | 13 +--- src/observed/missing.jl | 16 ++--- src/optimizer/Empty.jl | 6 -- src/optimizer/abstract.jl | 14 +++- src/optimizer/optim.jl | 4 +- src/package_extensions/SEMNLOptExt.jl | 4 -- src/package_extensions/SEMProximalOptExt.jl | 10 ++- src/types.jl | 6 +- test/unit_tests/bootstrap.jl | 5 -- test/unit_tests/partable_dataframe.jl | 3 - test/unit_tests/ridge_id.jl | 8 --- test/unit_tests/sorting.jl | 17 ----- test/unit_tests/start_val.jl | 1 - test/unit_tests/unit_tests.jl | 8 --- 54 files changed, 167 insertions(+), 347 deletions(-) delete mode 100644 test/unit_tests/bootstrap.jl delete mode 100644 test/unit_tests/partable_dataframe.jl delete mode 100644 test/unit_tests/ridge_id.jl delete mode 100644 test/unit_tests/sorting.jl delete mode 100644 test/unit_tests/start_val.jl diff --git a/Project.toml b/Project.toml index 9a5bc0916..ae0fa5ba4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StructuralEquationModels" uuid = "383ca8c5-e4ff-4104-b0a9-f7b279deed53" authors = ["Maximilian Ernst", "Aaron Peikert"] -version = "0.4.0" +version = "0.4.1" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" diff --git a/README.md b/README.md index 7fc19d0f8..1c3714546 100644 --- a/README.md +++ b/README.md @@ -39,7 +39,7 @@ The package makes use of - SparseArrays.jl to speed up symbolic computations. - Optim.jl and NLopt.jl to provide a range of different Optimizers/Linesearches. - ProximalAlgorithms.jl for regularization. -- FiniteDiff.jl and ForwardDiff.jl to provide gradients for user-defined loss functions. +- FiniteDiff.jl and to provide gradient approximations for user-defined loss functions. # At the moment, we are still working on: - optimizing performance for big models (with hundreds of parameters) diff --git a/docs/make.jl b/docs/make.jl index 4542cf48f..1bb68c4da 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -52,7 +52,6 @@ makedocs( "files" => "internals/files.md", "types" => "internals/types.md", ], - "Complementary material" => ["Mathematical appendix" => "complementary/maths.md"], ], format = Documenter.HTML( prettyurls = get(ENV, "CI", nothing) == "true", diff --git a/docs/src/developer/implied.md b/docs/src/developer/implied.md index bea824a94..056cd6638 100644 --- a/docs/src/developer/implied.md +++ b/docs/src/developer/implied.md @@ -62,11 +62,8 @@ model per group and an additional model with `ImpliedEmpty` and `SemRidge` for t # Extended help ## Interfaces -- `params(::RAMSymbolic) `-> Vector of parameter labels -- `nparams(::RAMSymbolic)` -> Number of parameters - -## Implementation -Subtype of `SemImplied`. +- `param_labels(::ImpliedEmpty) `-> Vector of parameter labels +- `nparams(::ImpliedEmpty)` -> Number of parameters """ struct ImpliedEmpty{A, B, C} <: SemImplied hessianeval::A @@ -78,7 +75,12 @@ end ### Constructors ############################################################################################ -function ImpliedEmpty(;specification, meanstruct = NoMeanStruct(), hessianeval = ExactHessian(), kwargs...) +function ImpliedEmpty(; + specification, + meanstruct = NoMeanStruct(), + hessianeval = ExactHessian(), + kwargs..., +) return ImpliedEmpty(hessianeval, meanstruct, convert(RAMMatrices, specification)) end diff --git a/docs/src/developer/loss.md b/docs/src/developer/loss.md index 931c2d0e5..d6949842b 100644 --- a/docs/src/developer/loss.md +++ b/docs/src/developer/loss.md @@ -204,7 +204,7 @@ julia>? help?> RAM -help?> SemObservedCommon +help?> SemObservedData ``` We see that the model implied covariance matrix can be assessed as `Σ(implied)` and the observed covariance matrix as `obs_cov(observed)`. diff --git a/docs/src/developer/observed.md b/docs/src/developer/observed.md index 240c1c34f..ee795cd3e 100644 --- a/docs/src/developer/observed.md +++ b/docs/src/developer/observed.md @@ -28,8 +28,8 @@ nsamples(observed::MyObserved) = ... nobserved_vars(observed::MyObserved) = ... ``` -As always, you can add additional methods for properties that implied types and loss function want to access, for example (from the `SemObservedCommon` implementation): +As always, you can add additional methods for properties that implied types and loss function want to access, for example (from the `SemObservedData` implementation): ```julia -obs_cov(observed::SemObservedCommon) = observed.obs_cov +obs_cov(observed::SemObservedData) = observed.obs_cov ``` \ No newline at end of file diff --git a/docs/src/developer/optimizer.md b/docs/src/developer/optimizer.md index a651ec636..0a85eb6a6 100644 --- a/docs/src/developer/optimizer.md +++ b/docs/src/developer/optimizer.md @@ -37,13 +37,13 @@ options(optimizer::SemOptimizerName) = optimizer.options Note that your optimizer is a subtype of `SemOptimizer{:Name}`, where you can choose a `:Name` that can later be used as a keyword argument to `fit(engine = :Name)`. Similarly, `SemOptimizer{:Name}(args...; kwargs...) = SemOptimizerName(args...; kwargs...)` should be defined as well as a constructor that uses only keyword arguments: -´´´julia +```julia SemOptimizerName(; algorithm = LBFGS(), options = Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8), kwargs..., ) = SemOptimizerName(algorithm, options) -´´´ +``` A method for `update_observed` and additional methods might be usefull, but are not necessary. Now comes the substantive part: We need to provide a method for `fit`: diff --git a/docs/src/internals/files.md b/docs/src/internals/files.md index 90ceceaaf..4c2338393 100644 --- a/docs/src/internals/files.md +++ b/docs/src/internals/files.md @@ -10,7 +10,7 @@ Source code is in the `"src"` folder: - `"StructuralEquationModels.jl"` defines the module and the exported objects - `"types.jl"` defines all abstract types and the basic type hierarchy - `"objective_gradient_hessian.jl"` contains methods for computing objective, gradient and hessian values for different model types as well as generic fallback methods -- The four folders `"observed"`, `"implied"`, `"loss"` and `"diff"` contain implementations of specific subtypes (for example, the `"loss"` folder contains a file `"ML.jl"` that implements the `SemML` loss function). +- The folders `"observed"`, `"implied"`, and `"loss"` contain implementations of specific subtypes (for example, the `"loss"` folder contains a file `"ML.jl"` that implements the `SemML` loss function). - `"optimizer"` contains connections to different optimization backends (aka methods for `fit`) - `"optim.jl"`: connection to the `Optim.jl` package - `"frontend"` contains user-facing functions diff --git a/docs/src/internals/internals.md b/docs/src/internals/internals.md index 77d1085c8..18f82dbca 100644 --- a/docs/src/internals/internals.md +++ b/docs/src/internals/internals.md @@ -1,3 +1,3 @@ # Internals and Design -On the following pages, we document the internals and design of the package. Those informations are no prerequisite for extending the package (as decribed in the developer documentation)!, but they may be useful and hopefully interesting. \ No newline at end of file +On the following pages, we document some technical information about the package. Those informations are no prerequisite for extending the package (as decribed in the developer documentation)!, but they may be useful. \ No newline at end of file diff --git a/docs/src/performance/simulation.md b/docs/src/performance/simulation.md index 0cb2ea25d..d268853f6 100644 --- a/docs/src/performance/simulation.md +++ b/docs/src/performance/simulation.md @@ -1,9 +1,5 @@ # Simulation studies -!!! note "Simulation study interface" - We are currently working on an interface for simulation studies. - Until we are finished with this, this page is just a collection of tips. - ## Replace observed data In simulation studies, a common task is fitting the same model to many different datasets. It would be a waste of resources to reconstruct the complete model for each dataset. diff --git a/docs/src/tutorials/backends/nlopt.md b/docs/src/tutorials/backends/nlopt.md index 2afa5e547..feb5c8f48 100644 --- a/docs/src/tutorials/backends/nlopt.md +++ b/docs/src/tutorials/backends/nlopt.md @@ -1,7 +1,7 @@ # Using NLopt.jl [`SemOptimizerNLopt`](@ref) implements the connection to `NLopt.jl`. -It is only available if the `NLopt` package is loaded alongside `StructuralEquationModel.jl` in the running Julia session. +It is only available if the `NLopt` package is loaded alongside `StructuralEquationModels.jl` in the running Julia session. It takes a bunch of arguments: ```julia diff --git a/docs/src/tutorials/construction/build_by_parts.md b/docs/src/tutorials/construction/build_by_parts.md index 45d2a2ea1..6b6b59ac9 100644 --- a/docs/src/tutorials/construction/build_by_parts.md +++ b/docs/src/tutorials/construction/build_by_parts.md @@ -47,21 +47,21 @@ partable = ParameterTable( Now, we construct the different parts: ```@example build -# observed --------------------------------------------------------------------------------- +# observed ----------------------------------------------------------------------------- observed = SemObservedData(specification = partable, data = data) -# implied ------------------------------------------------------------------------------------ +# implied ------------------------------------------------------------------------------ implied_ram = RAM(specification = partable) -# loss ------------------------------------------------------------------------------------- +# loss --------------------------------------------------------------------------------- ml = SemML(observed = observed) loss_ml = SemLoss(ml) -# optimizer ------------------------------------------------------------------------------------- +# optimizer ---------------------------------------------------------------------------- optimizer = SemOptimizerOptim() -# model ------------------------------------------------------------------------------------ +# model -------------------------------------------------------------------------------- model_ml = Sem(observed, implied_ram, loss_ml) diff --git a/docs/src/tutorials/construction/outer_constructor.md b/docs/src/tutorials/construction/outer_constructor.md index a1c0b8ad3..e27724307 100644 --- a/docs/src/tutorials/construction/outer_constructor.md +++ b/docs/src/tutorials/construction/outer_constructor.md @@ -14,9 +14,8 @@ Structural Equation Model - Loss Functions SemML - Fields - observed: SemObservedCommon - implied: RAM - optimizer: SemOptimizerOptim + observed: SemObservedData + implied: RAM ``` The output of this call tells you exactly what model you just constructed (i.e. what the loss functions, observed, implied and optimizer parts are). diff --git a/docs/src/tutorials/regularization/regularization.md b/docs/src/tutorials/regularization/regularization.md index b1fc5c157..f743ac79c 100644 --- a/docs/src/tutorials/regularization/regularization.md +++ b/docs/src/tutorials/regularization/regularization.md @@ -30,13 +30,12 @@ It can be used as ```julia SemOptimizerProximal( algorithm = ProximalAlgorithms.PANOC(), - options = Dict{Symbol, Any}(), operator_g, operator_h = nothing ) ``` -The proximal operator (aka the regularization function) can be passed as `operator_g`, available options are listed [here](https://juliafirstorder.github.io/ProximalOperators.jl/stable/functions/). +The proximal operator (aka the regularization function) can be passed as `operator_g`. The available Algorithms are listed [here](https://juliafirstorder.github.io/ProximalAlgorithms.jl/stable/guide/implemented_algorithms/). ## First example - lasso diff --git a/src/additional_functions/helper.jl b/src/additional_functions/helper.jl index 5559034e0..6cbcb0573 100644 --- a/src/additional_functions/helper.jl +++ b/src/additional_functions/helper.jl @@ -14,84 +14,29 @@ function neumann_series(mat::SparseMatrixCSC; maxiter::Integer = size(mat, 1)) return inverse end -#= -function make_onelement_array(A) - isa(A, Array) ? nothing : (A = [A]) - return A -end - =# - -function semvec(observed, implied, loss, optimizer) - observed = make_onelement_array(observed) - implied = make_onelement_array(implied) - loss = make_onelement_array(loss) - optimizer = make_onelement_array(optimizer) - - #sem_vec = Array{AbstractSem}(undef, maximum(length.([observed, implied, loss, optimizer]))) - sem_vec = Sem.(observed, implied, loss, optimizer) - - return sem_vec -end - -skipmissing_mean(mat::AbstractMatrix) = - [mean(skipmissing(coldata)) for coldata in eachcol(mat)] - -function F_one_person(imp_mean, meandiff, inverse, data, logdet) - F = logdet - @. meandiff = data - imp_mean - F += dot(meandiff, inverse, meandiff) - return F -end - -function remove_all_missing(data::AbstractMatrix) - keep = Vector{Int64}() - for (i, coldata) in zip(axes(data, 1), eachrow(data)) - if any(!ismissing, coldata) - push!(keep, i) - end - end - return data[keep, :], keep -end - function batch_inv!(fun, model) for i in 1:size(fun.inverses, 1) fun.inverses[i] .= LinearAlgebra.inv!(fun.choleskys[i]) end end -#= -function batch_sym_inv_update!(fun::Union{LossFunction, DiffFunction}, model) - M_inv = inv(fun.choleskys[1]) - for i = 1:size(fun.inverses, 1) - if size(model.observed.patterns_not[i]) == 0 - fun.inverses[i] .= M_inv - else - ind_not = model.observed.patterns_not[i] - ind = model.observed.patterns[i] - - A = M_inv[ind_not, ind] - H = cholesky(M_inv[ind_not, ind_not]) - D = H \ A - out = M_inv[ind, ind] - LinearAlgebra.BLAS.gemm('T', 'N', 1.0, A, D) - fun.inverses[i] .= out - end - end -end =# - -function sparse_outer_mul!(C, A, B, ind) #computes A*S*B -> C, where ind gives the entries of S that are 1 +# computes A*S*B -> C, where ind gives the entries of S that are 1 +function sparse_outer_mul!(C, A, B, ind) fill!(C, 0.0) for i in 1:length(ind) BLAS.ger!(1.0, A[:, ind[i][1]], B[ind[i][2], :], C) end end -function sparse_outer_mul!(C, A, ind) #computes A*∇m, where ∇m ind gives the entries of ∇m that are 1 +# computes A*∇m, where ∇m ind gives the entries of ∇m that are 1 +function sparse_outer_mul!(C, A, ind) fill!(C, 0.0) @views C .= sum(A[:, ind], dims = 2) return C end -function sparse_outer_mul!(C, A, B::Vector, ind) #computes A*S*B -> C, where ind gives the entries of S that are 1 +# computes A*S*B -> C, where ind gives the entries of S that are 1 +function sparse_outer_mul!(C, A, B::Vector, ind) fill!(C, 0.0) @views @inbounds for i in 1:length(ind) C .+= B[ind[i][2]] .* A[:, ind[i][1]] diff --git a/src/additional_functions/simulation.jl b/src/additional_functions/simulation.jl index 27d58f93f..da3e6a906 100644 --- a/src/additional_functions/simulation.jl +++ b/src/additional_functions/simulation.jl @@ -7,7 +7,7 @@ Return a new model with swaped observed part. # Arguments - `model::AbstractSemSingle`: model to swap the observed part of. -- `kwargs`: additional keyword arguments; typically includes `data = ...` +- `kwargs`: additional keyword arguments; typically includes `data` and `specification` - `observed`: Either an object of subtype of `SemObserved` or a subtype of `SemObserved` # Examples diff --git a/src/additional_functions/start_val/start_fabin3.jl b/src/additional_functions/start_val/start_fabin3.jl index bd55f21d7..ab79d9ada 100644 --- a/src/additional_functions/start_val/start_fabin3.jl +++ b/src/additional_functions/start_val/start_fabin3.jl @@ -18,7 +18,7 @@ end # SemObservedMissing function start_fabin3(observed::SemObservedMissing, implied, args...; kwargs...) if !observed.em_model.fitted - em_mvn(observed; kwargs...) + em_mvn!(observed; kwargs...) end return start_fabin3(implied.ram_matrices, observed.em_model.Σ, observed.em_model.μ) @@ -45,24 +45,6 @@ function start_fabin3( ) @assert length(F_var2obs) == size(F, 1) - # check in which matrix each parameter appears - - #= in_S = length.(S_ind) .!= 0 - in_A = length.(A_ind) .!= 0 - A_ind_c = [linear2cartesian(ind, (n_var, n_var)) for ind in A_ind] - in_Λ = [any(ind[2] .∈ F_ind) for ind in A_ind_c] - - if !isnothing(M) - in_M = length.(M_ind) .!= 0 - in_any = in_A .| in_S .| in_M - else - in_any = in_A .| in_S - end - - if !all(in_any) - @warn "Could not determine fabin3 starting values for some parameters, default to 0." - end =# - # set undirected parameters in S S_indices = CartesianIndices(S) for j in 1:nparams(S) @@ -79,7 +61,6 @@ function start_fabin3( # set loadings A_indices = CartesianIndices(A) - # ind_Λ = findall([is_in_Λ(ind_vec, F_ind) for ind_vec in A_ind_c]) # collect latent variable indicators in A # maps latent parameter to the vector of dependent vars diff --git a/src/additional_functions/start_val/start_simple.jl b/src/additional_functions/start_val/start_simple.jl index ad5148e3f..4fbc8719c 100644 --- a/src/additional_functions/start_val/start_simple.jl +++ b/src/additional_functions/start_val/start_simple.jl @@ -63,7 +63,6 @@ function start_simple( nparams(ram_matrices) start_val = zeros(n_par) - n_obs = nobserved_vars(ram_matrices) n_var = nvars(ram_matrices) C_indices = CartesianIndices((n_var, n_var)) diff --git a/src/frontend/common.jl b/src/frontend/common.jl index e89a6cf8b..2734e8f2b 100644 --- a/src/frontend/common.jl +++ b/src/frontend/common.jl @@ -1,7 +1,7 @@ # API methods supported by multiple SEM.jl types """ - params(semobj) -> Vector{Symbol} + params(partable::ParameterTable) -> Vector{Symbol} Return the vector of SEM model parameter identifiers. """ @@ -42,7 +42,7 @@ nlatent_vars(semobj) = length(latent_vars(semobj)) """ param_indices(semobj) -Returns a dict of parameter names and their indices in `semobj`. +Returns a dict of parameter labels and their indices in `semobj`. # Examples ```julia diff --git a/src/frontend/fit/fitmeasures/fit_measures.jl b/src/frontend/fit/fitmeasures/fit_measures.jl index 2fc4dfba0..afdde173b 100644 --- a/src/frontend/fit/fitmeasures/fit_measures.jl +++ b/src/frontend/fit/fitmeasures/fit_measures.jl @@ -14,6 +14,6 @@ end """ fit_measures(sem_fit, args...) -Return a default set of fit measures or the fit measures passed as `arg...`. +Return a default set of fit measures or the fit measures passed as `args...`. """ function fit_measures end diff --git a/src/frontend/fit/fitmeasures/minus2ll.jl b/src/frontend/fit/fitmeasures/minus2ll.jl index 2cb87d79c..ab4d24e53 100644 --- a/src/frontend/fit/fitmeasures/minus2ll.jl +++ b/src/frontend/fit/fitmeasures/minus2ll.jl @@ -41,7 +41,7 @@ end function minus2ll(observed::SemObservedMissing) # fit EM-based mean and cov if not yet fitted # FIXME EM could be very computationally expensive - observed.em_model.fitted || em_mvn(observed) + observed.em_model.fitted || em_mvn!(observed) Σ = observed.em_model.Σ μ = observed.em_model.μ diff --git a/src/frontend/fit/summary.jl b/src/frontend/fit/summary.jl index 8ee134a9c..3071d5653 100644 --- a/src/frontend/fit/summary.jl +++ b/src/frontend/fit/summary.jl @@ -15,7 +15,7 @@ function details(sem_fit::SemFit; show_fitmeasures = false, color = :light_cyan, println("Number of data samples: $(nsamples(sem_fit))") print("\n") printstyled( - "----------------------------------- Model ----------------------------------- \n"; + "----------------------------------- Model ------------------------------------ \n"; color = color, ) print("\n") @@ -23,7 +23,7 @@ function details(sem_fit::SemFit; show_fitmeasures = false, color = :light_cyan, print("\n") if show_fitmeasures printstyled( - "--------------------------------- Fitmeasures --------------------------------- \n"; + "-------------------------------- Fitmeasures --------------------------------- \n"; color = color, ) print("\n") @@ -51,7 +51,7 @@ function details( if show_variables print("\n") printstyled( - "--------------------------------- Variables --------------------------------- \n"; + "---------------------------------- Variables --------------------------------- \n"; color = color, ) print("\n") @@ -242,9 +242,6 @@ function details( ) print("\n") end - - #printstyled("""No need to copy and paste results, you can use CSV.write(DataFrame(my_partable), "myfile.csv")"""; hidden = true) - end function details( @@ -297,9 +294,6 @@ function details( show_columns = show_columns, ) end - - # printstyled("""No need to copy and paste results, you can use CSV.write(DataFrame(my_partable), "myfile.csv")"""; hidden = true) - end function check_round(vec; digits) diff --git a/src/frontend/pretty_printing.jl b/src/frontend/pretty_printing.jl index c1cd72c2f..2fa970f24 100644 --- a/src/frontend/pretty_printing.jl +++ b/src/frontend/pretty_printing.jl @@ -1,3 +1,7 @@ +############################################################## +# Some helpers to implement show methods for SEM.jl objects +############################################################## + function print_field_types(io::IO, struct_instance) fields = fieldnames(typeof(struct_instance)) types = [typeof(getproperty(struct_instance, field)) for field in fields] @@ -25,7 +29,7 @@ function print_type(io::IO, struct_instance) end ############################################################## -# Loss Functions, Implied, +# Loss Function, Implied, Observed, Optimizer ############################################################## function Base.show(io::IO, struct_inst::SemLossFunction) diff --git a/src/frontend/specification/RAMMatrices.jl b/src/frontend/specification/RAMMatrices.jl index 75175a87d..d430e9c01 100644 --- a/src/frontend/specification/RAMMatrices.jl +++ b/src/frontend/specification/RAMMatrices.jl @@ -9,7 +9,7 @@ struct RAMMatrices <: SemSpecification F::SparseMatrixCSC{Float64} M::Union{ParamsVector{Float64}, Nothing} param_labels::Vector{Symbol} - vars::Union{Vector{Symbol}, Nothing} # better call it "variables": it's a mixture of observed and latent (and it gets confusing with get_vars()) + vars::Union{Vector{Symbol}, Nothing} end nparams(ram::RAMMatrices) = nparams(ram.A) diff --git a/src/frontend/specification/Sem.jl b/src/frontend/specification/Sem.jl index 7ba8f7fb7..53858abd4 100644 --- a/src/frontend/specification/Sem.jl +++ b/src/frontend/specification/Sem.jl @@ -45,11 +45,6 @@ Returns the [*observed*](@ref SemObserved) part of a model. """ observed(model::AbstractSemSingle) = model.observed -""" - nsamples(model::AbstractSem) -> Int - -Returns the number of samples from the [*observed*](@ref SemObserved) part of a model. -""" nsamples(model::AbstractSemSingle) = nsamples(observed(model)) """ @@ -156,11 +151,6 @@ end # pretty printing ############################################################## -#= function Base.show(io::IO, sem::Sem{O, I, L, D}) where {O, I, L, D} - lossfuntypes = @. nameof(typeof(sem.loss.functions)) - print(io, "Sem{$(nameof(O)), $(nameof(I)), $lossfuntypes, $(nameof(D))}") -end =# - function Base.show(io::IO, sem::Sem{O, I, L}) where {O, I, L} lossfuntypes = @. string(nameof(typeof(sem.loss.functions))) lossfuntypes = " " .* lossfuntypes .* ("\n") diff --git a/src/frontend/specification/StenoGraphs.jl b/src/frontend/specification/StenoGraphs.jl index 314abcc35..79e17f719 100644 --- a/src/frontend/specification/StenoGraphs.jl +++ b/src/frontend/specification/StenoGraphs.jl @@ -55,7 +55,6 @@ function ParameterTable( value_fixed = columns[:value_fixed] start = columns[:start] param_refs = columns[:label] - # group = Vector{Symbol}(undef, n) for (i, element) in enumerate(graph) edge = element isa ModifiedEdge ? element.edge : element diff --git a/src/frontend/specification/documentation.jl b/src/frontend/specification/documentation.jl index 7697a63c3..e46620fbc 100644 --- a/src/frontend/specification/documentation.jl +++ b/src/frontend/specification/documentation.jl @@ -1,10 +1,3 @@ -""" - param_labels(semobj) -> Vector{Symbol} - -Return the vector of parameter labels (in the same order as [`params`](@ref)). -""" -param_labels(spec::SemSpecification) = spec.param_labels - """ vars(semobj) -> Vector{Symbol} @@ -37,14 +30,22 @@ function latent_vars end latent_vars(spec::SemSpecification) = error("latent_vars(spec::$(typeof(spec))) is not implemented") +""" + param_labels(semobj) -> Vector{Symbol} + +Return the vector of parameter labels (in the same order as [`params`](@ref)). +""" +param_labels(spec::SemSpecification) = spec.param_labels + + """ `ParameterTable`s contain the specification of a structural equation model. # Constructor - (1) ParameterTable(;graph, observed_vars, latent_vars, ...) + (1) ParameterTable(graph; observed_vars, latent_vars, ...) - (2) ParameterTable(ram_matrices) + (2) ParameterTable(ram_matrices; ...) Return a `ParameterTable` constructed from (1) a graph or (2) RAM matrices. @@ -55,7 +56,7 @@ Return a `ParameterTable` constructed from (1) a graph or (2) RAM matrices. - `ram_matrices::RAMMatrices`: a `RAMMatrices` object # Examples -See the online documentation on [Model specification](@ref) and the [ParameterTable interface](@ref). +See the online documentation on [Model specification](@ref) and the [Graph interface](@ref). # Extended help ## Additional keyword arguments @@ -68,7 +69,7 @@ function ParameterTable end # Constructor - (1) EnsembleParameterTable(;graph, observed_vars, latent_vars, groups) + (1) EnsembleParameterTable(graph; observed_vars, latent_vars, groups) (2) EnsembleParameterTable(ps::Pair...; param_labels = nothing) @@ -91,9 +92,9 @@ function EnsembleParameterTable end # Constructor - (1) RAMMatrices(partable::ParameterTable) + (1) RAMMatrices(partable::ParameterTable; param_labels = nothing) - (2) RAMMatrices(;A, S, F, M = nothing, param_labels, vars) + (2) RAMMatrices(;A, S, F, M = nothing, param_labels, vars = nothing) (3) RAMMatrices(partable::EnsembleParameterTable) diff --git a/src/implied/RAM/generic.jl b/src/implied/RAM/generic.jl index 301c455e9..fd2ef61b5 100644 --- a/src/implied/RAM/generic.jl +++ b/src/implied/RAM/generic.jl @@ -6,8 +6,7 @@ Model implied covariance and means via RAM notation. # Constructor - RAM(; - specification, + RAM(;specification, meanstructure = false, gradient = true, kwargs...) @@ -19,9 +18,6 @@ Model implied covariance and means via RAM notation. # Extended help -## Implementation -Subtype of `SemImplied`. - ## RAM notation The model implied covariance matrix is computed as @@ -37,33 +33,32 @@ and for models with a meanstructure, the model implied means are computed as - `param_labels(::RAM) `-> vector of parameter labels - `nparams(::RAM)` -> number of parameters -- `Σ(::RAM)` -> model implied covariance matrix -- `μ(::RAM)` -> model implied mean vector +- `ram.Σ` -> model implied covariance matrix +- `ram.μ` -> model implied mean vector RAM matrices for the current parameter values: -- `A(::RAM)` -- `S(::RAM)` -- `F(::RAM)` -- `M(::RAM)` +- `ram.A` +- `ram.S` +- `ram.F` +- `ram.M` Jacobians of RAM matrices w.r.t to the parameter vector `θ` -- `∇A(::RAM)` -> ``∂vec(A)/∂θᵀ`` -- `∇S(::RAM)` -> ``∂vec(S)/∂θᵀ`` -- `∇M(::RAM)` = ``∂M/∂θᵀ`` +- `ram.∇A` -> ``∂vec(A)/∂θᵀ`` +- `ram.∇S` -> ``∂vec(S)/∂θᵀ`` +- `ram.∇M` = ``∂M/∂θᵀ`` Vector of indices of each parameter in the respective RAM matrix: -- `A_indices(::RAM)` -- `S_indices(::RAM)` -- `M_indices(::RAM)` +- `ram.A_indices` +- `ram.S_indices` +- `ram.M_indices` Additional interfaces -- `F⨉I_A⁻¹(::RAM)` -> ``F(I-A)^{-1}`` -- `F⨉I_A⁻¹S(::RAM)` -> ``F(I-A)^{-1}S`` -- `I_A(::RAM)` -> ``I-A`` -- `has_meanstructure(::RAM)` -> `Val{Bool}` does the model have a meanstructure? +- `ram.F⨉I_A⁻¹` -> ``F(I-A)^{-1}`` +- `ram.F⨉I_A⁻¹S` -> ``F(I-A)^{-1}S`` +- `ram.I_A` -> ``I-A`` Only available in gradient! calls: -- `I_A⁻¹(::RAM)` -> ``(I-A)^{-1}`` +- `ram.I_A⁻¹` -> ``(I-A)^{-1}`` """ mutable struct RAM{MS, A1, A2, A3, A4, A5, A6, V2, M1, M2, M3, M4, S1, S2, S3} <: SemImplied meanstruct::MS @@ -97,7 +92,6 @@ end function RAM(; specification::SemSpecification, - #vech = false, gradient_required = true, meanstructure = false, kwargs..., diff --git a/src/implied/RAM/symbolic.jl b/src/implied/RAM/symbolic.jl index eff193c17..9634bfa89 100644 --- a/src/implied/RAM/symbolic.jl +++ b/src/implied/RAM/symbolic.jl @@ -6,7 +6,8 @@ Subtype of `SemImplied` that implements the RAM notation with symbolic precomput # Constructor - RAMSymbolic(;specification, + RAMSymbolic(; + specification, vech = false, gradient = true, hessian = false, @@ -25,32 +26,25 @@ Subtype of `SemImplied` that implements the RAM notation with symbolic precomput # Extended help -## Implementation -Subtype of `SemImplied`. - ## Interfaces - `param_labels(::RAMSymbolic) `-> vector of parameter ids - `nparams(::RAMSymbolic)` -> number of parameters -- `Σ(::RAMSymbolic)` -> model implied covariance matrix -- `μ(::RAMSymbolic)` -> model implied mean vector +- `ram.Σ` -> model implied covariance matrix +- `ram.μ` -> model implied mean vector Jacobians (only available in gradient! calls) -- `∇Σ(::RAMSymbolic)` -> ``∂vec(Σ)/∂θᵀ`` -- `∇μ(::RAMSymbolic)` -> ``∂μ/∂θᵀ`` +- `ram.∇Σ` -> ``∂vec(Σ)/∂θᵀ`` +- `ram.∇μ` -> ``∂μ/∂θᵀ`` -- `∇Σ_function(::RAMSymbolic)` -> function to overwrite `∇Σ` in place, - i.e. `∇Σ_function(∇Σ, θ)`. Normally, you do not want to use this but simply - query `∇Σ(::RAMSymbolic)`. +- `ram.∇Σ_function` -> function to overwrite `∇Σ` in place, + i.e. `∇Σ_function(∇Σ, θ)`. Typically, you do not want to use this but simply + query `ram.∇Σ`. Hessians -The computation of hessians is more involved, and uses the "chain rule for -hessian matrices". -Therefore, we desribe it at length in the mathematical appendix of the online documentation, -and the relevant interfaces are omitted here. - -Additional interfaces -- `has_meanstructure(::RAMSymbolic)` -> `Val{Bool}` does the model have a meanstructure? +The computation of hessians is more involved. +Therefore, we desribe it in the online documentation, +and the respective interfaces are omitted here. ## RAM notation The model implied covariance matrix is computed as diff --git a/src/implied/abstract.jl b/src/implied/abstract.jl index af51440c6..d4868d746 100644 --- a/src/implied/abstract.jl +++ b/src/implied/abstract.jl @@ -1,4 +1,3 @@ - # vars and params API methods for SemImplied vars(implied::SemImplied) = vars(implied.ram_matrices) observed_vars(implied::SemImplied) = observed_vars(implied.ram_matrices) diff --git a/src/implied/empty.jl b/src/implied/empty.jl index 3b0292e73..82a6c9469 100644 --- a/src/implied/empty.jl +++ b/src/implied/empty.jl @@ -19,11 +19,8 @@ model per group and an additional model with `ImpliedEmpty` and `SemRidge` for t # Extended help ## Interfaces -- `param_labels(::RAMSymbolic) `-> Vector of parameter labels -- `nparams(::RAMSymbolic)` -> Number of parameters - -## Implementation -Subtype of `SemImplied`. +- `param_labels(::ImpliedEmpty) `-> Vector of parameter labels +- `nparams(::ImpliedEmpty)` -> Number of parameters """ struct ImpliedEmpty{A, B, C} <: SemImplied hessianeval::A diff --git a/src/loss/ML/FIML.jl b/src/loss/ML/FIML.jl index ca23ded97..ea8da6b37 100644 --- a/src/loss/ML/FIML.jl +++ b/src/loss/ML/FIML.jl @@ -19,10 +19,6 @@ my_fiml = SemFIML(observed = my_observed, specification = my_parameter_table) # Interfaces Analytic gradients are available. - -# Extended help -## Implementation -Subtype of `SemLossFunction`. """ mutable struct SemFIML{INV, C, L, O, M, IM, I, T, W} <: SemLossFunction hessianeval::ExactHessian diff --git a/src/loss/ML/ML.jl b/src/loss/ML/ML.jl index d14af648c..ec5eb997c 100644 --- a/src/loss/ML/ML.jl +++ b/src/loss/ML/ML.jl @@ -21,11 +21,8 @@ my_ml = SemML(observed = my_observed) ``` # Interfaces -Analytic gradients are available, and for models without a meanstructure, also analytic hessians. - -# Extended help -## Implementation -Subtype of `SemLossFunction`. +Analytic gradients are available, and for models without a meanstructure +and RAMSymbolic implied type, also analytic hessians. """ struct SemML{HE <: HessianEval, INV, M, M2} <: SemLossFunction hessianeval::HE diff --git a/src/loss/WLS/WLS.jl b/src/loss/WLS/WLS.jl index 0fe2c9b3c..dd5be4874 100644 --- a/src/loss/WLS/WLS.jl +++ b/src/loss/WLS/WLS.jl @@ -5,6 +5,7 @@ ############################################################################################ """ Weighted least squares estimation. +At the moment only available with the `RAMSymbolic` implied type. # Constructor @@ -32,11 +33,7 @@ my_wls = SemWLS(observed = my_observed) ``` # Interfaces -Analytic gradients are available, and for models without a meanstructure, also analytic hessians. - -# Extended help -## Implementation -Subtype of `SemLossFunction`. +Analytic gradients are available, and for models without a meanstructure also analytic hessians. """ struct SemWLS{HE <: HessianEval, Vt, St, C} <: SemLossFunction hessianeval::HE @@ -124,7 +121,7 @@ function evaluate!( if issparse(∇σ) gradient .= (σ₋' * V * ∇σ)' else # save one allocation - mul!(gradient, σ₋' * V, ∇σ) # actually transposed, but should be fine for vectors + mul!(gradient, σ₋' * V, ∇σ) end gradient .*= -2 end diff --git a/src/loss/constant/constant.jl b/src/loss/constant/constant.jl index cb5157346..3aed5e27c 100644 --- a/src/loss/constant/constant.jl +++ b/src/loss/constant/constant.jl @@ -20,10 +20,6 @@ Constant loss term. Can be used for comparability to other packages. # Interfaces Analytic gradients and hessians are available. - -# Extended help -## Implementation -Subtype of `SemLossFunction`. """ struct SemConstant{C} <: SemLossFunction hessianeval::ExactHessian diff --git a/src/loss/regularization/ridge.jl b/src/loss/regularization/ridge.jl index aee521624..90cbcc231 100644 --- a/src/loss/regularization/ridge.jl +++ b/src/loss/regularization/ridge.jl @@ -24,10 +24,6 @@ my_ridge = SemRidge(;α_ridge = 0.02, which_ridge = [:λ₁, :λ₂, :ω₂₃], # Interfaces Analytic gradients and hessians are available. - -# Extended help -## Implementation -Subtype of `SemLossFunction`. """ struct SemRidge{P, W1, W2, GT, HT} <: SemLossFunction hessianeval::ExactHessian diff --git a/src/objective_gradient_hessian.jl b/src/objective_gradient_hessian.jl index 4aafe4235..06f39329f 100644 --- a/src/objective_gradient_hessian.jl +++ b/src/objective_gradient_hessian.jl @@ -95,7 +95,8 @@ function evaluate!(objective, gradient, hessian, model::AbstractSemSingle, param end ############################################################################################ -# methods for SemFiniteDiff (approximate gradient and hessian with finite differences of objective) +# methods for SemFiniteDiff +# (approximate gradient and hessian with finite differences of objective) ############################################################################################ function evaluate!(objective, gradient, hessian, model::SemFiniteDiff, params) @@ -157,14 +158,6 @@ function evaluate!(objective, gradient, hessian, ensemble::SemEnsemble, params) return objective end -# throw an error by default if gradient! and hessian! are not implemented - -#= gradient!(lossfun::SemLossFunction, par, model) = - throw(ArgumentError("gradient for $(typeof(lossfun).name.wrapper) is not available")) - -hessian!(lossfun::SemLossFunction, par, model) = - throw(ArgumentError("hessian for $(typeof(lossfun).name.wrapper) is not available")) =# - ############################################################################################ # Documentation ############################################################################################ diff --git a/src/observed/EM.jl b/src/observed/EM.jl index beac45ca8..46d0622be 100644 --- a/src/observed/EM.jl +++ b/src/observed/EM.jl @@ -3,16 +3,32 @@ ############################################################################################ # An EM Algorithm for MVN-distributed Data with missing values -# Adapted from supplementary Material to the book Machine Learning: A Probabilistic Perspective -# Copyright (2010) Kevin Murphy and Matt Dunham -# found at https://github.com/probml/pmtk3/blob/master/toolbox/BasicModels/gauss/sub/gaussMissingFitEm.m -# and at https://github.com/probml/pmtk3/blob/master/toolbox/Algorithms/optimization/emAlgo.m +# Adapted from https://github.com/probml/pmtk3, licensed as +#= The MIT License -# what about random restarts? +Copyright (2010) Kevin Murphy and Matt Dunham + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. =# # outer function --------------------------------------------------------------------------- """ - em_mvn(; + em_mvn!(; observed::SemObservedMissing, start_em = start_em_observed, max_iter_em = 100, @@ -22,7 +38,7 @@ Estimates the covariance matrix and mean vector of the normal distribution via expectation maximization for `observed`. Overwrites the statistics stored in `observed`. """ -function em_mvn( +function em_mvn!( observed::SemObservedMissing; start_em = start_em_observed, max_iter_em = 100, @@ -32,7 +48,7 @@ function em_mvn( nvars = nobserved_vars(observed) nsamps = nsamples(observed) - # preallocate stuff? + # preallocate stuff 𝔼x_pre = zeros(nvars) 𝔼xxᵀ_pre = zeros(nvars, nvars) @@ -45,9 +61,6 @@ function em_mvn( end end - # ess = 𝔼x, 𝔼xxᵀ, ismissing, missingRows, nsamps - # estepFn = (em_model, data) -> estep(em_model, data, EXsum, EXXsum, ismissing, missingRows, nsamps) - # initialize em_model = start_em(observed; kwargs...) em_model_prev = EmMVNModel(zeros(nvars, nvars), zeros(nvars), false) @@ -65,13 +78,11 @@ function em_mvn( @warn "EM Algorithm for MVN missing data did not converge. Likelihood for FIML is not interpretable. Maybe try passing different starting values via 'start_em = ...' " elseif iter > 1 - # done = isapprox(ll, ll_prev; rtol = rtol) done = isapprox(em_model_prev.μ, em_model.μ; rtol = rtol_em) & isapprox(em_model_prev.Σ, em_model.Σ; rtol = rtol_em) end - # print("$iter \n") iter = iter + 1 em_model_prev.μ, em_model_prev.Σ = em_model.μ, em_model.Σ end @@ -135,23 +146,7 @@ end function em_mvn_Mstep!(em_model, nsamples, 𝔼x, 𝔼xxᵀ) em_model.μ = 𝔼x / nsamples Σ = Symmetric(𝔼xxᵀ / nsamples - em_model.μ * em_model.μ') - - # ridge Σ - # while !isposdef(Σ) - # Σ += 0.5I - # end - em_model.Σ = Σ - - # diagonalization - #if !isposdef(Σ) - # print("Matrix not positive definite") - # em_model.Σ .= 0 - # em_model.Σ[diagind(em_model.Σ)] .= diag(Σ) - #else - # em_model.Σ = Σ - #end - return nothing end @@ -178,7 +173,6 @@ function start_em_simple(observed::SemObservedMissing; kwargs...) μ = zeros(nvars) Σ = rand(nvars, nvars) Σ = Σ * Σ' - # Σ = Matrix(1.0I, nvars, nvars) return EmMVNModel(Σ, μ, false) end diff --git a/src/observed/abstract.jl b/src/observed/abstract.jl index bb92ea12e..cf5000e4f 100644 --- a/src/observed/abstract.jl +++ b/src/observed/abstract.jl @@ -48,7 +48,7 @@ end # function to prepare input data shared by SemObserved implementations # returns tuple of -# 1) the matrix of data +# 1) the data matrix # 2) the observed variable symbols that match matrix columns # 3) the permutation of the original observed_vars (nothing if no reordering) # If observed_vars is not specified, the vars order is taken from the specification. diff --git a/src/observed/covariance.jl b/src/observed/covariance.jl index 221ef5ca3..f81fe8e57 100644 --- a/src/observed/covariance.jl +++ b/src/observed/covariance.jl @@ -9,10 +9,9 @@ const SemObservedCovariance{S} = SemObservedData{Nothing, S} SemObservedCovariance(; specification, obs_cov, - obs_colnames = nothing, - meanstructure = false, + nsamples, obs_mean = nothing, - nsamples::Integer, + observed_vars = nothing, kwargs...) Construct [`SemObserved`](@ref) without providing the observations data, @@ -21,12 +20,13 @@ but with the covariations (`obs_cov`) and the means (`obs_means`) of the observe Returns [`SemObservedCovariance`](@ref) object. # Arguments -- `obs_cov`: pre-computed covariations of the observed variables +- `obs_cov`: pre-computed covariance matrix of the observed variables +- `nsamples::Integer`: number of samples (observed data points) used to compute `obs_cov` and `obs_means`, + used for calculating fit statistics - `obs_mean`: optional pre-computed means of the observed variables - `observed_vars::AbstractVector`: IDs of the observed variables (rows and columns of the `obs_cov` matrix) - `specification`: optional SEM specification ([`SemSpecification`](@ref)) -- `nsamples::Number`: number of samples (observed data points) used to compute `obs_cov` and `obs_means` - necessary for calculating fit statistics + """ function SemObservedCovariance(; obs_cov::AbstractMatrix, diff --git a/src/observed/data.jl b/src/observed/data.jl index b6ddaa43d..7ba38edf5 100644 --- a/src/observed/data.jl +++ b/src/observed/data.jl @@ -10,9 +10,9 @@ For observed data without missings. kwargs...) # Arguments -- `specification`: optional SEM specification ([`SemSpecification`](@ref)) - `data`: observed data -- *DataFrame* or *Matrix* - `observed_vars::Vector{Symbol}`: column names of the data (if the object passed as data does not have column names, i.e. is not a data frame) +- `specification`: optional SEM specification ([`SemSpecification`](@ref)) # Extended help ## Interfaces @@ -20,11 +20,8 @@ For observed data without missings. - `nobserved_vars(::SemObservedData)` -> number of observed (manifested) variables - `samples(::SemObservedData)` -> observed data -- `obs_cov(::SemObservedData)` -> observed.obs_cov -- `obs_mean(::SemObservedData)` -> observed.obs_mean - -## Implementation -Subtype of `SemObserved` +- `obs_cov(::SemObservedData)` -> observed covariance matrix +- `obs_mean(::SemObservedData)` -> observed mean vector """ struct SemObservedData{D <: Union{Nothing, AbstractMatrix}, S <: Number} <: SemObserved data::D @@ -48,10 +45,6 @@ function SemObservedData(; return SemObservedData(data, obs_vars, obs_cov, vec(obs_mean), size(data, 1)) end -############################################################################################ -### Recommended methods -############################################################################################ - ############################################################################################ ### additional methods ############################################################################################ diff --git a/src/observed/missing.jl b/src/observed/missing.jl index cf699252e..ac8b7ea5f 100644 --- a/src/observed/missing.jl +++ b/src/observed/missing.jl @@ -25,9 +25,9 @@ For observed data with missing values. kwargs...) # Arguments -- `specification`: optional SEM model specification ([`SemSpecification`](@ref)) - `data`: observed data - `observed_vars::Vector{Symbol}`: column names of the data (if the object passed as data does not have column names, i.e. is not a data frame) +- `specification`: optional SEM model specification ([`SemSpecification`](@ref)) # Extended help ## Interfaces @@ -35,10 +35,14 @@ For observed data with missing values. - `nobserved_vars(::SemObservedMissing)` -> number of observed variables - `samples(::SemObservedMissing)` -> data matrix (contains both measured and missing values) -- `em_model(::SemObservedMissing)` -> `EmMVNModel` that contains the covariance matrix and mean vector found via expectation maximization -## Implementation -Subtype of `SemObserved` +## Expectation maximization +`em_mvn!(::SemObservedMissing)` can be called to fit a covariance matrix and mean vector to the data +using an expectation maximization (EM) algorithm under the assumption of multivariate normality. +After, the following methods are available: +- `em_model(::SemObservedMissing)` -> `EmMVNModel` that contains the covariance matrix and mean vector found via EM +- `obs_cov(::SemObservedData)` -> EM covariance matrix +- `obs_mean(::SemObservedData)` -> EM mean vector """ struct SemObservedMissing{T <: Real, S <: Real, E <: EmMVNModel} <: SemObserved data::Matrix{Union{T, Missing}} @@ -91,10 +95,6 @@ function SemObservedMissing(; ) end -############################################################################################ -### Recommended methods -############################################################################################ - ############################################################################################ ### Additional methods ############################################################################################ diff --git a/src/optimizer/Empty.jl b/src/optimizer/Empty.jl index 45a20db55..1bf0c30ac 100644 --- a/src/optimizer/Empty.jl +++ b/src/optimizer/Empty.jl @@ -8,12 +8,6 @@ an optimizer part. # Constructor SemOptimizerEmpty() - -# Extended help - -## Implementation - -Subtype of `SemOptimizer`. """ struct SemOptimizerEmpty <: SemOptimizer{:Empty} end diff --git a/src/optimizer/abstract.jl b/src/optimizer/abstract.jl index 2487b7c52..3b1e98842 100644 --- a/src/optimizer/abstract.jl +++ b/src/optimizer/abstract.jl @@ -20,11 +20,17 @@ the online documentation on [Starting values](@ref). # Examples ```julia -fit( - my_model; +fit(my_model; start_val = start_simple, start_covariances_latent = 0.5) ``` + +```julia +using Optim + +fit(my_model; + algorithm = BFGS()) +``` """ function fit(optim::SemOptimizer, model::AbstractSem; start_val = nothing, kwargs...) start_params = prepare_start_params(start_val, model; kwargs...) @@ -49,6 +55,10 @@ prepare_start_params(start_val::Nothing, model::AbstractSemSingle; kwargs...) = prepare_start_params(start_val::Nothing, model::AbstractSem; kwargs...) = start_simple(model; kwargs...) +# first argument is a function +prepare_start_params(start_val, model::AbstractSem; kwargs...) = + start_val(model; kwargs...) + function prepare_start_params(start_val::AbstractVector, model::AbstractSem; kwargs...) (length(start_val) == nparams(model)) || throw( DimensionMismatch( diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl index 8f5404bc2..e01a606d4 100644 --- a/src/optimizer/optim.jl +++ b/src/optimizer/optim.jl @@ -16,7 +16,7 @@ Connects to `Optim.jl` as the optimization backend. kwargs...) # Arguments -- `algorithm`: optimization algorithm. +- `algorithm`: optimization algorithm from `Optim.jl` - `options::Optim.Options`: options for the optimization algorithm # Usage @@ -110,7 +110,7 @@ function fit( upper_bounds::Union{AbstractVector, AbstractDict, Nothing} = nothing, lower_bound = -Inf, upper_bound = Inf, - variance_lower_bound::Number = 0.0, + variance_lower_bound::Number = -Inf, variance_upper_bound::Number = Inf, kwargs..., ) diff --git a/src/package_extensions/SEMNLOptExt.jl b/src/package_extensions/SEMNLOptExt.jl index 69721ac94..64c4cff04 100644 --- a/src/package_extensions/SEMNLOptExt.jl +++ b/src/package_extensions/SEMNLOptExt.jl @@ -49,10 +49,6 @@ see [Constrained optimization](@ref) in our online documentation. - `local_options(::SemOptimizerNLopt)` - `equality_constraints(::SemOptimizerNLopt)` - `inequality_constraints(::SemOptimizerNLopt)` - -## Implementation - -Subtype of `SemOptimizer`. """ struct SemOptimizerNLopt{A, A2, B, B2, C} <: SemOptimizer{:NLopt} algorithm::A diff --git a/src/package_extensions/SEMProximalOptExt.jl b/src/package_extensions/SEMProximalOptExt.jl index 5d4007504..ad4c2da2a 100644 --- a/src/package_extensions/SEMProximalOptExt.jl +++ b/src/package_extensions/SEMProximalOptExt.jl @@ -1,5 +1,6 @@ """ Connects to `ProximalAlgorithms.jl` as the optimization backend. +Can be used for regularized SEM, for a tutorial see the online docs on [Regularization](@ref). # Constructor @@ -11,8 +12,13 @@ Connects to `ProximalAlgorithms.jl` as the optimization backend. # Arguments - `algorithm`: optimization algorithm. -- `operator_g`: gradient of the objective function -- `operator_h`: optional hessian of the objective function +- `operator_g`: proximal operator (e.g., regularization penalty) +- `operator_h`: optional second proximal operator + +# Usage +All algorithms and operators from `ProximalAlgorithms.jl` are available, +for more information see the online docs on [Regularization](@ref) and +the documentation of `ProximalAlgorithms.jl` / `ProximalOperators.jl`. """ mutable struct SemOptimizerProximal{A, B, C} <: SemOptimizer{:Proximal} algorithm::A diff --git a/src/types.jl b/src/types.jl index 64a4acbac..44d472ebf 100644 --- a/src/types.jl +++ b/src/types.jl @@ -106,7 +106,7 @@ abstract type SemObserved end """ Supertype of all objects that can serve as the implied field of a SEM. -Computed model-implied values that should be compared with the observed data to find parameter estimates, +Computes model-implied values that should be compared with the observed data to find parameter estimates, e. g. the model implied covariance or mean. If you would like to implement a different notation, e.g. LISREL, you should implement a subtype of SemImplied. """ @@ -168,7 +168,7 @@ end # ensemble models ############################################################################################ """ - (1) SemEnsemble(models..., weights = nothing, kwargs...) + (1) SemEnsemble(models...; weights = nothing, kwargs...) (2) SemEnsemble(;specification, data, groups, column = :group, kwargs...) @@ -189,6 +189,8 @@ Returns a SemEnsemble with fields - `sems::Tuple`: `AbstractSem`s. - `weights::Vector`: Weights for each model. - `param_labels::Vector`: Stores parameter labels and their position. + +For instructions on multigroup models, see the online documentation. """ struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, I} <: AbstractSemCollection n::N diff --git a/test/unit_tests/bootstrap.jl b/test/unit_tests/bootstrap.jl deleted file mode 100644 index a2d5b6832..000000000 --- a/test/unit_tests/bootstrap.jl +++ /dev/null @@ -1,5 +0,0 @@ -solution_ml = fit(model_ml) -bs = se_bootstrap(solution_ml; n_boot = 20) - -update_se_hessian!(partable, solution_ml) -update_partable!(partable, solution_ml, bs, :se_boot) diff --git a/test/unit_tests/partable_dataframe.jl b/test/unit_tests/partable_dataframe.jl deleted file mode 100644 index 7786e2ec1..000000000 --- a/test/unit_tests/partable_dataframe.jl +++ /dev/null @@ -1,3 +0,0 @@ -# Convert PramaterTable to DataFrame - -# Convert EnsembleParameterTable to DataFrame diff --git a/test/unit_tests/ridge_id.jl b/test/unit_tests/ridge_id.jl deleted file mode 100644 index 13366bf01..000000000 --- a/test/unit_tests/ridge_id.jl +++ /dev/null @@ -1,8 +0,0 @@ -model_ridge_id = Sem( - specification = spec, - data = dat, - loss = (SemML, SemRidge), - α_ridge = 0.001, - which_ridge = [:x16, :x17, :x18, :x19, :x20], - optimizer = semoptimizer, -) diff --git a/test/unit_tests/sorting.jl b/test/unit_tests/sorting.jl deleted file mode 100644 index 3c61e13c4..000000000 --- a/test/unit_tests/sorting.jl +++ /dev/null @@ -1,17 +0,0 @@ -############################################################################ -### test variables sorting -############################################################################ - -sort_vars!(partable) - -model_ml_sorted = Sem(specification = partable, data = dat) - -@testset "graph sorting" begin - @test model_ml_sorted.implied.I_A isa LowerTriangular -end - -@testset "ml_solution_sorted" begin - solution_ml_sorted = fit(model_ml_sorted) - update_estimate!(partable, solution_ml_sorted) - @test test_estimates(par_ml, partable, 0.01) -end diff --git a/test/unit_tests/start_val.jl b/test/unit_tests/start_val.jl deleted file mode 100644 index 8b1378917..000000000 --- a/test/unit_tests/start_val.jl +++ /dev/null @@ -1 +0,0 @@ - diff --git a/test/unit_tests/unit_tests.jl b/test/unit_tests/unit_tests.jl index 883098707..7189addd4 100644 --- a/test/unit_tests/unit_tests.jl +++ b/test/unit_tests/unit_tests.jl @@ -33,11 +33,3 @@ selected_tests = isempty(ARGS) ? collect(keys(available_tests)) : ARGS end end end - -@safetestset "SemSpecification" begin - include("specification.jl") -end - -@safetestset "Sem model" begin - include("model.jl") -end From e9351dde572da1646b2cff23c1a785a050a5f9ce Mon Sep 17 00:00:00 2001 From: Maximilian Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> Date: Wed, 16 Apr 2025 14:15:46 +0200 Subject: [PATCH 2/4] Hotfix/optim (#262) --- docs/src/developer/optimizer.md | 4 ++-- src/optimizer/optim.jl | 4 ++-- .../examples/recover_parameters/recover_parameters_twofact.jl | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/developer/optimizer.md b/docs/src/developer/optimizer.md index 0a85eb6a6..9e01ac87c 100644 --- a/docs/src/developer/optimizer.md +++ b/docs/src/developer/optimizer.md @@ -16,7 +16,7 @@ SemOptimizer{:Name}(args...; kwargs...) = SemOptimizerName(args...; kwargs...) SemOptimizerName(; algorithm = LBFGS(), - options = Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8), + options = Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8), kwargs..., ) = SemOptimizerName(algorithm, options) @@ -40,7 +40,7 @@ Similarly, `SemOptimizer{:Name}(args...; kwargs...) = SemOptimizerName(args...; ```julia SemOptimizerName(; algorithm = LBFGS(), - options = Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8), + options = Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8), kwargs..., ) = SemOptimizerName(algorithm, options) ``` diff --git a/src/optimizer/optim.jl b/src/optimizer/optim.jl index e01a606d4..bd57942d8 100644 --- a/src/optimizer/optim.jl +++ b/src/optimizer/optim.jl @@ -12,7 +12,7 @@ Connects to `Optim.jl` as the optimization backend. SemOptimizerOptim(; algorithm = LBFGS(), - options = Optim.Options(;f_tol = 1e-10, x_tol = 1.5e-8), + options = Optim.Options(;f_reltol = 1e-10, x_abstol = 1.5e-8), kwargs...) # Arguments @@ -67,7 +67,7 @@ SemOptimizer{:Optim}(args...; kwargs...) = SemOptimizerOptim(args...; kwargs...) SemOptimizerOptim(; algorithm = LBFGS(), - options = Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8), + options = Optim.Options(;f_reltol = 1e-10, x_abstol = 1.5e-8), kwargs..., ) = SemOptimizerOptim(algorithm, options) diff --git a/test/examples/recover_parameters/recover_parameters_twofact.jl b/test/examples/recover_parameters/recover_parameters_twofact.jl index a3e426cbc..a7b4cec9a 100644 --- a/test/examples/recover_parameters/recover_parameters_twofact.jl +++ b/test/examples/recover_parameters/recover_parameters_twofact.jl @@ -70,7 +70,7 @@ objective!(model_ml, true_val) optimizer = SemOptimizerOptim( BFGS(; linesearch = BackTracking(order = 3), alphaguess = InitialHagerZhang()),# m = 100), - Optim.Options(; f_tol = 1e-10, x_tol = 1.5e-8), + Optim.Options(; f_reltol = 1e-10, x_abstol = 1.5e-8), ) solution_ml = fit(optimizer, model_ml) From d0ff25b70ae22945cdc4315679a10efa97be64ad Mon Sep 17 00:00:00 2001 From: Maximilian Ernst Date: Wed, 16 Apr 2025 14:31:51 +0200 Subject: [PATCH 3/4] version number --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ae0fa5ba4..376347083 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "StructuralEquationModels" uuid = "383ca8c5-e4ff-4104-b0a9-f7b279deed53" authors = ["Maximilian Ernst", "Aaron Peikert"] -version = "0.4.1" +version = "0.4.2" [deps] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" From 7358daac1c656e7eb178db344e31e5313135f904 Mon Sep 17 00:00:00 2001 From: Maximilian Ernst <34346372+Maximilian-Stefan-Ernst@users.noreply.github.com> Date: Sat, 19 Apr 2025 11:42:56 +0200 Subject: [PATCH 4/4] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 1c3714546..68b7bcb71 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ | [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://structuralequationmodels.github.io/StructuralEquationModels.jl/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://structuralequationmodels.github.io/StructuralEquationModels.jl/dev/) | [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![Github Action CI](https://github.com/StructuralEquationModels/StructuralEquationModels.jl/workflows/CI_extended/badge.svg)](https://github.com/StructuralEquationModels/StructuralEquationModels.jl/actions/) [![codecov](https://codecov.io/gh/StructuralEquationModels/StructuralEquationModels.jl/branch/main/graph/badge.svg?token=P2kjzpvM4V)](https://codecov.io/gh/StructuralEquationModels/StructuralEquationModels.jl) | [![DOI](https://zenodo.org/badge/228649704.svg)](https://zenodo.org/badge/latestdoi/228649704) | > [!NOTE] -> Check out our [preprint](https://formal-methods-mpi.github.io/pkgmanuscript/manuscript.pdf) on the package! +> Check out our [preprint](https://doi.org/10.31234/osf.io/zwe8g_v1) on the package! # What is this Package for?