Skip to content

Monte Carlo

Uncertainty and sensitivity analysis via Monte Carlo simulation. Wrap any motor parameter with MonteCarloParameter to assign a statistical distribution (normal or uniform) and spread. MonteCarloSimulation then generates N randomized scenarios, runs the internal ballistics for each, and provides statistical aggregates (mean, std, percentiles, skew, kurtosis) over the results.

Useful for characterizing thrust variability, pressure envelope, and impulse confidence intervals due to manufacturing tolerances and propellant batch variation.

machwave.montecarlo

MonteCarloParameter dataclass

Bases: float

Float value annotated with a distribution and spread for sampling.

Source code in machwave/montecarlo/base.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
@dataclass(init=False, slots=True)
class MonteCarloParameter(float):
    """Float value annotated with a distribution and spread for sampling."""

    _spread: float | tuple[float, float]
    _distribution: str

    def __new__(cls, value: float, *, spread=0.0, distribution="normal"):
        """Construct a parameter with a center value, spread, and distribution."""
        obj = float.__new__(cls, value)
        obj._spread = spread
        obj._distribution = distribution
        return obj

    @property
    def random_generator(self) -> random.RandomGenerator:
        """Return a generator that samples from this parameter's distribution."""
        return random.get_random_generator(
            self._distribution, float(self), self._spread
        )

    def get_random_value(self) -> float:
        """Sample a single random value from this parameter's distribution."""
        return self.random_generator.get_value()

    def __repr__(self) -> str:
        """Return a debugging representation of the parameter."""
        return (
            f"MonteCarloParameter({float(self):.6g}, spread={self._spread}, "
            f"dist='{self._distribution}')"
        )

random_generator property

Return a generator that samples from this parameter's distribution.

__new__(value, *, spread=0.0, distribution='normal')

Construct a parameter with a center value, spread, and distribution.

Source code in machwave/montecarlo/base.py 
24
25
26
27
28
29
def __new__(cls, value: float, *, spread=0.0, distribution="normal"):
    """Construct a parameter with a center value, spread, and distribution."""
    obj = float.__new__(cls, value)
    obj._spread = spread
    obj._distribution = distribution
    return obj

__repr__()

Return a debugging representation of the parameter.

Source code in machwave/montecarlo/base.py 
42
43
44
45
46
47
def __repr__(self) -> str:
    """Return a debugging representation of the parameter."""
    return (
        f"MonteCarloParameter({float(self):.6g}, spread={self._spread}, "
        f"dist='{self._distribution}')"
    )

get_random_value()

Sample a single random value from this parameter's distribution.

Source code in machwave/montecarlo/base.py 
38
39
40
def get_random_value(self) -> float:
    """Sample a single random value from this parameter's distribution."""
    return self.random_generator.get_value()

MonteCarloSimulation

Monte Carlo driver that runs scenarios and stores their results.

Stores parameters and simulation type, executes the configured number of scenarios, and delegates plotting to plot_service.py.

