Distributions - design points for discussion

[TomMonks]

A few points all related to the use of the distributions that we should discuss:

Distributions are also parameters

The design you have used (where distributions are part of the model class) is definitely valid (especially to get reproducible results), but in general the type of distribution used is a parameter that can be varied in a simulation model run. For example, we might switch from a Lognormal to a Exponential or Fixed distribution as part of model validation. Or we might switch from Lognormal to an Empirical distribution if lognormals are not suitable for different stroke pathway at a different trust (e.g when lognormal is a poor fit to the data).

At the moment it is possible to vary the parameters by monkey patching, but its not ideal to have to modify the model class.

I personally think the distributions should be decoupled from the model class itself. That way you can can modify the parameters however you need (one example is a subclassed Param class, but see below) and pass to the model which remains as intended with no dependencies on the type of distribution used (although potentially with routing, but maybe a half way how is to keep that type fixed as its model logic).

sim-tools distribution registry might help

I created the distribution registry for this very purpose. Although I've not yet had time to implement it in a model in anger. Take a look here https://tommonks.github.io/sim-tools/01_sampling/03_distributions_registry.html

That would allow us to have a more flexible model with less dependencies within the model logic.

@amyheather What do you think? I could if you wanted provide a simple example before we go any further?

[amyheather]

For my future reference, a note about changes in results:

(This is not an issue, just a change that will happen, and I wanted to record why).

In the original code, Model() had a create_distributions method:

def create_distributions(self, asu_param, rehab_param, distribution_type):
        """
        Create a nested dictionary with two items: "asu" and "rehab". Each
        then contains a set of distributions by patient type.

        Parameters
        ----------
        asu_param: Class
            Acute stroke unit (ASU) parameters.
        rehab_param: Class
            Rehabilitation unit parameters.
        distribution_type: str
            Name of the distribution to use ("exponential", "discrete",
            "lognormal").

        Returns
        -------
        distributions: dictionary
            Dictionary with "asu" and "rehab", which each then contains a
            dictionary of sampling distributions by patient type.
        """
        # Create dictionary to store distributions
        distributions = {}

        # Loop through each unit
        for unit, unit_param in [("asu", asu_param), ("rehab", rehab_param)]:

            # Make sub-dictionary to store that unit's distributions
            distributions[unit] = {}

            # Get a list of the patients in that unit (ignore other attributes)
            patient_types = [attr for attr in dir(unit_param) if attr in
                             ["stroke", "stroke_esd", "stroke_noesd",
                              "stroke_mortality", "tia", "neuro", "other"]]

            # For each patient type...
            for patient_type in patient_types:

                # Get the entry for that patient
                patient_param = getattr(unit_param, patient_type)

                # Create distributions of specified type
                if distribution_type == "exponential":
                    distributions[unit][patient_type] = Exponential(
                        mean=patient_param,
                        random_seed=next(self.seed_generator)
                    )
                elif distribution_type == "discrete":
                    distributions[unit][patient_type] = DiscreteEmpirical(
                        values=list(patient_param.keys()),
                        freq=list(patient_param.values()),
                        random_seed=next(self.seed_generator)
                    )
                elif distribution_type == "lognormal":
                    distributions[unit][patient_type] = Lognormal(
                        mean=patient_param["mean"],
                        stdev=patient_param["sd"],
                        random_seed=next(self.seed_generator)
                    )
        return distributions

When finding the relevant patient types for a given distribution, it ran dir(unit_param) and extracted patient names. dir() will return all the class attributes in alphabetical order, so the patient_types would be e.g. ["neuro", "other", "stroke", "tia"].

This would then dictate the order in which the distributions are created - and, hence, which seeds are used for each.

When switching to DistributionRegistry, the distributions are created in the order they are defined in parameters.json. Hence, back tests fail, as the seeds used are different.

The backtesting succeeds if you reorder param.dist_config to match the old code. Created with the help of Perplexity, replace:

# Create all the distributions
        flat_dist = DistributionRegistry.create_batch(
            config=self.param.dist_config, main_seed=self.run_number
        )

With:

# Create all the distributions
        import re

        # Orderings from your old code logic
        distribution_types = ["arrival", "routing", "los"]
        units = ["asu", "rehab"]
        # Patient type canonical order (as you did with dir(...)), or sorted list:
        patient_types = ["neuro", "other", "stroke", "stroke_esd", "stroke_mortality", "stroke_noesd", "tia"]

        # Function to parse and order a key, e.g. "asu_arrival_stroke"
        def old_code_order_key(key):
            # Example key format: "asu_arrival_stroke_noesd"
            m = re.match(r'^(asu|rehab)_(arrival|routing|los)_(.+)$', key)
            if not m:
                raise ValueError(f"Key {key} does not match expected format.")
            unit, dist_type, patient = m.groups()
            # Map to their indexes
            return (
                distribution_types.index(dist_type),
                units.index(unit),
                patient_types.index(patient) if patient in patient_types else 999
            )

        sorted_keys = sorted(self.param.dist_config.keys(), key=old_code_order_key)
        sorted_config = {k: self.param.dist_config[k] for k in sorted_keys}
        flat_dist = DistributionRegistry.create_batch(
            config=sorted_config, main_seed=self.run_number
        )

I will switch to the new order and update the results for back testing, satisifed that there is no underlying issues, and this is the cause of the change.

[amyheather]

Remaining to do:

• Address potential issue that changing the order of parameters.json changes seeds and so changes results. Could lead to accidentally attributing differences to a change that was actually down to seeds. Solution could be to have e.g. `create_batch(sort=True)". : See pull request DistributionRegistry and linting TomMonks/sim-tools#57
• Add validation - previously, I used RestrictAttributes to ensure there aren't silent errors from mis-naming when replacing parameters. Need to add similar for this. I also had validation of some of the values themselves (e.g. routing between 0 and 1)... but that might already be implemented by DiscreteEmpirical - to check.
• Add tests for any validation introduced (similar to the tests I used to have for the old validation methods).

Some validation locally (LockedDict), others in sim-tools: TomMonks/sim-tools#57