diff --git a/ChangeLog.md b/ChangeLog.md index 5e4500be7176396cebc480f119ab31b0409954f2..56af74a738e9a7bd32381764b5c3b1f713dc86d3 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -2,6 +2,9 @@ ## 0.1.3.0 +- The aggregation simulation study now accepts a random seed from the + configuration file and assumes there will be a positive rate for unscheduled + unsequenced sampling. - **Important** The aggregation simulation study now ignores any occurrence events that occur after the last unscheduled sequenced sample. - **Important** Change the behaviour of the aggregation functions to remove any diff --git a/apps/simulation-study-aggregated-observations/Main.hs b/apps/simulation-study-aggregated-observations/Main.hs index 42aa9dd5e11e946e76281eec4fbf9d17ca6e3572..7a14871cba36b6c9aaf09abb6b64a0a34b02fb14 100644 --- a/apps/simulation-study-aggregated-observations/Main.hs +++ b/apps/simulation-study-aggregated-observations/Main.hs @@ -59,12 +59,16 @@ import Epidemic.Types.Parameter (Probability, Rate, Time, Timed(..)) import Epidemic.Types.Population (Person(..)) import qualified Epidemic.Utility as SimUtil import GHC.Generics +import GHC.Word (Word32(..)) import Numeric.GSL.Minimization (MinimizeMethod(NMSimplex2), minimizeV) import Numeric.LinearAlgebra (dot) import Numeric.LinearAlgebra.Data (Vector(..), fromList, linspace, toList) import System.Environment (getArgs) import System.Random.MWC (initialize) +-- | Alias for the type used to seed the MWC PRNG. +type MWCSeed = Word32 + -- | These objects define the specifics of the evaluation of LLHD profiles. If a -- point estimate is given, then that is the central point of the profiles, -- otherwise the parameters are estimated first. In every case the natural death @@ -106,6 +110,7 @@ data Configuration = , InferenceConfiguration , InferenceConfiguration) , isVerbose :: Bool + , mwcSeed :: MWCSeed } deriving (Show, Generic) @@ -137,7 +142,7 @@ data AnnotatedParameter bdscodConfiguration = do simParams@(Parameters (pLambda, pMu, pPsi, Timed pRhos, pOmega, Timed pNus)) <- asks simulationParameters - if pLambda > 0 && pMu > 0 && pPsi > 0 && null pRhos && null pNus + if pLambda > 0 && pMu > 0 && pPsi > 0 && pOmega > 0 && null pRhos && null pNus then do simDur <- asks simulationDuration let bdscodConfig = @@ -386,8 +391,8 @@ estimateAggregatedParameters deathRate (rhoTimes,nuTimes) obs = simulationStudy :: Simulation () simulationStudy = do bdscodConfig <- bdscodConfiguration - let seedInt = 42 - epiSim <- simulateEpidemic seedInt bdscodConfig + prngSeed <- asks mwcSeed + epiSim <- simulateEpidemic prngSeed bdscodConfig (regObs, regObs', aggObs) <- observeEpidemicThrice epiSim uncurry evaluateLLHD regObs uncurry estimateLLHD regObs' diff --git a/examples/simulation-study-aggregated-observations/README.org b/examples/simulation-study-aggregated-observations/README.org index 12e55dc070cdef195c611ed53d1d5d0be6113dcd..35188fb4d9fccfc8b5a694f7dce6960586d8a6bd 100644 --- a/examples/simulation-study-aggregated-observations/README.org +++ b/examples/simulation-study-aggregated-observations/README.org @@ -8,7 +8,9 @@ aggregated into scheduled observations. The running of this simulation study is coordinated by the script =run.sh=. The contents of that file are given below. The simulation study is specified in a JSON file which is given as a command line argument. The R script -=make-json-config.R= produces this so you don't need to edit it manually. +=make-json-config.R= produces this so you don't need to edit it manually. Note +that this also contains the seed for the random number generator so you control +this from outside of the application. =run.sh= #+BEGIN_SRC sh :tangle run.sh @@ -23,3 +25,7 @@ Rscript src/plot.R Note that there is an assumption here that every occurrence event that occurs after the last unscheduled sequenced sample is removed before the aggregation step. + +The =src/plot.R= script generates all of the figures and writes them to =out/=. +The most interesting figure is =regular-and-aggregated-data.png= which contrasts +the prevalence estimate using the regular and aggregated data. diff --git a/examples/simulation-study-aggregated-observations/src/make-json-config.R b/examples/simulation-study-aggregated-observations/src/make-json-config.R index 4c9f4edef6b3a28944fcf2ffd56616fec9a4522f..355e7e0e6fb9d534a9f57e4abb7acbfb92a056fa 100644 --- a/examples/simulation-study-aggregated-observations/src/make-json-config.R +++ b/examples/simulation-study-aggregated-observations/src/make-json-config.R @@ -13,7 +13,7 @@ simulation_duration <- 8.5 - 1e-6 birth_rate <- 1.7 death_rate <- 0.5 -sampling_rate <- 0.15 +sampling_rate <- 0.2 occurrence_rate <- 0.3 disaster_params <- list() @@ -22,7 +22,7 @@ catastrophe_params <- list() ## For the aggregation, these are the times at which we carry out the ## aggregation. seq_agg_times <- as.list(seq(from = 2.5, to = 8.5, by = 1)) -unseq_agg_times <- as.list(seq(from = 2.0, to = 8.0, by = 1)) +unseq_agg_times <- as.list(seq(from = 2.4, to = 8.4, by = 1)) @@ -82,11 +82,16 @@ result <- list( inferenceConfigurations = list( inference_configuration("true-params-regular-data", NULL), inference_configuration("est-params-regular-data", NULL), - inference_configuration("est-params-agg-data", - list(seq_agg_times, - unseq_agg_times)) + inference_configuration( + "est-params-agg-data", + list( + seq_agg_times, + unseq_agg_times + ) + ) ), - isVerbose = TRUE + isVerbose = TRUE, + mwcSeed = 56 ) jsonlite::write_json(result, @@ -95,21 +100,3 @@ jsonlite::write_json(result, auto_unbox = TRUE, digits = 7 ) - - - -## save a copy of the true parameters so they can be read out later rather than -## hardcoded. -## -## TODO convert this to using a JSON output rather than a CSV because it is -## neater. -## -true_parameters <- data.frame( - parameter = c("lambda", "mu"), - value = c(birth_rate, death_rate) -) -write.table( - x = true_parameters, - file = "out/true-parameters.csv", - row.names = FALSE -)