Source code in machwave/montecarlo/base.py 
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
class MonteCarloSimulation:
    """
    Monte Carlo driver that runs scenarios and stores their results.

    Stores parameters and simulation type, executes the configured number of
    scenarios, and delegates plotting to `plot_service.py`.
    """

    def __init__(
        self,
        parameters: list[typing.Any],
        number_of_scenarios: int,
        simulation: type[machwave_simulation.InternalBallisticsSimulation],
    ) -> None:
        """
        Initialize a Monte Carlo driver.

        Args:
            parameters: Input parameters for a simulation class.
            number_of_scenarios: Number of scenarios to simulate.
            simulation: Simulation class reference.
        """
        self.parameters = parameters
        self.number_of_scenarios = number_of_scenarios
        self.simulation = simulation

        self.scenarios: list = []
        self.results: list = []

        self._object_store = dict()  # for nested parameter processing

    def generate_scenario(self) -> list[typing.Any]:
        """Generates a Monte Carlo scenario in the form of a list of parameters."""
        new_scenario = []
        parameters_copy = deepcopy(self.parameters)

        for parameter in parameters_copy:
            if isinstance(parameter, MonteCarloParameter):
                parameter = parameter.get_random_value()
            else:
                self._process_nested_parameters(parameter)

            new_scenario.append(parameter)

        self.scenarios.append(new_scenario)
        return new_scenario

    def _process_nested_parameters(self, parameter: typing.Any) -> None:
        """
        Replace nested `MonteCarloParameter` instances with randomized values.

        Recursively walks the object's attributes; replacements are performed
        in place and intermediate objects are tracked using UUIDs.

        Args:
            parameter: Object whose attributes will be processed.
        """
        parameter_uuid = uuid.uuid4()
        self._object_store[parameter_uuid] = parameter
        search_tree = {parameter_uuid: common_objects.get_object_dict(parameter)}

        i = 0
        while search_tree and i < SEARCH_TREE_DEPTH_LIMIT:
            i += 1
            new_search_tree = {}
            for param_uuid, sub_params in search_tree.items():
                param = self._object_store[param_uuid]
                for name, attr in sub_params.items():
                    object_uuid = uuid.uuid4()

                    if isinstance(attr, MonteCarloParameter):
                        setattr(param, name, attr.get_random_value())
                    elif isinstance(attr, list):
                        for item in attr:
                            if isinstance(item, dict):
                                continue
                            self._object_store[object_uuid] = item
                            new_search_tree[object_uuid] = (
                                common_objects.get_object_dict(item)
                            )
                    else:
                        object_uuid = uuid.uuid4()
                        self._object_store[object_uuid] = attr
                        new_search_tree[object_uuid] = common_objects.get_object_dict(
                            attr
                        )

            search_tree = new_search_tree

    def run(self) -> None:
        """Executes `number_of_scenarios` runs of the underlying Simulation."""
        self.results = []
        for _ in range(self.number_of_scenarios):
            scenario = self.generate_scenario()
            self.results.append(self.simulation(*scenario).run())

    def retrieve_values_from_result(self, property_name: str) -> np.ndarray:
        """
        Retrieve a scalar property across every Monte Carlo result.

        Args:
            property_name: Attribute name on the `SimulationResult`.

        Returns:
            Array of length `number_of_scenarios`.
        """
        return np.array(
            [getattr(sim_result, property_name) for sim_result in self.results]
        )

    def get_property_stats(self, property_name: str) -> dict[str, float]:
        """
        Return descriptive statistics for a scalar property across results.

        Metrics returned: mean, median, variance, std_dev, mode, skew (Fisher,
        unbiased), kurtosis (excess, unbiased), p5 (5th percentile), p95 (95th
        percentile).

        Args:
            property_name: Attribute name on the `SimulationResult`.

        Returns:
            Mapping of metric name to value.
        """
        values = np.asarray(self.retrieve_values_from_result(property_name))

        mean_val = np.mean(values)
        median_val = np.median(values)
        var_val = np.var(values)
        std_val = np.std(values)

        mode_val = float(scipy_stats.mode(values, nan_policy="omit").mode)
        skew_val = scipy_stats.skew(values, bias=False)  # unbiased Fisher skew
        kurt_val = scipy_stats.kurtosis(values, fisher=True, bias=False)
        p5, p95 = np.percentile(values, [5, 95])

        return {
            "mean": float(mean_val),
            "median": float(median_val),
            "variance": float(var_val),
            "std_dev": float(std_val),
            "mode": float(mode_val),
            "skew": float(skew_val),
            "kurtosis": float(kurt_val),
            "p5": float(p5),
            "p95": float(p95),
        }

    def plot_histogram(
        self,
        property_name: str,
        x_axes_title: str = "x",
        **plotly_kwargs,
    ) -> None:
        """Plot a histogram of a scalar property across all results."""
        plot_service.plot_histogram(
            self.results, property_name, x_axes_title, **plotly_kwargs
        )

    def plot_histogram_with_kde(
        self,
        property_name: str,
        x_axes_title: str = "x",
        nbins: int = 30,
        kde_points: int = 200,
        **plotly_kwargs,
    ) -> None:
        """Plot a histogram of a scalar property with a KDE overlay."""
        plot_service.plot_histogram_with_kde(
            self.results,
            property_name,
            x_axes_title,
            nbins,
            kde_points,
            **plotly_kwargs,
        )

    def plot_cdf(
        self,
        property_name: str,
        x_axes_title: str = "x",
        **plotly_kwargs,
    ) -> None:
        """Plot the empirical CDF of a scalar property across all results."""
        plot_service.plot_cdf(
            self.results, property_name, x_axes_title, **plotly_kwargs
        )

    def plot_time_series_extremes(
        self,
        time_property: str,
        series_property: str,
        x_axes_title: str = "time",
        title: str | None = None,
        **plotly_kwargs,
    ) -> None:
        """Plot min and max envelopes of a time series across all results."""
        plot_service.plot_time_series_extremes(
            self.results,
            time_property,
            series_property,
            x_axes_title,
            title,
            **plotly_kwargs,
        )

