Fix DiffEqBayesLotkaVolterra: merge split imports, fix sample_kwargs type, switch to @elapsed

benchmarks/BayesianInference/DiffEqBayesLotkaVolterra.jmd

@@ -8,9 +8,6 @@ author: Vaibhav Dixit, Chris Rackauckas
 
 ```julia
 using DiffEqBayes, StanSample, DynamicHMC, Turing
-```
-
-```julia
 using Distributions, BenchmarkTools, StaticArrays
 using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions
 using Plots, LinearAlgebra
@@ -107,9 +104,9 @@ We use `adapt_delta = 0.85` (Stan's default) consistently across all backends fo
 Stan infers a separate noise parameter (`sigma`) per data dimension via the `vars` specification.
 
 ```julia
-bayesian_result_stan = @time stan_inference(
+elapsed_stan = @elapsed bayesian_result_stan = stan_inference(
     prob, :rk45, t, data, priors; print_summary = false,
-    sample_kwargs = Dict(:delta => 0.85, :num_samples => 10_000),
+    sample_kwargs = (delta = 0.85, num_samples = 10_000),
     vars = (DiffEqBayes.StanODEData(), InverseGamma(2, 3)))
 ```
 

benchmarks/BayesianInference/DiffEqBayesLotkaVolterraSDE.jmd

@@ -0,0 +1,169 @@
+---
+title: Lotka-Volterra SDE Bayesian Parameter Estimation Benchmarks
+author: Arpan Chakraborty, Chris Rackauckas
+---
+## Parameter Estimation of Stochastic Lotka-Volterra using Turing.jl
+For ODEs, the likelihood given parameters is exact. For SDEs, the true
+transition density p(X_{t+1} | X_t, θ) is generally intractable. We use
+the Euler-Maruyama pseudo-likelihood as an approximation:
+    p_EM(X_{t+1} | X_t, θ) ≈ N(X_t + f(X_t,θ)Δt, g(X_t,θ)²Δt)
+This benchmark compares Bayesian MCMC (Turing) against MLE optimization on
+the same pseudo-likelihood, asking whether the Bayesian approach gives better
+uncertainty quantification for stochastic dynamical systems.
+```julia
+using Turing, StochasticDiffEq
+using Distributions, BenchmarkTools, StaticArrays
+using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions
+using Plots, LinearAlgebra, Statistics, Printf
+gr(fmt = :png)
+```
+```julia
+"""Display ESS/s (effective samples per second) from a Turing chain."""
+function display_ess_per_sec(chain, elapsed)
+    stats = summarystats(chain)
+    ess_bulk = stats[:, :ess_bulk]
+    println("Elapsed time: $(round(elapsed; digits=2)) seconds\n")
+    println("ESS/s (effective samples per second, bulk):")
+    for (i, param) in enumerate(stats[:, :parameters])
+        println("  $param: $(round(ess_bulk[i] / elapsed; digits=1))")
+    end
+    println("\nMinimum ESS/s (bulk): $(round(minimum(ess_bulk) / elapsed; digits=1))")
+end
+```
+#### Initializing the problem
+The drift term is the standard Lotka-Volterra ODE. The diffusion term adds
+proportional noise to each species, scaled by σ (the 5th parameter).
+```julia
+f = @ode_def LotkaVolterraSDE begin
+    dx = a*x - b*x*y
+    dy = -c*y + d*x*y
+end a b c d
+function g(du, u, p, t)
+    du[1] = p[5] * u[1]
+    du[2] = p[5] * u[2]
+end
+```
+```julia
+u0    = [1.0, 1.0]
+tspan = (0.0, 10.0)
+p_true = [1.5, 1.0, 3.0, 1.0, 0.1]   # a, b, c, d, σ
+```
+```julia
+prob = SDEProblem(f, g, u0, tspan, p_true)
+sol  = solve(prob, SRIW1(), seed = 14)
+plot(sol, title = "Stochastic Lotka-Volterra (true parameters)",
+     label = ["prey" "predator"])
+```
+#### Generate noisy observations
+```julia
+t    = collect(range(1, stop = 10, length = 10))
+sig  = 0.49
+data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))
+```
+```julia
+scatter(t, data[1, :], lab = "#prey (data)")
+scatter!(t, data[2, :], lab = "#predator (data)")
+plot!(sol)
+```
+```julia
+priors = [truncated(Normal(1.5, 0.5), 0.5, 2.5),
+          truncated(Normal(1.2, 0.5), 0, 2),
+          truncated(Normal(3.0, 0.5), 1, 4),
+          truncated(Normal(1.0, 0.5), 0, 2)]
+```
+### Bayesian inference with Euler-Maruyama pseudo-likelihood
+Because the SDE transition density is intractable, we embed the SDE solver
+inside a Turing model and observe data through a Gaussian likelihood. The
+posterior reflects both parameter uncertainty and the intrinsic stochasticity of
+the process.
+```julia
+@model function fitlv_sde(data, t, prob)
+    a ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
+    b ~ truncated(Normal(1.2, 0.5), 0.0, 2.0)
+    c ~ truncated(Normal(3.0, 0.5), 1.0, 4.0)
+    d ~ truncated(Normal(1.0, 0.5), 0.0, 2.0)
+    σ ~ truncated(Normal(0.1, 0.05), 0.01, 0.5)
+    σ_obs ~ filldist(InverseGamma(2, 3), 2)
+    p = [a, b, c, d, σ]
+    _prob = remake(prob; p = p)
+    sol = solve(_prob, SRIW1(); saveat = t)
+    if sol.retcode != ReturnCode.Success || size(Array(sol), 2) != length(t)
+        Turing.@addlogprob! -Inf
+        return
+    end
+    for i in 1:length(t)
+        data[:, i] ~ MvNormal(sol[i], Diagonal(σ_obs .^ 2))
+    end
+end
+model = fitlv_sde(data, t, prob)
+sample(model, NUTS(0.65), 10; progress = false)   # warmup
+```
+```julia
+bayesian_result_turing = @btime sample($model, NUTS(0.65), 1000; progress = false)
+```
+```julia
+elapsed_bayes = @elapsed chain = sample(model, NUTS(0.65), 1000; progress = false)
+display_ess_per_sec(chain, elapsed_bayes)
+```
+### MLE optimization baseline
+The same Euler-Maruyama pseudo-likelihood maximized to a point estimate.
+No uncertainty is expressed — the bias from the discretization is locked in.
+```julia
+using Optimization, OptimizationOptimJL
+function neg_log_pseudo_likelihood(θ, _)
+    a, b, c, d, σ, σ1, σ2 = θ
+    any(x -> x <= 0, θ) && return Inf
+    p     = [a, b, c, d, σ]
+    _prob = remake(prob; p = p)
+    sol   = solve(_prob, SRIW1(); saveat = t)
+    (sol.retcode != ReturnCode.Success || size(Array(sol), 2) != length(t)) && return Inf
+    ll = 0.0
+    for i in 1:length(t)
+        ll += logpdf(MvNormal(sol[i], Diagonal([σ1, σ2] .^ 2)), data[:, i])
+    end
+    return -ll
+end
+θ0       = [1.5, 1.2, 3.0, 1.0, 0.1, 0.5, 0.5]
+opt      = OptimizationFunction(neg_log_pseudo_likelihood, Optimization.AutoFiniteDiff())
+opt_prob = OptimizationProblem(opt, θ0)
+elapsed_mle = @elapsed opt_result = solve(opt_prob, NelderMead())
+mle_params  = opt_result.u
+@printf "MLE: a=%.3f  b=%.3f  c=%.3f  d=%.3f\n" mle_params[1] mle_params[2] mle_params[3] mle_params[4]
+```
+### Comparison: posterior vs point estimate
+```julia
+param_names = [:a, :b, :c, :d]
+true_vals   = [1.5, 1.0, 3.0, 1.0]
+bayes_means = [mean(chain[p]) for p in param_names]
+bayes_stds  = [std(chain[p])  for p in param_names]
+mle_vals    = mle_params[1:4]
+@printf "\n%-6s  %-5s  %-22s  %-6s\n" "Param" "True" "Bayes mean ± std" "MLE"
+@printf "%s\n" "-"^48
+for i in 1:4
+    @printf "%-6s  %-5.2f  %.3f ± %.3f %14s  %.3f\n" \
+        string(param_names[i]) true_vals[i] bayes_means[i] bayes_stds[i] "" mle_vals[i]
+end
+@printf "\nBayesian MCMC : %.1f s\nMLE           : %.3f s\n" elapsed_bayes elapsed_mle
+```
+```julia
+pl = plot(layout = (2, 2), size = (800, 600))
+for (i, param) in enumerate(param_names)
+    histogram!(pl[i], vec(Array(chain[param])), normalize = true,
+               label = "Posterior", alpha = 0.6, title = string(param))
+    vline!(pl[i], [true_vals[i]], label = "True",  lw = 2, color = :green)
+    vline!(pl[i], [mle_vals[i]], label = "MLE",   lw = 2, color = :red, linestyle = :dash)
+end
+plot!(pl, legend = :topright)
+```
+## Conclusion
+For SDE parameter estimation, the Bayesian pseudo-likelihood approach
+provides full posterior distributions over parameters, naturally quantifying
+uncertainty from both measurement noise and the intrinsic stochasticity of
+the process. The MLE baseline recovers similar point estimates but cannot
+express this uncertainty — and as σ increases or the dataset shrinks, the
+Bayesian posterior widens appropriately, a property point estimates
+fundamentally lack.
+```julia, echo = false
+using SciMLBenchmarks
+SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
+```