Cleanup FIML, enhance em_mvn()

Project.toml

@@ -16,6 +16,7 @@ NLSolversBase = "d41bc354-129a-5804-8e4c-c37616107c6c"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/StructuralEquationModels.jl

@@ -14,7 +14,8 @@ using LinearAlgebra,
     StenoGraphs,
     LazyArtifacts,
     DelimitedFiles,
-    DataFrames
+    DataFrames,
+    ProgressMeter
 
 import StatsAPI: params, coef, coefnames, dof, fit, nobs, coeftable
 

src/additional_functions/helper.jl

@@ -14,12 +14,6 @@ function neumann_series(mat::SparseMatrixCSC; maxiter::Integer = size(mat, 1))
     return inverse
 end
 
-function batch_inv!(fun, model)
-    for i in 1:size(fun.inverses, 1)
-        fun.inverses[i] .= LinearAlgebra.inv!(fun.choleskys[i])
-    end
-end
-
 # 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)
@@ -75,6 +69,38 @@ function elimination_matrix(n::Integer)
     return L
 end
 
+# truncate eigenvalues of a symmetric matrix and return the result
+function trunc_eigvals(
+    mtx::AbstractMatrix{T},
+    min_eigval::Number;
+    mtx_label::AbstractString = "matrix",
+    verbose::Bool = false,
+) where {T}
+    size(mtx, 1) == size(mtx, 2) ||
+        throw(ArgumentError("Matrix must be square, $(size(mtx)) given"))
+    issymmetric(mtx) || throw(ArgumentError("Matrix must be symmetric"))
+
+    # eigen decomposition of the mtx
+    mtx_eig = eigen(convert(Matrix{T}, mtx))
+    verbose &&
+        @info "min(eigvals($mtx_label))=$(Base.minimum(mtx_eig.values)), N(eigvals < $min_eigval) = $(sum(<(min_eigval), mtx_eig.values))"
+
+    eigmin = Base.minimum(mtx_eig.values)
+    if eigmin < min_eigval
+        # substitute small eigvals with min_eigval
+        eigvals_mtx = Diagonal(max.(mtx_eig.values, min_eigval))
+        newmtx = mtx_eig.vectors * eigvals_mtx * mtx_eig.vectors'
+        StatsBase._symmetrize!(newmtx)
+        if verbose
+            Δmtx = newmtx .- mtx
+            @info "Δ($mtx_label, posdef)=$(norm(Δmtx, 2)), min,max(Δᵢ)=$(extrema(Δmtx))"
+        end
+        return newmtx
+    else
+        return mtx
+    end
+end
+
 # returns the vector of non-unique values in the order of appearance
 # each non-unique values is reported once
 function nonunique(values::AbstractVector)

src/additional_functions/start_val/start_fabin3.jl

@@ -11,19 +11,10 @@ function start_fabin3(model::AbstractSemSingle; kwargs...)
     return start_fabin3(model.observed, model.implied, model.loss.functions..., kwargs...)
 end
 
-function start_fabin3(observed, implied, args...; kwargs...)
+function start_fabin3(observed::SemObserved, implied::SemImplied, args...; kwargs...)
     return start_fabin3(implied.ram_matrices, obs_cov(observed), obs_mean(observed))
 end
 
-# SemObservedMissing
-function start_fabin3(observed::SemObservedMissing, implied, args...; kwargs...)
-    if !observed.em_model.fitted
-        em_mvn!(observed; kwargs...)
-    end
-
-    return start_fabin3(implied.ram_matrices, observed.em_model.Σ, observed.em_model.μ)
-end
-
 function start_fabin3(
     ram_matrices::RAMMatrices,
     Σ::AbstractMatrix,

src/frontend/fit/fitmeasures/minus2ll.jl

@@ -40,22 +40,22 @@ end
 # compute likelihood for missing data - H1 -------------------------------------------------
 # -2ll =  ∑ log(2π)*(nᵢ + mᵢ) + ln(Σᵢ) + (mᵢ - μᵢ)ᵀ Σᵢ⁻¹ (mᵢ - μᵢ)) + tr(SᵢΣᵢ)
 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.Σ
-    μ = observed.em_model.μ
+    Σ, μ = obs_cov(observed), obs_mean(observed)
 
+    # FIXME: this code is duplicate to objective(fiml, ...)
     F = sum(observed.patterns) do pat
         # implied covariance/mean
         Σᵢ = Σ[pat.measured_mask, pat.measured_mask]
         Σᵢ_chol = cholesky!(Σᵢ)
         ld = logdet(Σᵢ_chol)
         Σᵢ⁻¹ = LinearAlgebra.inv!(Σᵢ_chol)
-        meandiffᵢ = pat.measured_mean - μ[pat.measured_mask]
+        μ_diffᵢ = pat.measured_mean - μ[pat.measured_mask]
 
-        F_one_pattern(meandiffᵢ, Σᵢ⁻¹, pat.measured_cov, ld, nsamples(pat))
+        F_pat = ld + dot(μ_diffᵢ, Σᵢ⁻¹, μ_diffᵢ)
+        if nsamples(pat) > 1
+            F_pat += dot(pat.measured_cov, Σᵢ⁻¹)
+        end
+        F_pat * nsamples(pat)
     end
 
     F += log(2π) * sum(pat -> nsamples(pat) * nmeasured_vars(pat), observed.patterns)