__init__(parameters, number_of_scenarios, simulation)

Initialize a Monte Carlo driver.

Parameters:

Name Type Description Default
parameters list[Any]

Input parameters for a simulation class.

 required
number_of_scenarios int

Number of scenarios to simulate.

 required
simulation type[InternalBallisticsSimulation]

Simulation class reference.

 required
Source code in machwave/montecarlo/base.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def __init__(
    self,
    parameters: list[typing.Any],
    number_of_scenarios: int,
    simulation: type[machwave_simulation.InternalBallisticsSimulation],
) -> None:
    """
    Initialize a Monte Carlo driver.

    Args:
        parameters: Input parameters for a simulation class.
        number_of_scenarios: Number of scenarios to simulate.
        simulation: Simulation class reference.
    """
    self.parameters = parameters
    self.number_of_scenarios = number_of_scenarios
    self.simulation = simulation

    self.scenarios: list = []
    self.results: list = []

    self._object_store = dict()  # for nested parameter processing

generate_scenario()

Generates a Monte Carlo scenario in the form of a list of parameters.

Source code in machwave/montecarlo/base.py 
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def generate_scenario(self) -> list[typing.Any]:
    """Generates a Monte Carlo scenario in the form of a list of parameters."""
    new_scenario = []
    parameters_copy = deepcopy(self.parameters)

    for parameter in parameters_copy:
        if isinstance(parameter, MonteCarloParameter):
            parameter = parameter.get_random_value()
        else:
            self._process_nested_parameters(parameter)

        new_scenario.append(parameter)

    self.scenarios.append(new_scenario)
    return new_scenario

get_property_stats(property_name)

Return descriptive statistics for a scalar property across results.

Metrics returned: mean, median, variance, std_dev, mode, skew (Fisher, unbiased), kurtosis (excess, unbiased), p5 (5th percentile), p95 (95th percentile).

Parameters:

Name Type Description Default
property_name str

Attribute name on the SimulationResult.

 required

Returns:

Type Description
dict[str, float]

Mapping of metric name to value.
Source code in machwave/montecarlo/base.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
def get_property_stats(self, property_name: str) -> dict[str, float]:
    """
    Return descriptive statistics for a scalar property across results.

    Metrics returned: mean, median, variance, std_dev, mode, skew (Fisher,
    unbiased), kurtosis (excess, unbiased), p5 (5th percentile), p95 (95th
    percentile).

    Args:
        property_name: Attribute name on the `SimulationResult`.

    Returns:
        Mapping of metric name to value.
    """
    values = np.asarray(self.retrieve_values_from_result(property_name))

    mean_val = np.mean(values)
    median_val = np.median(values)
    var_val = np.var(values)
    std_val = np.std(values)

    mode_val = float(scipy_stats.mode(values, nan_policy="omit").mode)
    skew_val = scipy_stats.skew(values, bias=False)  # unbiased Fisher skew
    kurt_val = scipy_stats.kurtosis(values, fisher=True, bias=False)
    p5, p95 = np.percentile(values, [5, 95])

    return {
        "mean": float(mean_val),
        "median": float(median_val),
        "variance": float(var_val),
        "std_dev": float(std_val),
        "mode": float(mode_val),
        "skew": float(skew_val),
        "kurtosis": float(kurt_val),
        "p5": float(p5),
        "p95": float(p95),
    }

