diff --git a/Project.toml b/Project.toml
index cfbecee3..6690ecbc 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.5.0"
+version = "0.5.1"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -27,7 +27,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
-julia = "1.10, 1.11, 1.12"
+julia = "1.10, 1.11, 1.12, 1.13"
 StenoGraphs = "0.5"
 DataFrames = "1"
 Distributions = "0.25"
diff --git a/src/additional_functions/simulation.jl b/src/additional_functions/simulation.jl
index e85e9d5c..d5bdbe91 100644
--- a/src/additional_functions/simulation.jl
+++ b/src/additional_functions/simulation.jl
@@ -34,6 +34,9 @@ Distributions.rand(loss::SemLoss, params, n::Integer) = rand(SEM.implied(loss),
 
 Distributions.rand(model::Sem, params, n::Integer) = rand(sem_term(model), params, n)
 
+Distributions.rand(wrapper::SemFiniteDiff, params, n::Integer) =
+    rand(wrapper.model, params, n)
+
 # rand() overloads without SEM params
 Distributions.rand(implied::Union{SemImplied, SemLoss, Sem}, n::Integer) =
     Distributions.rand(implied, nothing, n)
diff --git a/src/frontend/fit/fitmeasures/fit_measures.jl b/src/frontend/fit/fitmeasures/fit_measures.jl
index 7fabc950..c35ec9ca 100644
--- a/src/frontend/fit/fitmeasures/fit_measures.jl
+++ b/src/frontend/fit/fitmeasures/fit_measures.jl
@@ -1,4 +1,4 @@
-const DEFAULT_FIT_MEASURES = [AIC, BIC, dof, χ², p_value, nparams, RMSEA, CFI]
+const DEFAULT_FIT_MEASURES = [minus2ll, AIC, BIC, dof, χ², p_value, nparams, RMSEA, CFI]
 
 fit_measures(fit, measures::AbstractVector) = Dict(Symbol(fn) => fn(fit) for fn in measures)
 fit_measures(fit, measures...) = fit_measures(fit, measures)
diff --git a/src/frontend/fit/standard_errors/hessian.jl b/src/frontend/fit/standard_errors/hessian.jl
index 80b96d33..7d67b59e 100644
--- a/src/frontend/fit/standard_errors/hessian.jl
+++ b/src/frontend/fit/standard_errors/hessian.jl
@@ -34,22 +34,19 @@ function se_hessian(fit::SemFit; method = :finitediff)
     return [sqrt(c * H_inv[i]) for i in diagind(H_inv)]
 end
 
-# Addition functions -------------------------------------------------------------
-H_scaling(loss::SemML) = 2 / (nsamples(loss) - 1)
+# Additional functions -------------------------------------------------------------
 
-function H_scaling(loss::SemWLS)
-    @warn "Standard errors for WLS are only correct if a GLS weight matrix (the default) is used."
-    return 2 / (nsamples(loss) - 1)
-end
-
-H_scaling(loss::SemFIML) = 2 / nsamples(loss)
+H_nsamples(loss::SemML) = nsamples(loss) - 1
+H_nsamples(loss::SemWLS) = nsamples(loss) - 1
+H_nsamples(loss::SemFIML) = nsamples(loss)
+H_nsamples(loss::SemLoss) = nsamples(loss)
+H_nsamples(wrapper::SemLossFiniteDiff) = H_nsamples(_unwrap(wrapper))
 
 function H_scaling(model::AbstractSem)
     semterms = SEM.sem_terms(model)
-    if length(semterms) > 1
-        #@warn "Hessian scaling for multiple loss functions is not implemented yet"
-        return 2 / nsamples(model)
-    else
-        return length(semterms) >= 1 ? H_scaling(loss(semterms[1])) : 1.0
+    isempty(semterms) && return 1.0
+    if any(term -> _unwrap(loss(term)) isa SemWLS, semterms)
+        @warn "Standard errors for WLS are only correct if a GLS weight matrix (the default) is used."
     end
+    return 2 / sum(term -> H_nsamples(loss(term)), semterms)
 end
diff --git a/src/frontend/specification/ParameterTable.jl b/src/frontend/specification/ParameterTable.jl
index c9b9dc24..72158d1b 100644
--- a/src/frontend/specification/ParameterTable.jl
+++ b/src/frontend/specification/ParameterTable.jl
@@ -277,7 +277,7 @@ function update_partable!(
 
     for (i, par) in enumerate(partable.columns[:label])
         if par == :const
-            coldata[i] = !isnothing(default) ? (isvec_def ? default[i] : default) : zero(T)
+            coldata[i] = !isnothing(default) ? (isvec_def ? default[i] : default) : T(NaN)
         elseif haskey(params, par)
             coldata[i] = params[par]
         else
diff --git a/test/examples/helper.jl b/test/examples/helper.jl
index fed95f3c..d461c225 100644
--- a/test/examples/helper.jl
+++ b/test/examples/helper.jl
@@ -59,6 +59,7 @@ fitmeasure_semjl_to_lavaan = Dict(
     :nparams => "npar",
     :RMSEA => "rmsea",
     :CFI => "cfi",
+    :minus2ll => "logl",
 )
 
 function test_fitmeasures(
@@ -81,11 +82,16 @@ function test_fitmeasures(
             @test ismissing(measure)
         else
             measure_lav = measures_lav.x[lav_ix]
+            measure_lav = name == :minus2ll ? -2measure_lav : measure_lav
             @test measure ≈ measure_lav rtol = rtol atol = atol
         end
     end
 end
 
+# LinearAlgebra v1.13 ignores`norm` keyword for isapprox on arrays (issue #1675)
+isapprox_infnorm(x::AbstractArray, y::AbstractArray; atol::Real = 0, rtol::Real = 0) =
+    norm(x - y, Inf) <= max(atol, rtol * max(norm(x, Inf), norm(y, Inf)))
+
 function test_estimates(
     partable::ParameterTable,
     partable_lav;
@@ -103,10 +109,9 @@ function test_estimates(
     @test !any(isnan, expected)
 
     if skip # workaround skip=false not supported in earlier versions
-        @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf) skip =
-            skip
+        @test isapprox_infnorm(actual, expected; atol = atol, rtol = rtol) skip = skip
     else
-        @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf)
+        @test isapprox_infnorm(actual, expected; atol = atol, rtol = rtol)
     end
 end
 
@@ -136,10 +141,9 @@ function test_estimates(
     @test !any(isnan, expected)
 
     if skip # workaround skip=false not supported in earlier versions
-        @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf) skip =
-            skip
+        @test isapprox_infnorm(actual, expected; atol = atol, rtol = rtol) skip = skip
     else
-        @test actual ≈ expected rtol = rtol atol = atol norm = Base.Fix2(norm, Inf)
+        @test isapprox_infnorm(actual, expected; atol = atol, rtol = rtol)
     end
 end
 
diff --git a/test/examples/multigroup/build_models.jl b/test/examples/multigroup/build_models.jl
index 47bdea22..a706e966 100644
--- a/test/examples/multigroup/build_models.jl
+++ b/test/examples/multigroup/build_models.jl
@@ -58,7 +58,7 @@ end
     test_estimates(
         partable,
         solution_lav[:parameter_estimates_ml];
-        atol = 1e-3,
+        atol = 1e-4,
         col = :se,
         lav_col = :se,
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
@@ -118,7 +118,7 @@ end
     test_estimates(
         partable_s,
         solution_lav[:parameter_estimates_ml];
-        atol = 1e-3,
+        atol = 1e-4,
         col = :se,
         lav_col = :se,
         lav_groups = Dict(:Pasteur => 1, :Grant_White => 2),
diff --git a/test/examples/political_democracy/constructor.jl b/test/examples/political_democracy/constructor.jl
index 2efa5abe..1376cfe3 100644
--- a/test/examples/political_democracy/constructor.jl
+++ b/test/examples/political_democracy/constructor.jl
@@ -191,8 +191,13 @@ if opt_engine == :Optim
     end
 
     @testset "ml_solution_hessian" begin
-        solution = fit(SemOptimizer(engine = :Optim, algorithm = Newton()), model_ml)
-
+        solution = fit(
+            SemOptimizer(
+                engine = :Optim,
+                algorithm = Newton(linesearch = BackTracking(order = 3)),
+            ),
+            model_ml,
+        )
         update_estimate!(partable, solution)
         test_estimates(partable, solution_lav[:parameter_estimates_ml]; atol = 1e-2)
     end