Merge pull request #101 from fjebaker/fergus/fixes

Variance fixes in simulated spectra
fjebaker · May 14, 2024 · 233cb5d · 233cb5d
2 parents f54f999 + 3e0a635
commit 233cb5d
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 13 deletions.
diff --git a/docs/src/walkthrough.md b/docs/src/walkthrough.md
@@ -130,7 +130,7 @@ update_model!(model, result)
 To estimate the goodness of our fit, we can mimic the `goodness` command from XSPEC. This will use the [`simulate`](@ref) function to simulate spectra for a dataset (here determined by the result), and fit the model to the simulated dataset. The fit statistic for each fit is then appended to an array, which we can use to plot a histogram:
 
 ```@example walk
-spread = goodness(result; N = 1000, seed = 42)
+spread = goodness(result; N = 1000, seed = 42, exposure_time = data.data.spectrum.exposure_time)
 histogram(spread, ylims = (0, 300), label = "Simulated")
 vline!([result.χ2], label = "Best fit")
 ```

diff --git a/src/simulate.jl b/src/simulate.jl
@@ -3,7 +3,7 @@ mutable struct SimulatedSpectrum{T,F} <: AbstractDataset
     model_domain::Vector{T}
     output_domain::Vector{T}
     data::Vector{T}
-    errors::Vector{T}
+    variance::Vector{T}
     units::Union{Nothing,SpectralUnits.RateOrCount}
     transformer!!::F
     seed::Int
@@ -18,7 +18,7 @@ end
 
 function make_objective_variance(::ContiguouslyBinned, dataset::SimulatedSpectrum)
     check_units_warning(dataset.units)
-    dataset.errors .^ 2
+    dataset.variance
 end
 
 make_model_domain(::ContiguouslyBinned, dataset::SimulatedSpectrum) = dataset.model_domain
@@ -51,13 +51,18 @@ function simulate!(
     simulate_distribution = Distributions.Poisson,
     rng = Random.default_rng(),
     exposure_time = 1e5,
+    counting_variance = true,
 )
     config.objective .= _invoke_and_transform!(config.cache, config.model_domain, p)
     for (i, m) in enumerate(config.objective)
         distr = simulate_distribution(m * exposure_time)
         count = rand(rng, distr)
         config.objective[i] = count / exposure_time
-        config.variance[i] = sqrt(abs(count)) / exposure_time
+
+        # how to propagate the variance
+        if counting_variance
+            config.variance[i] = count / exposure_time^2
+        end
     end
 end
 
@@ -69,7 +74,6 @@ function simulate!(conf::FittingConfig; seed = abs(rand(Int)), kwargs...)
         conf.model_domain,
         conf.output_domain,
         conf.objective,
-        # variance has already been sqrt'd
         conf.variance,
         u"counts / (s * keV)",
         conf.cache.transformer!!,

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -93,21 +93,21 @@ end
 end
 
 @testset "fitting" verbose = true begin
-    include("fitting/test-fit-simple-dataset.jl")
-    include("fitting/test-cache.jl")
-    include("fitting/test-binding.jl")
-    include("fitting/test-results.jl")
+    @time include("fitting/test-fit-simple-dataset.jl")
+    @time include("fitting/test-cache.jl")
+    @time include("fitting/test-binding.jl")
+    @time include("fitting/test-results.jl")
 
     @testset "powerlaws" begin
-        include("fitting/test-fit-powerlaw.jl")
+        @time include("fitting/test-fit-powerlaw.jl")
     end
     @testset "multifits" begin
-        include("fitting/test-fit-multi.jl")
-        include("fitting/test-fit-optim.jl")
+        @time include("fitting/test-fit-multi.jl")
+        @time include("fitting/test-fit-optim.jl")
     end
     if has_test_dir
         @testset "sample-data" begin
-            include("fitting/test-sample-data.jl")
+            @time include("fitting/test-sample-data.jl")
         end
     else
         @warn "Skipping dataset tests."