plot_cdf(property_name, x_axes_title='x', **plotly_kwargs)

Plot the empirical CDF of a scalar property across all results.

Source code in machwave/montecarlo/base.py 
227
228
229
230
231
232
233
234
235
236
def plot_cdf(
    self,
    property_name: str,
    x_axes_title: str = "x",
    **plotly_kwargs,
) -> None:
    """Plot the empirical CDF of a scalar property across all results."""
    plot_service.plot_cdf(
        self.results, property_name, x_axes_title, **plotly_kwargs
    )

plot_histogram(property_name, x_axes_title='x', **plotly_kwargs)

Plot a histogram of a scalar property across all results.

Source code in machwave/montecarlo/base.py 
198
199
200
201
202
203
204
205
206
207
def plot_histogram(
    self,
    property_name: str,
    x_axes_title: str = "x",
    **plotly_kwargs,
) -> None:
    """Plot a histogram of a scalar property across all results."""
    plot_service.plot_histogram(
        self.results, property_name, x_axes_title, **plotly_kwargs
    )

plot_histogram_with_kde(property_name, x_axes_title='x', nbins=30, kde_points=200, **plotly_kwargs)

Plot a histogram of a scalar property with a KDE overlay.

Source code in machwave/montecarlo/base.py 
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
def plot_histogram_with_kde(
    self,
    property_name: str,
    x_axes_title: str = "x",
    nbins: int = 30,
    kde_points: int = 200,
    **plotly_kwargs,
) -> None:
    """Plot a histogram of a scalar property with a KDE overlay."""
    plot_service.plot_histogram_with_kde(
        self.results,
        property_name,
        x_axes_title,
        nbins,
        kde_points,
        **plotly_kwargs,
    )

plot_time_series_extremes(time_property, series_property, x_axes_title='time', title=None, **plotly_kwargs)

Plot min and max envelopes of a time series across all results.

Source code in machwave/montecarlo/base.py 
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
def plot_time_series_extremes(
    self,
    time_property: str,
    series_property: str,
    x_axes_title: str = "time",
    title: str | None = None,
    **plotly_kwargs,
) -> None:
    """Plot min and max envelopes of a time series across all results."""
    plot_service.plot_time_series_extremes(
        self.results,
        time_property,
        series_property,
        x_axes_title,
        title,
        **plotly_kwargs,
    )

retrieve_values_from_result(property_name)

Retrieve a scalar property across every Monte Carlo result.

Parameters:

Name Type Description Default
property_name str

Attribute name on the SimulationResult.

 required

Returns:

Type Description
ndarray

Array of length number_of_scenarios.
Source code in machwave/montecarlo/base.py 
146
147
148
149
150
151
152
153
154
155
156
157
158
def retrieve_values_from_result(self, property_name: str) -> np.ndarray:
    """
    Retrieve a scalar property across every Monte Carlo result.

    Args:
        property_name: Attribute name on the `SimulationResult`.

    Returns:
        Array of length `number_of_scenarios`.
    """
    return np.array(
        [getattr(sim_result, property_name) for sim_result in self.results]
    )

run()

Executes number_of_scenarios runs of the underlying Simulation.

Source code in machwave/montecarlo/base.py 
139
140
141
142
143
144
def run(self) -> None:
    """Executes `number_of_scenarios` runs of the underlying Simulation."""
    self.results = []
    for _ in range(self.number_of_scenarios):
        scenario = self.generate_scenario()
        self.results.append(self.simulation(*scenario).run())