Skip to content

API Reference

This API reference provides detailed documentation for the classes and functions in the pybmc package. It is automatically generated from the docstrings in the source code.

Dataset

A general-purpose dataset class for loading and managing model data for Bayesian model combination workflows.

Supports .h5 and .csv files, and provides data splitting functionality.

Source code in pybmc/data.py 
 10
 11
 12
 13
 14
 15
 16
 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
 48
 49
 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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
class Dataset:
    """
    A general-purpose dataset class for loading and managing model data
    for Bayesian model combination workflows.

    Supports .h5 and .csv files, and provides data splitting functionality.
    """

    def __init__(self, data_source=None, verbose=True):
        """
        Initialize the Dataset object.

        :param data_source: Path to the data file (.h5 or .csv).
        :param verbose: If True, display warnings and informational messages. Default is True.
        """
        self.data_source = data_source
        self.data = {}  # Dictionary of model to DataFrame
        self.verbose = verbose
        self.logger = logging.getLogger(__name__)
        if not self.logger.handlers:
            handler = logging.StreamHandler()
            handler.setFormatter(logging.Formatter('%(message)s'))
            self.logger.addHandler(handler)
        self.logger.setLevel(logging.INFO if verbose else logging.WARNING)

    def load_data(self, models, keys=None, domain_keys=None, model_column='model', truth_column_name=None):
        """
        Load data for each property and return a dictionary of synchronized DataFrames.
        Each DataFrame has columns: domain_keys + one column per model for that property.

        Parameters:
            models (list): List of model names (for HDF5 keys or filtering CSV).
            keys (list): List of property names to extract (each will be a key in the output dict).
            domain_keys (list, optional): List of columns used to define the common domain (default ['N', 'Z']).
            model_column (str, optional): Name of the column in the CSV that identifies which model each row belongs to.
                                          Only used for CSV files; ignored for HDF5 files.
            truth_column_name (str, optional): Name of the truth model. If provided, the truth data will be 
                                               left-joined to the common domain of the other models, allowing 
                                               the truth data to have a smaller domain than the models.

        Returns:
            dict: Dictionary where each key is a property name and each value is a DataFrame with columns:
                  domain_keys + one column per model for that property.
                  The DataFrames are synchronized to the intersection of the domains for all models.
                  If truth_column_name is provided, truth data is left-joined (may have NaN values).

        Supports both .h5 and .csv files.
        """
        self.domain_keys = domain_keys 

        if self.data_source is None:
            raise ValueError("Data source must be specified.")
        if not os.path.exists(self.data_source):
            raise FileNotFoundError(f"Data source '{self.data_source}' not found.")
        if keys is None:
            raise ValueError("You must specify which properties to extract via 'keys'.")

        result = {}

        for prop in keys:
            dfs = []
            truth_df = None
            skipped_models = []

            # Separate regular models from truth model
            regular_models = [m for m in models if m != truth_column_name]

            if self.data_source.endswith('.h5'):
                for model in models:
                    df = pd.read_hdf(self.data_source, key=model)
                    # Check required columns
                    missing_cols = [col for col in domain_keys + [prop] if col not in df.columns]
                    if missing_cols:
                        self.logger.info(f"[Skipped] Model '{model}' missing columns {missing_cols} for property '{prop}'.")
                        skipped_models.append(model)
                        continue
                    temp = df[domain_keys + [prop]].copy()
                    temp.rename(columns={prop: model}, inplace=True) # type: ignore

                    # Store truth data separately if truth_column_name is provided
                    if truth_column_name and model == truth_column_name:
                        truth_df = temp
                    else:
                        dfs.append(temp)

            elif self.data_source.endswith('.csv'):
                df = pd.read_csv(self.data_source)
                for model in models:
                    if model_column not in df.columns:
                        raise ValueError(f"Expected column '{model_column}' not found in CSV.")
                    model_df = df[df[model_column] == model]
                    missing_cols = [col for col in domain_keys + [prop] if col not in model_df.columns]
                    if missing_cols:
                        self.logger.info(f"[Skipped] Model '{model}' missing columns {missing_cols} for key '{prop}'.")
                        skipped_models.append(model)
                        continue
                    temp = model_df[domain_keys + [prop]].copy()
                    temp.rename(columns={prop: model}, inplace=True)

                    # Store truth data separately if truth_column_name is provided
                    if truth_column_name and model == truth_column_name:
                        truth_df = temp
                    else:
                        dfs.append(temp)
            else:
                raise ValueError("Unsupported file format. Only .h5 and .csv are supported.")

            if not dfs:
                self.logger.info(f"[Warning] No models with property '{prop}'. Resulting DataFrame will be empty.")
                result[prop] = pd.DataFrame(columns=domain_keys + [m for m in models if m not in skipped_models])
                continue

            # Intersect domain for regular models only
            common_df = dfs[0]
            for other_df in dfs[1:]:
                common_df = pd.merge(common_df, other_df, on=domain_keys, how="inner")
            # Drop rows with NaN in any required column (for regular models)
            common_df = common_df.dropna()

            # Left join truth data if it exists and was specified
            if truth_df is not None:
                common_df = pd.merge(common_df, truth_df, on=domain_keys, how="left")

            result[prop] = common_df
            self.data = result
        return result

    def view_data(self, property_name=None, model_name=None):
        """
        View data flexibly based on input parameters.

        - No arguments: returns available property names and model names.
        - property_name only: returns the full DataFrame for that property.
        - model_name only: Return model values across all properties.
        - property_name + model_name: returns a Series of values for the model.

        :param property_name: Optional property name 
        :param model_name: Optional model name 
        :return: dict, DataFrame, or Series depending on input.
        """

        if not self.data:
            raise RuntimeError("No data loaded. Run `load_data(...)` first.")

        if property_name is None and model_name is None:
            props = list(self.data.keys())
            models = sorted(set(col for prop_df in self.data.values() for col in prop_df.columns if col not in self.domain_keys))

            return {
                "available_properties": props,
                "available_models": models
            }

        if model_name is not None and property_name is None:
            # Return a dictionary: {property: Series of model values}
            result = {}
            for prop, df in self.data.items():
                if model_name in df.columns:
                    cols = self.domain_keys + [model_name]
                    result[prop] = df[cols]
                else:
                    result[prop] = f"[Model '{model_name}' not available]"
            return result

        if property_name is not None:
            if property_name not in self.data:
                raise KeyError(f"Property '{property_name}' not found.")

            df = self.data[property_name]

            if model_name is None:
                return df  # Full property DataFrame

            if model_name not in df.columns:
                raise KeyError(f"Model '{model_name}' not found in property '{property_name}'.")

            return df[model_name]



    def separate_points_distance_allSets(self, list1, list2, distance1, distance2):
            """
            Separates points in list1 into three groups based on their proximity to any point in list2.

            :param list1: List of (x, y) tuples.
            :param list2: List of (x, y) tuples.
            :param distance: The threshold distance to determine proximity.
            :return: Two lists - close_points and distant_points.
            """
            train = []
            validation=[]
            test = []

            train_list_coordinates=[]
            validation_list_coordinates=[]
            test_list_coordinates=[]

            for i in range(len(list1)):
                point1=list1[i]
                close = False
                for point2 in list2:
                    if np.linalg.norm(np.array(point1) - np.array(point2)) <= distance1:
                        close = True
                        break
                if close:
                    train.append(point1)
                    train_list_coordinates.append(i)
                else:
                    close2=False
                    for point2 in list2:
                        if np.linalg.norm(np.array(point1) - np.array(point2)) <= distance2:
                            close2 = True
                            break
                    if close2:
                        validation.append(point1)
                        validation_list_coordinates.append(i)
                    else:
                        test.append(point1)
                        test_list_coordinates.append(i)                

            return train_list_coordinates, validation_list_coordinates, test_list_coordinates

    def split_data(self, data_dict, property_name, splitting_algorithm="random", **kwargs):
        """
        Split data into training, validation, and testing sets using random or inside-to-outside logic.

        :param data_dict: Dictionary output from `load_data`, where keys are property names and values are DataFrames.
        :param property_name: The key in `data_dict` specifying which DataFrame to use for splitting.
        :param splitting_algorithm: 'random' (default) or 'inside_to_outside'.
        :param kwargs: Additional arguments depending on the chosen algorithm.
            For 'random': train_size, val_size, test_size
            For 'inside_to_outside': stable_points (list of (x, y)), distance1, distance2
        :return: Tuple of train, validation, test datasets as DataFrames.
        """
        if property_name not in data_dict:
            raise ValueError(f"Property '{property_name}' not found in the provided data dictionary.")

        data = data_dict[property_name]

        if isinstance(data, pd.DataFrame):
            indexable_data = data.reset_index(drop=True)
            point_list = list(indexable_data.itertuples(index=False, name=None))
        else:
            raise TypeError("Data for the specified property must be a pandas DataFrame.")

        if splitting_algorithm == "random":
            required = ['train_size', 'val_size', 'test_size']
            if not all(k in kwargs for k in required):
                raise ValueError(f"Missing required kwargs for 'random': {required}")

            train_size = kwargs['train_size']
            val_size = kwargs['val_size']
            test_size = kwargs['test_size']

            if not np.isclose(train_size + val_size + test_size, 1.0):
                raise ValueError("train_size + val_size + test_size must equal 1.0")

            # Random split using indexes
            train_idx, temp_idx = train_test_split(indexable_data.index, train_size=train_size, random_state=1)
            val_rel = val_size / (val_size + test_size)
            val_idx, test_idx = train_test_split(temp_idx, test_size=1 - val_rel, random_state=1)

        elif splitting_algorithm == "inside_to_outside":
            required = ['stable_points', 'distance1', 'distance2']
            if not all(k in kwargs for k in required):
                raise ValueError(f"Missing required kwargs for 'inside_to_outside': {required}")

            stable_points = kwargs['stable_points']
            distance1 = kwargs['distance1']
            distance2 = kwargs['distance2']

            train_idx, val_idx, test_idx = self.separate_points_distance_allSets(
                point_list, stable_points, distance1, distance2
            )
        else:
            raise ValueError("splitting_algorithm must be either 'random' or 'inside_to_outside'")

        train_data = indexable_data.iloc[train_idx]
        val_data = indexable_data.iloc[val_idx]
        test_data = indexable_data.iloc[test_idx]

        return train_data, val_data, test_data


    def get_subset(self, property_name, filters=None, models_to_include=None):
        """
        Return a filtered, wide-format DataFrame for a given property.

        :param property_name: Name of the property (e.g., "BE", "ChRad").
        :param filters: Dictionary of filtering rules applied to the domain columns (e.g., {"Z": (26, 28)}).
        :param models_to_include: Optional list of model names to retain in the output.
                                If None, all model columns are retained.
        :return: Filtered wide-format DataFrame with columns: domain keys + model columns.
        """
        if property_name not in self.data:
            raise ValueError(f"Property '{property_name}' not found in dataset.")

        df = self.data[property_name].copy()

        # Apply row-level filters (domain-based)
        if filters:
            for column, condition in filters.items():
                if column == 'multi' and callable(condition):
                    df = df[df.apply(condition, axis=1)]
                elif callable(condition):
                    df = df[condition(df[column])]
                elif isinstance(condition, tuple) and len(condition) == 2:
                    df = df[(df[column] >= condition[0]) & (df[column] <= condition[1])]
                elif isinstance(condition, list):
                    df = df[df[column].isin(condition)]
                else:
                    df = df[df[column] == condition]

        # Optionally restrict to a subset of models
        if models_to_include is not None:
            domain_keys = [col for col in ['N', 'Z'] if col in df.columns]
            allowed_cols = domain_keys + [m for m in models_to_include if m in df.columns]
            df = df[allowed_cols]

        return df

__init__(data_source=None, verbose=True)

Initialize the Dataset object.

:param data_source: Path to the data file (.h5 or .csv). :param verbose: If True, display warnings and informational messages. Default is True.

Source code in pybmc/data.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def __init__(self, data_source=None, verbose=True):
    """
    Initialize the Dataset object.

    :param data_source: Path to the data file (.h5 or .csv).
    :param verbose: If True, display warnings and informational messages. Default is True.
    """
    self.data_source = data_source
    self.data = {}  # Dictionary of model to DataFrame
    self.verbose = verbose
    self.logger = logging.getLogger(__name__)
    if not self.logger.handlers:
        handler = logging.StreamHandler()
        handler.setFormatter(logging.Formatter('%(message)s'))
        self.logger.addHandler(handler)
    self.logger.setLevel(logging.INFO if verbose else logging.WARNING)

get_subset(property_name, filters=None, models_to_include=None)

Return a filtered, wide-format DataFrame for a given property.

:param property_name: Name of the property (e.g., "BE", "ChRad"). :param filters: Dictionary of filtering rules applied to the domain columns (e.g., {"Z": (26, 28)}). :param models_to_include: Optional list of model names to retain in the output. If None, all model columns are retained. :return: Filtered wide-format DataFrame with columns: domain keys + model columns.

Source code in pybmc/data.py 
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
def get_subset(self, property_name, filters=None, models_to_include=None):
    """
    Return a filtered, wide-format DataFrame for a given property.

    :param property_name: Name of the property (e.g., "BE", "ChRad").
    :param filters: Dictionary of filtering rules applied to the domain columns (e.g., {"Z": (26, 28)}).
    :param models_to_include: Optional list of model names to retain in the output.
                            If None, all model columns are retained.
    :return: Filtered wide-format DataFrame with columns: domain keys + model columns.
    """
    if property_name not in self.data:
        raise ValueError(f"Property '{property_name}' not found in dataset.")

    df = self.data[property_name].copy()

    # Apply row-level filters (domain-based)
    if filters:
        for column, condition in filters.items():
            if column == 'multi' and callable(condition):
                df = df[df.apply(condition, axis=1)]
            elif callable(condition):
                df = df[condition(df[column])]
            elif isinstance(condition, tuple) and len(condition) == 2:
                df = df[(df[column] >= condition[0]) & (df[column] <= condition[1])]
            elif isinstance(condition, list):
                df = df[df[column].isin(condition)]
            else:
                df = df[df[column] == condition]

    # Optionally restrict to a subset of models
    if models_to_include is not None:
        domain_keys = [col for col in ['N', 'Z'] if col in df.columns]
        allowed_cols = domain_keys + [m for m in models_to_include if m in df.columns]
        df = df[allowed_cols]

    return df

load_data(models, keys=None, domain_keys=None, model_column='model', truth_column_name=None)

Load data for each property and return a dictionary of synchronized DataFrames. Each DataFrame has columns: domain_keys + one column per model for that property.

Parameters:

Name Type Description Default
models list

List of model names (for HDF5 keys or filtering CSV).

 required
keys list

List of property names to extract (each will be a key in the output dict).

 None
domain_keys list

List of columns used to define the common domain (default ['N', 'Z']).

 None
model_column str

Name of the column in the CSV that identifies which model each row belongs to. Only used for CSV files; ignored for HDF5 files.

 'model'
truth_column_name str

Name of the truth model. If provided, the truth data will be left-joined to the common domain of the other models, allowing the truth data to have a smaller domain than the models.

 None

Returns:

Name Type Description
dict

Dictionary where each key is a property name and each value is a DataFrame with columns: domain_keys + one column per model for that property. The DataFrames are synchronized to the intersection of the domains for all models. If truth_column_name is provided, truth data is left-joined (may have NaN values).

Supports both .h5 and .csv files.

Source code in pybmc/data.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 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
def load_data(self, models, keys=None, domain_keys=None, model_column='model', truth_column_name=None):
    """
    Load data for each property and return a dictionary of synchronized DataFrames.
    Each DataFrame has columns: domain_keys + one column per model for that property.

    Parameters:
        models (list): List of model names (for HDF5 keys or filtering CSV).
        keys (list): List of property names to extract (each will be a key in the output dict).
        domain_keys (list, optional): List of columns used to define the common domain (default ['N', 'Z']).
        model_column (str, optional): Name of the column in the CSV that identifies which model each row belongs to.
                                      Only used for CSV files; ignored for HDF5 files.
        truth_column_name (str, optional): Name of the truth model. If provided, the truth data will be 
                                           left-joined to the common domain of the other models, allowing 
                                           the truth data to have a smaller domain than the models.

    Returns:
        dict: Dictionary where each key is a property name and each value is a DataFrame with columns:
              domain_keys + one column per model for that property.
              The DataFrames are synchronized to the intersection of the domains for all models.
              If truth_column_name is provided, truth data is left-joined (may have NaN values).

    Supports both .h5 and .csv files.
    """
    self.domain_keys = domain_keys 

    if self.data_source is None:
        raise ValueError("Data source must be specified.")
    if not os.path.exists(self.data_source):
        raise FileNotFoundError(f"Data source '{self.data_source}' not found.")
    if keys is None:
        raise ValueError("You must specify which properties to extract via 'keys'.")

    result = {}

    for prop in keys:
        dfs = []
        truth_df = None
        skipped_models = []

        # Separate regular models from truth model
        regular_models = [m for m in models if m != truth_column_name]

        if self.data_source.endswith('.h5'):
            for model in models:
                df = pd.read_hdf(self.data_source, key=model)
                # Check required columns
                missing_cols = [col for col in domain_keys + [prop] if col not in df.columns]
                if missing_cols:
                    self.logger.info(f"[Skipped] Model '{model}' missing columns {missing_cols} for property '{prop}'.")
                    skipped_models.append(model)
                    continue
                temp = df[domain_keys + [prop]].copy()
                temp.rename(columns={prop: model}, inplace=True) # type: ignore

                # Store truth data separately if truth_column_name is provided
                if truth_column_name and model == truth_column_name:
                    truth_df = temp
                else:
                    dfs.append(temp)

        elif self.data_source.endswith('.csv'):
            df = pd.read_csv(self.data_source)
            for model in models:
                if model_column not in df.columns:
                    raise ValueError(f"Expected column '{model_column}' not found in CSV.")
                model_df = df[df[model_column] == model]
                missing_cols = [col for col in domain_keys + [prop] if col not in model_df.columns]
                if missing_cols:
                    self.logger.info(f"[Skipped] Model '{model}' missing columns {missing_cols} for key '{prop}'.")
                    skipped_models.append(model)
                    continue
                temp = model_df[domain_keys + [prop]].copy()
                temp.rename(columns={prop: model}, inplace=True)

                # Store truth data separately if truth_column_name is provided
                if truth_column_name and model == truth_column_name:
                    truth_df = temp
                else:
                    dfs.append(temp)
        else:
            raise ValueError("Unsupported file format. Only .h5 and .csv are supported.")

        if not dfs:
            self.logger.info(f"[Warning] No models with property '{prop}'. Resulting DataFrame will be empty.")
            result[prop] = pd.DataFrame(columns=domain_keys + [m for m in models if m not in skipped_models])
            continue

        # Intersect domain for regular models only
        common_df = dfs[0]
        for other_df in dfs[1:]:
            common_df = pd.merge(common_df, other_df, on=domain_keys, how="inner")
        # Drop rows with NaN in any required column (for regular models)
        common_df = common_df.dropna()

        # Left join truth data if it exists and was specified
        if truth_df is not None:
            common_df = pd.merge(common_df, truth_df, on=domain_keys, how="left")

        result[prop] = common_df
        self.data = result
    return result

separate_points_distance_allSets(list1, list2, distance1, distance2)

Separates points in list1 into three groups based on their proximity to any point in list2.

:param list1: List of (x, y) tuples. :param list2: List of (x, y) tuples. :param distance: The threshold distance to determine proximity. :return: Two lists - close_points and distant_points.

Source code in pybmc/data.py 
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
def separate_points_distance_allSets(self, list1, list2, distance1, distance2):
        """
        Separates points in list1 into three groups based on their proximity to any point in list2.

        :param list1: List of (x, y) tuples.
        :param list2: List of (x, y) tuples.
        :param distance: The threshold distance to determine proximity.
        :return: Two lists - close_points and distant_points.
        """
        train = []
        validation=[]
        test = []

        train_list_coordinates=[]
        validation_list_coordinates=[]
        test_list_coordinates=[]

        for i in range(len(list1)):
            point1=list1[i]
            close = False
            for point2 in list2:
                if np.linalg.norm(np.array(point1) - np.array(point2)) <= distance1:
                    close = True
                    break
            if close:
                train.append(point1)
                train_list_coordinates.append(i)
            else:
                close2=False
                for point2 in list2:
                    if np.linalg.norm(np.array(point1) - np.array(point2)) <= distance2:
                        close2 = True
                        break
                if close2:
                    validation.append(point1)
                    validation_list_coordinates.append(i)
                else:
                    test.append(point1)
                    test_list_coordinates.append(i)                

        return train_list_coordinates, validation_list_coordinates, test_list_coordinates

split_data(data_dict, property_name, splitting_algorithm='random', **kwargs)

Split data into training, validation, and testing sets using random or inside-to-outside logic.

:param data_dict: Dictionary output from load_data, where keys are property names and values are DataFrames. :param property_name: The key in data_dict specifying which DataFrame to use for splitting. :param splitting_algorithm: 'random' (default) or 'inside_to_outside'. :param kwargs: Additional arguments depending on the chosen algorithm. For 'random': train_size, val_size, test_size For 'inside_to_outside': stable_points (list of (x, y)), distance1, distance2 :return: Tuple of train, validation, test datasets as DataFrames.

Source code in pybmc/data.py 
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
def split_data(self, data_dict, property_name, splitting_algorithm="random", **kwargs):
    """
    Split data into training, validation, and testing sets using random or inside-to-outside logic.

    :param data_dict: Dictionary output from `load_data`, where keys are property names and values are DataFrames.
    :param property_name: The key in `data_dict` specifying which DataFrame to use for splitting.
    :param splitting_algorithm: 'random' (default) or 'inside_to_outside'.
    :param kwargs: Additional arguments depending on the chosen algorithm.
        For 'random': train_size, val_size, test_size
        For 'inside_to_outside': stable_points (list of (x, y)), distance1, distance2
    :return: Tuple of train, validation, test datasets as DataFrames.
    """
    if property_name not in data_dict:
        raise ValueError(f"Property '{property_name}' not found in the provided data dictionary.")

    data = data_dict[property_name]

    if isinstance(data, pd.DataFrame):
        indexable_data = data.reset_index(drop=True)
        point_list = list(indexable_data.itertuples(index=False, name=None))
    else:
        raise TypeError("Data for the specified property must be a pandas DataFrame.")

    if splitting_algorithm == "random":
        required = ['train_size', 'val_size', 'test_size']
        if not all(k in kwargs for k in required):
            raise ValueError(f"Missing required kwargs for 'random': {required}")

        train_size = kwargs['train_size']
        val_size = kwargs['val_size']
        test_size = kwargs['test_size']

        if not np.isclose(train_size + val_size + test_size, 1.0):
            raise ValueError("train_size + val_size + test_size must equal 1.0")

        # Random split using indexes
        train_idx, temp_idx = train_test_split(indexable_data.index, train_size=train_size, random_state=1)
        val_rel = val_size / (val_size + test_size)
        val_idx, test_idx = train_test_split(temp_idx, test_size=1 - val_rel, random_state=1)

    elif splitting_algorithm == "inside_to_outside":
        required = ['stable_points', 'distance1', 'distance2']
        if not all(k in kwargs for k in required):
            raise ValueError(f"Missing required kwargs for 'inside_to_outside': {required}")

        stable_points = kwargs['stable_points']
        distance1 = kwargs['distance1']
        distance2 = kwargs['distance2']

        train_idx, val_idx, test_idx = self.separate_points_distance_allSets(
            point_list, stable_points, distance1, distance2
        )
    else:
        raise ValueError("splitting_algorithm must be either 'random' or 'inside_to_outside'")

    train_data = indexable_data.iloc[train_idx]
    val_data = indexable_data.iloc[val_idx]
    test_data = indexable_data.iloc[test_idx]

    return train_data, val_data, test_data

view_data(property_name=None, model_name=None)

View data flexibly based on input parameters.

  • No arguments: returns available property names and model names.
  • property_name only: returns the full DataFrame for that property.
  • model_name only: Return model values across all properties.
  • property_name + model_name: returns a Series of values for the model.

:param property_name: Optional property name :param model_name: Optional model name :return: dict, DataFrame, or Series depending on input.

Source code in pybmc/data.py 
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
def view_data(self, property_name=None, model_name=None):
    """
    View data flexibly based on input parameters.

    - No arguments: returns available property names and model names.
    - property_name only: returns the full DataFrame for that property.
    - model_name only: Return model values across all properties.
    - property_name + model_name: returns a Series of values for the model.

    :param property_name: Optional property name 
    :param model_name: Optional model name 
    :return: dict, DataFrame, or Series depending on input.
    """

    if not self.data:
        raise RuntimeError("No data loaded. Run `load_data(...)` first.")

    if property_name is None and model_name is None:
        props = list(self.data.keys())
        models = sorted(set(col for prop_df in self.data.values() for col in prop_df.columns if col not in self.domain_keys))

        return {
            "available_properties": props,
            "available_models": models
        }

    if model_name is not None and property_name is None:
        # Return a dictionary: {property: Series of model values}
        result = {}
        for prop, df in self.data.items():
            if model_name in df.columns:
                cols = self.domain_keys + [model_name]
                result[prop] = df[cols]
            else:
                result[prop] = f"[Model '{model_name}' not available]"
        return result

    if property_name is not None:
        if property_name not in self.data:
            raise KeyError(f"Property '{property_name}' not found.")

        df = self.data[property_name]

        if model_name is None:
            return df  # Full property DataFrame

        if model_name not in df.columns:
            raise KeyError(f"Model '{model_name}' not found in property '{property_name}'.")

        return df[model_name]

BayesianModelCombination

The main idea of this class is to perform Bayesian Model Combination (BMC) on the set of models that we choose from the dataset class. What should this class contain: + Orthogonalization step. + Perform Bayesian inference on the training data that we extract from the Dataset class. + Predictions for certain isotopes.

Source code in pybmc/bmc.py 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
class BayesianModelCombination:
    """
    The main idea of this class is to perform Bayesian Model Combination (BMC) on the set of models that we choose
    from the dataset class. What should this class contain:
    + Orthogonalization step.
    + Perform Bayesian inference on the training data that we extract from the Dataset class.
    + Predictions for certain isotopes.
    """

    VALID_CONSTRAINTS = ("unconstrained", "simplex")
    VALID_ERROR_MODELS = tuple(VARIANCE_MODELS)

    def __init__(self, models_list, data_dict, truth_column_name, weights=None, constraint="unconstrained", error_model="homoscedastic"):
        """
        Initialize the BayesianModelCombination class.

        :param models_list: List of model names
        :param data_dict: Dictionary from `load_data()` where each key is a model name and each value is a DataFrame of properties
        :param truth_column_name: Name of the column containing the truth values.
        :param weights: Optional initial weights for the models.
        :param constraint: Weight constraint mode. Options:
            - ``"unconstrained"`` (default): No constraints on model weights.
            - ``"simplex"``: Forces weights to lie on the probability simplex
              (each weight between 0 and 1, weights sum to 1). Uses a
              Metropolis-within-Gibbs sampler to enforce the constraint.
        :param error_model: Noise structure of the combination. Options
            (see :mod:`pybmc.error_models` for the variance forms):
            - ``"homoscedastic"`` (default): A single constant variance.
            - ``"hetero_pc_dist"`` / ``"hetero_pc_dist_quad"``: Variance
              linear/quadratic in the distance from the training centroid
              in principal-component space.
            - ``"hetero_model_var"`` / ``"hetero_model_var_quad"``: Variance
              linear/quadratic in the spread among model predictions.
            - ``"hetero_combined_linear"`` / ``"hetero_combined_quadratic"``:
              Variance depending on both metrics.
            Heteroscedastic error models currently require the
            ``"unconstrained"`` weight mode.
        """

        if not isinstance(models_list, list) or not all(isinstance(model, str) for model in models_list):
            raise ValueError("The 'models' should be a list of model names (strings) for Bayesian Combination.")
        if not isinstance(data_dict, dict) or not all(isinstance(df, pd.DataFrame) for df in data_dict.values()):
            raise ValueError("The 'data_dict' should be a dictionary of pandas DataFrames, one per property.")
        if constraint not in self.VALID_CONSTRAINTS:
            raise ValueError(
                f"Invalid constraint '{constraint}'. "
                f"Must be one of {self.VALID_CONSTRAINTS}."
            )
        if error_model not in self.VALID_ERROR_MODELS:
            raise ValueError(
                f"Invalid error model '{error_model}'. "
                f"Must be one of {self.VALID_ERROR_MODELS}."
            )
        if error_model != "homoscedastic" and constraint == "simplex":
            raise ValueError(
                "Heteroscedastic error models are not supported with the "
                "'simplex' constraint; use constraint='unconstrained'."
            )

        self.data_dict = data_dict
        self.models_list = models_list
        self.models = [m for m in models_list if m != 'truth']
        self.weights = weights if weights is not None else None
        self.truth_column_name = truth_column_name
        self.constraint = constraint
        self.error_model = error_model
        self.samples = None
        self.Vt_hat = None
        self.mh_acceptance_rate_ = None
        self._trained_error_model = None
        self._metrics_calculator = None


    def orthogonalize(self, property, train_df, components_kept):
        """
        Perform orthogonalization for the specified property using training data.

        :param property: The nuclear property to orthogonalize on (e.g., 'BE').
        :param train_index: Training data from split_data
        :param components_kept: Number of SVD components to retain.
        """
        # Store selected property
        self.current_property = property

        # Extract the relevant DataFrame for that property
        df = self.data_dict[property].copy()
        self.selected_models_dataset = df  # Store for train() and predict()

        # Extract model outputs (only the model columns)
        models_output_train = train_df[self.models]
        model_predictions_train = models_output_train.values

        # Mean prediction across models (per nucleus)
        predictions_mean_train = np.mean(model_predictions_train, axis=1)

        # Experimental truth values for the property
        centered_experiment_train = train_df[self.truth_column_name].values - predictions_mean_train

        # Center model predictions
        model_predictions_train_centered = model_predictions_train - predictions_mean_train[:, None]

        # Perform SVD
        U, S, Vt = np.linalg.svd(model_predictions_train_centered)

        # Dimensionality reduction
        U_hat, S_hat, Vt_hat, Vt_hat_normalized = USVt_hat_extraction(U, S, Vt, components_kept) #type: ignore

        # Save for training
        self.centered_experiment_train = centered_experiment_train
        self.U_hat = U_hat
        self.Vt_hat = Vt_hat
        self.S_hat = S_hat
        self.Vt_hat_normalized = Vt_hat_normalized
        self._predictions_mean_train = predictions_mean_train
        # Raw training predictions, needed to fit the heteroscedasticity
        # metrics (PC-space centroid and normalization bounds).
        self._train_model_predictions = model_predictions_train


    def train(self, training_options=None):
        """
        Train the model combination using training data and optional training parameters.

        All error models (the homoscedastic one included) share the
        likelihood ``y_i ~ N(X_i . b, sigma_i^2)`` and are trained with
        the same Gibbs-within-Metropolis sampler; the homoscedastic
        model is simply the case where the variance basis is the
        constant column alone, so ``sigma_i^2 = sigma^2``.

        :param training_options: Dictionary of training options. Keys:
            - 'iterations': (int) Number of retained Gibbs samples (default 50000)
            - 'sampler': (str) Override the constraint mode for this training run.
              ``"unconstrained"`` or ``"simplex"``. If not provided, uses the
              instance-level ``self.constraint`` set at initialization.
            - 'error_model': (str) Override the error model for this training
              run (see ``VALID_ERROR_MODELS``). If not provided, uses the
              instance-level ``self.error_model`` set at initialization.
            - 'seed': (int) Seed for the sampler. If not provided, the
              sampler draws from the shared package-wide generator
              (see :mod:`pybmc.rng`), which is seeded at import.
            - 'b_mean_prior': (np.ndarray) Prior mean vector (default zeros)
              *(unconstrained sampler)*
            - 'b_mean_cov': (np.ndarray) Prior covariance matrix (default diag(S_hat²))
              *(unconstrained sampler)*
            - 'nu0_chosen': (float) Degrees of freedom for variance prior (default 1.0)
              *(simplex sampler only)*
            - 'sigma20_chosen': (float) Prior variance (default 0.02)
              *(simplex sampler only)*
            - 'burn': (int) Burn-in iterations (default 10000 for simplex,
              5000 otherwise)
            - 'stepsize': (float) Proposal step size (default 0.001)
              *(simplex sampler only)*
            - 'proposal_scales': (list) Diagonal of the Metropolis-Hastings
              proposal covariance for the variance parameters; sensible
              per-model defaults are used if omitted
              *(unconstrained sampler)*
            - 'init_params': (list) Initial values for the non-constant
              variance parameters *(unconstrained sampler)*
            - 'prior_spec': (list of (shape, scale)) Gamma priors for the
              variance parameters *(unconstrained sampler)*
            - 'adapt_proposal': (bool) Rescale the proposal during burn-in
              toward 'target_acceptance' (default True)
              *(unconstrained sampler)*
            - 'target_acceptance': (float) Acceptance rate targeted by the
              burn-in adaptation (default 0.25)
              *(unconstrained sampler)*

        After training with the unconstrained sampler, the
        Metropolis-Hastings acceptance rate of the variance parameters
        is available in ``self.mh_acceptance_rate_``.
        """
        if training_options is None:
            training_options = {}

        # Determine which sampler to use: training_options overrides instance default
        sampler_mode = training_options.get('sampler', self.constraint)
        if sampler_mode not in self.VALID_CONSTRAINTS:
            raise ValueError(
                f"Invalid sampler '{sampler_mode}'. "
                f"Must be one of {self.VALID_CONSTRAINTS}."
            )

        # Same override pattern for the error model.
        error_model_mode = training_options.get('error_model', self.error_model)
        if error_model_mode not in self.VALID_ERROR_MODELS:
            raise ValueError(
                f"Invalid error model '{error_model_mode}'. "
                f"Must be one of {self.VALID_ERROR_MODELS}."
            )
        if error_model_mode != "homoscedastic" and sampler_mode == "simplex":
            raise ValueError(
                "Heteroscedastic error models are not supported with the "
                "'simplex' sampler; use the unconstrained sampler."
            )

        iterations = training_options.get('iterations', 50000)
        num_components = self.U_hat.shape[1]
        S_hat = self.S_hat
        seed = training_options.get('seed')

        if sampler_mode == "simplex":
            nu0_chosen = training_options.get('nu0_chosen', 1.0)
            sigma20_chosen = training_options.get('sigma20_chosen', 0.02)
            burn = training_options.get('burn', 10000)
            stepsize = training_options.get('stepsize', 0.001)
            self._metrics_calculator = None
            self.mh_acceptance_rate_ = None
            self.samples = gibbs_sampler_simplex(
                self.centered_experiment_train,
                self.U_hat,
                self.Vt_hat,
                self.S_hat,
                iterations,
                [nu0_chosen, sigma20_chosen],
                burn=burn,
                stepsize=stepsize,
                seed=seed,
            )
        else:
            settings = DEFAULT_SAMPLER_SETTINGS[error_model_mode]
            terms = VARIANCE_MODELS[error_model_mode]
            metric_names = required_metrics(error_model_mode)

            if metric_names:
                # Fit the metrics (PC centroid, normalization bounds) on the
                # training predictions, then build the variance design matrix.
                self._metrics_calculator = HeteroscedasticMetrics(
                    metric_names
                ).fit(self._train_model_predictions, self.Vt_hat)
                metrics_train = self._metrics_calculator.compute(
                    self._train_model_predictions
                )
                basis_train = variance_basis(metrics_train, terms)
            else:
                # Homoscedastic: the variance basis is the constant column
                # alone, so theta = [sigma^2].
                self._metrics_calculator = None
                basis_train = np.ones(
                    (len(self.centered_experiment_train), 1)
                )

            b_mean_prior = training_options.get('b_mean_prior', np.zeros(num_components))
            b_mean_cov = training_options.get('b_mean_cov', np.diag(S_hat**2))
            self.samples, self.mh_acceptance_rate_ = gibbs_sampler_heteroscedastic(
                self.centered_experiment_train,
                self.U_hat,
                basis_train,
                iterations,
                burn=training_options.get('burn', 5000),
                proposal_scales=training_options.get(
                    'proposal_scales', settings['proposal_scales']
                ),
                init_params=training_options.get(
                    'init_params', settings['init_params']
                ),
                prior_spec=training_options.get(
                    'prior_spec', settings['prior_spec']
                ),
                b_mean_prior=b_mean_prior,
                b_mean_cov=b_mean_cov,
                adapt_proposal=training_options.get('adapt_proposal', True),
                target_acceptance=training_options.get('target_acceptance', 0.25),
                seed=seed,
            )

        # Remember which error model produced self.samples so that
        # predict()/evaluate() use the matching predictive distribution.
        self._trained_error_model = error_model_mode



    def predict(self, property, seed=DEFAULT_PREDICTIVE_SEED):
        """
        Predict a specified property using the model weights learned during training.

        :param property: The property name to predict (e.g., 'ChRad').
        :param seed: Seed for the posterior predictive draws (subsampling
            of the posterior samples and the noise added on top).
            Defaults to a fixed constant so repeated calls are
            reproducible; pass a different value for independent draws.
        :return:
            - rndm_m: array of shape (n_samples, n_points), full posterior draws
            - lower_df: DataFrame with columns domain_keys + ['Predicted_Lower']
            - median_df: DataFrame with columns domain_keys + ['Predicted_Median']
            - upper_df: DataFrame with columns domain_keys + ['Predicted_Upper']
        """
        if self.samples is None or self.Vt_hat is None:
            raise ValueError("Must call `orthogonalize()` and `train()` before predicting.")

        if property not in self.data_dict:
            raise KeyError(f"Property '{property}' not found in data_dict.")

        df = self.data_dict[property].copy()

        # Infer domain and model columns
        full_model_cols = self.models
        domain_keys = [col for col in df.columns if col not in full_model_cols and col != self.truth_column_name]

        # Determine which models are present
        available_models = [m for m in full_model_cols if m in df.columns]

        if len(available_models) == 0:
            raise ValueError("No available trained models are present in prediction DataFrame.")

        # Filter predictions and model weights
        model_preds = df[available_models].values
        domain_df = df[domain_keys].reset_index(drop=True)

        rndm_m, (lower, median, upper) = self._posterior_predictive(
            model_preds, seed=seed
        )

        # Build output DataFrames
        lower_df = domain_df.copy()

        lower_df["Predicted_Lower"] = lower

        median_df = domain_df.copy()
        median_df["Predicted_Median"] = median

        upper_df = domain_df.copy()
        upper_df["Predicted_Upper"] = upper

        return rndm_m, lower_df, median_df, upper_df

    def evaluate(self, domain_filter=None, seed=DEFAULT_PREDICTIVE_SEED):
        """
        Evaluate the model combination using coverage calculation.

        :param domain_filter: dict with optional domain key ranges, e.g., {"Z": (20, 30), "N": (20, 40)}
        :param seed: Seed for the posterior predictive draws underlying
            the coverage calculation. Defaults to a fixed constant so
            repeated calls are reproducible.
        :return: coverage list for each percentile
        """
        df = self.data_dict[self.current_property]

        if domain_filter:
            # Inline optimized filtering
            for col, cond in domain_filter.items():
                if col == 'multi' and callable(cond):
                    df = df[df.apply(cond, axis=1)]
                elif callable(cond):
                    df = df[cond(df[col])]
                elif isinstance(cond, tuple) and len(cond) == 2:
                    df = df[df[col].between(*cond)]
                elif isinstance(cond, list):
                    df = df[df[col].isin(cond)]
                else:
                    df = df[df[col] == cond]

        # Coverage is only defined where truth values exist.
        df = df.dropna(subset=[self.truth_column_name])

        preds = df[self.models].to_numpy()
        rndm_m, (lower, median, upper) = self._posterior_predictive(preds, seed=seed)

        return coverage(np.arange(0, 101, 5), rndm_m, df, truth_column=self.truth_column_name)

    def _posterior_predictive(self, model_preds, seed=DEFAULT_PREDICTIVE_SEED):
        """
        Posterior predictive draws for the given model predictions, using
        the predictive distribution matching the trained error model.

        :param model_preds: Array of shape (n_points, n_models) with one
            column per model in ``self.models`` order.
        :param seed: Seed for the posterior predictive draws.
        :return: Tuple ``(rndm_m, (lower, median, upper))``.
        """
        error_model = self._trained_error_model or "homoscedastic"
        terms = VARIANCE_MODELS[error_model]
        if terms:
            metrics = self._metrics_calculator.compute(model_preds)
            basis = variance_basis(metrics, terms)
        else:
            # Homoscedastic: constant-only variance basis.
            basis = np.ones((model_preds.shape[0], 1))
        return rndm_m_heteroscedastic_calculator(
            model_preds, self.samples, self.Vt_hat, basis, seed=seed
        )

    def get_weights(self, summary=True):
        """
        Compute model weights from posterior samples.

        Converts the sampled coefficient vectors (beta) into model weights
        using the transformation ``omega = beta @ Vt_hat + 1/M``, where M is
        the number of models.  In simplex-constrained mode, all weights are
        guaranteed to be non-negative and sum to 1.

        :param summary: If True (default), return a dictionary with
            ``'mean'``, ``'std'``, ``'median'`` arrays keyed by statistic.
            If False, return the full ``(n_samples, n_models)`` weight matrix.
        :return: Weight summary dict or full weight matrix.
        :raises ValueError: If ``train()`` has not been called.
        """
        if self.samples is None or self.Vt_hat is None:
            raise ValueError("Must call `orthogonalize()` and `train()` before getting weights.")

        # The first k columns are the PC coefficients; the remaining
        # columns are variance parameters (a single sigma^2 for the
        # homoscedastic model, several for heteroscedastic models).
        betas = self.samples[:, : self.Vt_hat.shape[0]]
        n_models = self.Vt_hat.shape[1]
        default_weights = np.full(n_models, 1.0 / n_models)
        weight_matrix = betas @ self.Vt_hat + default_weights

        if summary:
            return {
                "mean": np.mean(weight_matrix, axis=0),
                "std": np.std(weight_matrix, axis=0),
                "median": np.median(weight_matrix, axis=0),
                "models": self.models,
            }
        return weight_matrix

__init__(models_list, data_dict, truth_column_name, weights=None, constraint='unconstrained', error_model='homoscedastic')

Initialize the BayesianModelCombination class.

:param models_list: List of model names :param data_dict: Dictionary from load_data() where each key is a model name and each value is a DataFrame of properties :param truth_column_name: Name of the column containing the truth values. :param weights: Optional initial weights for the models. :param constraint: Weight constraint mode. Options: - "unconstrained" (default): No constraints on model weights. - "simplex": Forces weights to lie on the probability simplex (each weight between 0 and 1, weights sum to 1). Uses a Metropolis-within-Gibbs sampler to enforce the constraint. :param error_model: Noise structure of the combination. Options (see :mod:pybmc.error_models for the variance forms): - "homoscedastic" (default): A single constant variance. - "hetero_pc_dist" / "hetero_pc_dist_quad": Variance linear/quadratic in the distance from the training centroid in principal-component space. - "hetero_model_var" / "hetero_model_var_quad": Variance linear/quadratic in the spread among model predictions. - "hetero_combined_linear" / "hetero_combined_quadratic": Variance depending on both metrics. Heteroscedastic error models currently require the "unconstrained" weight mode.

Source code in pybmc/bmc.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
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
def __init__(self, models_list, data_dict, truth_column_name, weights=None, constraint="unconstrained", error_model="homoscedastic"):
    """
    Initialize the BayesianModelCombination class.

    :param models_list: List of model names
    :param data_dict: Dictionary from `load_data()` where each key is a model name and each value is a DataFrame of properties
    :param truth_column_name: Name of the column containing the truth values.
    :param weights: Optional initial weights for the models.
    :param constraint: Weight constraint mode. Options:
        - ``"unconstrained"`` (default): No constraints on model weights.
        - ``"simplex"``: Forces weights to lie on the probability simplex
          (each weight between 0 and 1, weights sum to 1). Uses a
          Metropolis-within-Gibbs sampler to enforce the constraint.
    :param error_model: Noise structure of the combination. Options
        (see :mod:`pybmc.error_models` for the variance forms):
        - ``"homoscedastic"`` (default): A single constant variance.
        - ``"hetero_pc_dist"`` / ``"hetero_pc_dist_quad"``: Variance
          linear/quadratic in the distance from the training centroid
          in principal-component space.
        - ``"hetero_model_var"`` / ``"hetero_model_var_quad"``: Variance
          linear/quadratic in the spread among model predictions.
        - ``"hetero_combined_linear"`` / ``"hetero_combined_quadratic"``:
          Variance depending on both metrics.
        Heteroscedastic error models currently require the
        ``"unconstrained"`` weight mode.
    """

    if not isinstance(models_list, list) or not all(isinstance(model, str) for model in models_list):
        raise ValueError("The 'models' should be a list of model names (strings) for Bayesian Combination.")
    if not isinstance(data_dict, dict) or not all(isinstance(df, pd.DataFrame) for df in data_dict.values()):
        raise ValueError("The 'data_dict' should be a dictionary of pandas DataFrames, one per property.")
    if constraint not in self.VALID_CONSTRAINTS:
        raise ValueError(
            f"Invalid constraint '{constraint}'. "
            f"Must be one of {self.VALID_CONSTRAINTS}."
        )
    if error_model not in self.VALID_ERROR_MODELS:
        raise ValueError(
            f"Invalid error model '{error_model}'. "
            f"Must be one of {self.VALID_ERROR_MODELS}."
        )
    if error_model != "homoscedastic" and constraint == "simplex":
        raise ValueError(
            "Heteroscedastic error models are not supported with the "
            "'simplex' constraint; use constraint='unconstrained'."
        )

    self.data_dict = data_dict
    self.models_list = models_list
    self.models = [m for m in models_list if m != 'truth']
    self.weights = weights if weights is not None else None
    self.truth_column_name = truth_column_name
    self.constraint = constraint
    self.error_model = error_model
    self.samples = None
    self.Vt_hat = None
    self.mh_acceptance_rate_ = None
    self._trained_error_model = None
    self._metrics_calculator = None

evaluate(domain_filter=None, seed=DEFAULT_PREDICTIVE_SEED)

Evaluate the model combination using coverage calculation.

:param domain_filter: dict with optional domain key ranges, e.g., {"Z": (20, 30), "N": (20, 40)} :param seed: Seed for the posterior predictive draws underlying the coverage calculation. Defaults to a fixed constant so repeated calls are reproducible. :return: coverage list for each percentile

Source code in pybmc/bmc.py 
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
def evaluate(self, domain_filter=None, seed=DEFAULT_PREDICTIVE_SEED):
    """
    Evaluate the model combination using coverage calculation.

    :param domain_filter: dict with optional domain key ranges, e.g., {"Z": (20, 30), "N": (20, 40)}
    :param seed: Seed for the posterior predictive draws underlying
        the coverage calculation. Defaults to a fixed constant so
        repeated calls are reproducible.
    :return: coverage list for each percentile
    """
    df = self.data_dict[self.current_property]

    if domain_filter:
        # Inline optimized filtering
        for col, cond in domain_filter.items():
            if col == 'multi' and callable(cond):
                df = df[df.apply(cond, axis=1)]
            elif callable(cond):
                df = df[cond(df[col])]
            elif isinstance(cond, tuple) and len(cond) == 2:
                df = df[df[col].between(*cond)]
            elif isinstance(cond, list):
                df = df[df[col].isin(cond)]
            else:
                df = df[df[col] == cond]

    # Coverage is only defined where truth values exist.
    df = df.dropna(subset=[self.truth_column_name])

    preds = df[self.models].to_numpy()
    rndm_m, (lower, median, upper) = self._posterior_predictive(preds, seed=seed)

    return coverage(np.arange(0, 101, 5), rndm_m, df, truth_column=self.truth_column_name)

get_weights(summary=True)

Compute model weights from posterior samples.

Converts the sampled coefficient vectors (beta) into model weights using the transformation omega = beta @ Vt_hat + 1/M, where M is the number of models. In simplex-constrained mode, all weights are guaranteed to be non-negative and sum to 1.

:param summary: If True (default), return a dictionary with 'mean', 'std', 'median' arrays keyed by statistic. If False, return the full (n_samples, n_models) weight matrix. :return: Weight summary dict or full weight matrix. :raises ValueError: If train() has not been called.

Source code in pybmc/bmc.py 
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
def get_weights(self, summary=True):
    """
    Compute model weights from posterior samples.

    Converts the sampled coefficient vectors (beta) into model weights
    using the transformation ``omega = beta @ Vt_hat + 1/M``, where M is
    the number of models.  In simplex-constrained mode, all weights are
    guaranteed to be non-negative and sum to 1.

    :param summary: If True (default), return a dictionary with
        ``'mean'``, ``'std'``, ``'median'`` arrays keyed by statistic.
        If False, return the full ``(n_samples, n_models)`` weight matrix.
    :return: Weight summary dict or full weight matrix.
    :raises ValueError: If ``train()`` has not been called.
    """
    if self.samples is None or self.Vt_hat is None:
        raise ValueError("Must call `orthogonalize()` and `train()` before getting weights.")

    # The first k columns are the PC coefficients; the remaining
    # columns are variance parameters (a single sigma^2 for the
    # homoscedastic model, several for heteroscedastic models).
    betas = self.samples[:, : self.Vt_hat.shape[0]]
    n_models = self.Vt_hat.shape[1]
    default_weights = np.full(n_models, 1.0 / n_models)
    weight_matrix = betas @ self.Vt_hat + default_weights

    if summary:
        return {
            "mean": np.mean(weight_matrix, axis=0),
            "std": np.std(weight_matrix, axis=0),
            "median": np.median(weight_matrix, axis=0),
            "models": self.models,
        }
    return weight_matrix

orthogonalize(property, train_df, components_kept)

Perform orthogonalization for the specified property using training data.

:param property: The nuclear property to orthogonalize on (e.g., 'BE'). :param train_index: Training data from split_data :param components_kept: Number of SVD components to retain.

Source code in pybmc/bmc.py 
 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
def orthogonalize(self, property, train_df, components_kept):
    """
    Perform orthogonalization for the specified property using training data.

    :param property: The nuclear property to orthogonalize on (e.g., 'BE').
    :param train_index: Training data from split_data
    :param components_kept: Number of SVD components to retain.
    """
    # Store selected property
    self.current_property = property

    # Extract the relevant DataFrame for that property
    df = self.data_dict[property].copy()
    self.selected_models_dataset = df  # Store for train() and predict()

    # Extract model outputs (only the model columns)
    models_output_train = train_df[self.models]
    model_predictions_train = models_output_train.values

    # Mean prediction across models (per nucleus)
    predictions_mean_train = np.mean(model_predictions_train, axis=1)

    # Experimental truth values for the property
    centered_experiment_train = train_df[self.truth_column_name].values - predictions_mean_train

    # Center model predictions
    model_predictions_train_centered = model_predictions_train - predictions_mean_train[:, None]

    # Perform SVD
    U, S, Vt = np.linalg.svd(model_predictions_train_centered)

    # Dimensionality reduction
    U_hat, S_hat, Vt_hat, Vt_hat_normalized = USVt_hat_extraction(U, S, Vt, components_kept) #type: ignore

    # Save for training
    self.centered_experiment_train = centered_experiment_train
    self.U_hat = U_hat
    self.Vt_hat = Vt_hat
    self.S_hat = S_hat
    self.Vt_hat_normalized = Vt_hat_normalized
    self._predictions_mean_train = predictions_mean_train
    # Raw training predictions, needed to fit the heteroscedasticity
    # metrics (PC-space centroid and normalization bounds).
    self._train_model_predictions = model_predictions_train

predict(property, seed=DEFAULT_PREDICTIVE_SEED)

Predict a specified property using the model weights learned during training.

:param property: The property name to predict (e.g., 'ChRad'). :param seed: Seed for the posterior predictive draws (subsampling of the posterior samples and the noise added on top). Defaults to a fixed constant so repeated calls are reproducible; pass a different value for independent draws. :return: - rndm_m: array of shape (n_samples, n_points), full posterior draws - lower_df: DataFrame with columns domain_keys + ['Predicted_Lower'] - median_df: DataFrame with columns domain_keys + ['Predicted_Median'] - upper_df: DataFrame with columns domain_keys + ['Predicted_Upper']

Source code in pybmc/bmc.py 
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
def predict(self, property, seed=DEFAULT_PREDICTIVE_SEED):
    """
    Predict a specified property using the model weights learned during training.

    :param property: The property name to predict (e.g., 'ChRad').
    :param seed: Seed for the posterior predictive draws (subsampling
        of the posterior samples and the noise added on top).
        Defaults to a fixed constant so repeated calls are
        reproducible; pass a different value for independent draws.
    :return:
        - rndm_m: array of shape (n_samples, n_points), full posterior draws
        - lower_df: DataFrame with columns domain_keys + ['Predicted_Lower']
        - median_df: DataFrame with columns domain_keys + ['Predicted_Median']
        - upper_df: DataFrame with columns domain_keys + ['Predicted_Upper']
    """
    if self.samples is None or self.Vt_hat is None:
        raise ValueError("Must call `orthogonalize()` and `train()` before predicting.")

    if property not in self.data_dict:
        raise KeyError(f"Property '{property}' not found in data_dict.")

    df = self.data_dict[property].copy()

    # Infer domain and model columns
    full_model_cols = self.models
    domain_keys = [col for col in df.columns if col not in full_model_cols and col != self.truth_column_name]

    # Determine which models are present
    available_models = [m for m in full_model_cols if m in df.columns]

    if len(available_models) == 0:
        raise ValueError("No available trained models are present in prediction DataFrame.")

    # Filter predictions and model weights
    model_preds = df[available_models].values
    domain_df = df[domain_keys].reset_index(drop=True)

    rndm_m, (lower, median, upper) = self._posterior_predictive(
        model_preds, seed=seed
    )

    # Build output DataFrames
    lower_df = domain_df.copy()

    lower_df["Predicted_Lower"] = lower

    median_df = domain_df.copy()
    median_df["Predicted_Median"] = median

    upper_df = domain_df.copy()
    upper_df["Predicted_Upper"] = upper

    return rndm_m, lower_df, median_df, upper_df

train(training_options=None)

Train the model combination using training data and optional training parameters.

All error models (the homoscedastic one included) share the likelihood y_i ~ N(X_i . b, sigma_i^2) and are trained with the same Gibbs-within-Metropolis sampler; the homoscedastic model is simply the case where the variance basis is the constant column alone, so sigma_i^2 = sigma^2.

:param training_options: Dictionary of training options. Keys: - 'iterations': (int) Number of retained Gibbs samples (default 50000) - 'sampler': (str) Override the constraint mode for this training run. "unconstrained" or "simplex". If not provided, uses the instance-level self.constraint set at initialization. - 'error_model': (str) Override the error model for this training run (see VALID_ERROR_MODELS). If not provided, uses the instance-level self.error_model set at initialization. - 'seed': (int) Seed for the sampler. If not provided, the sampler draws from the shared package-wide generator (see :mod:pybmc.rng), which is seeded at import. - 'b_mean_prior': (np.ndarray) Prior mean vector (default zeros) (unconstrained sampler) - 'b_mean_cov': (np.ndarray) Prior covariance matrix (default diag(S_hat²)) (unconstrained sampler) - 'nu0_chosen': (float) Degrees of freedom for variance prior (default 1.0) (simplex sampler only) - 'sigma20_chosen': (float) Prior variance (default 0.02) (simplex sampler only) - 'burn': (int) Burn-in iterations (default 10000 for simplex, 5000 otherwise) - 'stepsize': (float) Proposal step size (default 0.001) (simplex sampler only) - 'proposal_scales': (list) Diagonal of the Metropolis-Hastings proposal covariance for the variance parameters; sensible per-model defaults are used if omitted (unconstrained sampler) - 'init_params': (list) Initial values for the non-constant variance parameters (unconstrained sampler) - 'prior_spec': (list of (shape, scale)) Gamma priors for the variance parameters (unconstrained sampler) - 'adapt_proposal': (bool) Rescale the proposal during burn-in toward 'target_acceptance' (default True) (unconstrained sampler) - 'target_acceptance': (float) Acceptance rate targeted by the burn-in adaptation (default 0.25) (unconstrained sampler)

After training with the unconstrained sampler, the Metropolis-Hastings acceptance rate of the variance parameters is available in self.mh_acceptance_rate_.

Source code in pybmc/bmc.py 
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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def train(self, training_options=None):
    """
    Train the model combination using training data and optional training parameters.

    All error models (the homoscedastic one included) share the
    likelihood ``y_i ~ N(X_i . b, sigma_i^2)`` and are trained with
    the same Gibbs-within-Metropolis sampler; the homoscedastic
    model is simply the case where the variance basis is the
    constant column alone, so ``sigma_i^2 = sigma^2``.

    :param training_options: Dictionary of training options. Keys:
        - 'iterations': (int) Number of retained Gibbs samples (default 50000)
        - 'sampler': (str) Override the constraint mode for this training run.
          ``"unconstrained"`` or ``"simplex"``. If not provided, uses the
          instance-level ``self.constraint`` set at initialization.
        - 'error_model': (str) Override the error model for this training
          run (see ``VALID_ERROR_MODELS``). If not provided, uses the
          instance-level ``self.error_model`` set at initialization.
        - 'seed': (int) Seed for the sampler. If not provided, the
          sampler draws from the shared package-wide generator
          (see :mod:`pybmc.rng`), which is seeded at import.
        - 'b_mean_prior': (np.ndarray) Prior mean vector (default zeros)
          *(unconstrained sampler)*
        - 'b_mean_cov': (np.ndarray) Prior covariance matrix (default diag(S_hat²))
          *(unconstrained sampler)*
        - 'nu0_chosen': (float) Degrees of freedom for variance prior (default 1.0)
          *(simplex sampler only)*
        - 'sigma20_chosen': (float) Prior variance (default 0.02)
          *(simplex sampler only)*
        - 'burn': (int) Burn-in iterations (default 10000 for simplex,
          5000 otherwise)
        - 'stepsize': (float) Proposal step size (default 0.001)
          *(simplex sampler only)*
        - 'proposal_scales': (list) Diagonal of the Metropolis-Hastings
          proposal covariance for the variance parameters; sensible
          per-model defaults are used if omitted
          *(unconstrained sampler)*
        - 'init_params': (list) Initial values for the non-constant
          variance parameters *(unconstrained sampler)*
        - 'prior_spec': (list of (shape, scale)) Gamma priors for the
          variance parameters *(unconstrained sampler)*
        - 'adapt_proposal': (bool) Rescale the proposal during burn-in
          toward 'target_acceptance' (default True)
          *(unconstrained sampler)*
        - 'target_acceptance': (float) Acceptance rate targeted by the
          burn-in adaptation (default 0.25)
          *(unconstrained sampler)*

    After training with the unconstrained sampler, the
    Metropolis-Hastings acceptance rate of the variance parameters
    is available in ``self.mh_acceptance_rate_``.
    """
    if training_options is None:
        training_options = {}

    # Determine which sampler to use: training_options overrides instance default
    sampler_mode = training_options.get('sampler', self.constraint)
    if sampler_mode not in self.VALID_CONSTRAINTS:
        raise ValueError(
            f"Invalid sampler '{sampler_mode}'. "
            f"Must be one of {self.VALID_CONSTRAINTS}."
        )

    # Same override pattern for the error model.
    error_model_mode = training_options.get('error_model', self.error_model)
    if error_model_mode not in self.VALID_ERROR_MODELS:
        raise ValueError(
            f"Invalid error model '{error_model_mode}'. "
            f"Must be one of {self.VALID_ERROR_MODELS}."
        )
    if error_model_mode != "homoscedastic" and sampler_mode == "simplex":
        raise ValueError(
            "Heteroscedastic error models are not supported with the "
            "'simplex' sampler; use the unconstrained sampler."
        )

    iterations = training_options.get('iterations', 50000)
    num_components = self.U_hat.shape[1]
    S_hat = self.S_hat
    seed = training_options.get('seed')

    if sampler_mode == "simplex":
        nu0_chosen = training_options.get('nu0_chosen', 1.0)
        sigma20_chosen = training_options.get('sigma20_chosen', 0.02)
        burn = training_options.get('burn', 10000)
        stepsize = training_options.get('stepsize', 0.001)
        self._metrics_calculator = None
        self.mh_acceptance_rate_ = None
        self.samples = gibbs_sampler_simplex(
            self.centered_experiment_train,
            self.U_hat,
            self.Vt_hat,
            self.S_hat,
            iterations,
            [nu0_chosen, sigma20_chosen],
            burn=burn,
            stepsize=stepsize,
            seed=seed,
        )
    else:
        settings = DEFAULT_SAMPLER_SETTINGS[error_model_mode]
        terms = VARIANCE_MODELS[error_model_mode]
        metric_names = required_metrics(error_model_mode)

        if metric_names:
            # Fit the metrics (PC centroid, normalization bounds) on the
            # training predictions, then build the variance design matrix.
            self._metrics_calculator = HeteroscedasticMetrics(
                metric_names
            ).fit(self._train_model_predictions, self.Vt_hat)
            metrics_train = self._metrics_calculator.compute(
                self._train_model_predictions
            )
            basis_train = variance_basis(metrics_train, terms)
        else:
            # Homoscedastic: the variance basis is the constant column
            # alone, so theta = [sigma^2].
            self._metrics_calculator = None
            basis_train = np.ones(
                (len(self.centered_experiment_train), 1)
            )

        b_mean_prior = training_options.get('b_mean_prior', np.zeros(num_components))
        b_mean_cov = training_options.get('b_mean_cov', np.diag(S_hat**2))
        self.samples, self.mh_acceptance_rate_ = gibbs_sampler_heteroscedastic(
            self.centered_experiment_train,
            self.U_hat,
            basis_train,
            iterations,
            burn=training_options.get('burn', 5000),
            proposal_scales=training_options.get(
                'proposal_scales', settings['proposal_scales']
            ),
            init_params=training_options.get(
                'init_params', settings['init_params']
            ),
            prior_spec=training_options.get(
                'prior_spec', settings['prior_spec']
            ),
            b_mean_prior=b_mean_prior,
            b_mean_cov=b_mean_cov,
            adapt_proposal=training_options.get('adapt_proposal', True),
            target_acceptance=training_options.get('target_acceptance', 0.25),
            seed=seed,
        )

    # Remember which error model produced self.samples so that
    # predict()/evaluate() use the matching predictive distribution.
    self._trained_error_model = error_model_mode

USVt_hat_extraction(U, S, Vt, components_kept)

Extracts reduced-dimensionality matrices from SVD results.

Parameters:

Name Type Description Default
U ndarray

Left singular vectors.

 required
S ndarray

Singular values.

 required
Vt ndarray

Right singular vectors (transposed).

 required
components_kept int

Number of components to retain.

 required

Returns:

Type Description

tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]: - U_hat (numpy.ndarray): Reduced left singular vectors. - S_hat (numpy.ndarray): Retained singular values. - Vt_hat (numpy.ndarray): Normalized right singular vectors. - Vt_hat_normalized (numpy.ndarray): Original right singular vectors.
Source code in pybmc/inference_utils.py 
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
def USVt_hat_extraction(U, S, Vt, components_kept):
    """
    Extracts reduced-dimensionality matrices from SVD results.

    Args:
        U (numpy.ndarray): Left singular vectors.
        S (numpy.ndarray): Singular values.
        Vt (numpy.ndarray): Right singular vectors (transposed).
        components_kept (int): Number of components to retain.

    Returns:
        tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]:
            - `U_hat` (numpy.ndarray): Reduced left singular vectors.
            - `S_hat` (numpy.ndarray): Retained singular values.
            - `Vt_hat` (numpy.ndarray): Normalized right singular vectors.
            - `Vt_hat_normalized` (numpy.ndarray): Original right singular vectors.
    """
    U_hat = np.array([U.T[i] for i in range(components_kept)]).T
    S_hat = S[:components_kept]
    Vt_hat = np.array([Vt[i] / S[i] for i in range(components_kept)])
    Vt_hat_normalized = np.array([Vt[i] for i in range(components_kept)])
    return U_hat, S_hat, Vt_hat, Vt_hat_normalized

gibbs_sampler(y, X, iterations, prior_info=None, burn=5000, seed=None)

Gibbs sampler for the homoscedastic (constant-variance) model.

The homoscedastic model is the constant-only special case of the heteroscedastic likelihood (sigma_i^2 = sigma^2 for every point), so this is a thin wrapper around gibbs_sampler_heteroscedastic with a constant-only variance basis. There is no separate homoscedastic likelihood implementation.

Parameters:

Name Type Description Default
y ndarray

Response vector (centered).

 required
X ndarray

Design matrix.

 required
iterations int

Number of retained posterior samples.

 required
prior_info dict

Optional priors with keys 'b_mean_prior', 'b_mean_cov' (Gaussian prior on the coefficients) and 'prior_spec' (Gamma (shape, scale) prior on sigma^2). See gibbs_sampler_heteroscedastic.

 None
burn int

Burn-in iterations (default: 5000).

 5000
seed int | Generator | None

Seed for the sampler. None (default) draws from the shared package-wide generator (see :mod:pybmc.rng).

 None

Returns:

Type Description

numpy.ndarray: Posterior samples [beta_1..beta_k, sigma^2].
Source code in pybmc/inference_utils.py 
 7
 8
 9
10
11
12
13
14
15
16
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
def gibbs_sampler(y, X, iterations, prior_info=None, burn=5000, seed=None):
    """
    Gibbs sampler for the homoscedastic (constant-variance) model.

    The homoscedastic model is the constant-only special case of the
    heteroscedastic likelihood (``sigma_i^2 = sigma^2`` for every point),
    so this is a thin wrapper around `gibbs_sampler_heteroscedastic`
    with a constant-only variance basis. There is no separate
    homoscedastic likelihood implementation.

    Args:
        y (numpy.ndarray): Response vector (centered).
        X (numpy.ndarray): Design matrix.
        iterations (int): Number of retained posterior samples.
        prior_info (dict, optional): Optional priors with keys
            ``'b_mean_prior'``, ``'b_mean_cov'`` (Gaussian prior on the
            coefficients) and ``'prior_spec'`` (Gamma ``(shape, scale)``
            prior on ``sigma^2``). See `gibbs_sampler_heteroscedastic`.
        burn (int, optional): Burn-in iterations (default: 5000).
        seed (int | numpy.random.Generator | None, optional): Seed for
            the sampler. None (default) draws from the shared
            package-wide generator (see :mod:`pybmc.rng`).

    Returns:
        numpy.ndarray: Posterior samples ``[beta_1..beta_k, sigma^2]``.
    """
    prior_info = prior_info or {}
    variance_basis = np.ones((len(y), 1))
    samples, _ = gibbs_sampler_heteroscedastic(
        y,
        X,
        variance_basis,
        iterations,
        burn=burn,
        prior_spec=prior_info.get("prior_spec"),
        b_mean_prior=prior_info.get("b_mean_prior"),
        b_mean_cov=prior_info.get("b_mean_cov"),
        seed=seed,
    )
    return samples

gibbs_sampler_heteroscedastic(y, X, variance_basis, iterations, burn=5000, proposal_scales=None, init_params=None, prior_spec=None, b_mean_prior=None, b_mean_cov=None, adapt_proposal=True, target_acceptance=0.25, seed=None)

Gibbs-within-Metropolis sampler for the shared error-model likelihood.

The regression model is y_i ~ N(X_i . b, sigma_i^2) with a per-point variance that is linear in basis functions of the heteroscedasticity metrics:

``sigma_i^2 = variance_basis[i] . theta``

The homoscedastic model is the special case where the basis is a single column of ones, so theta = [sigma^2] is the constant variance. All error models share this one likelihood.

The coefficients b are updated with a conjugate (weighted least-squares) Gibbs step; the variance parameters theta with a positivity-constrained Gaussian random-walk Metropolis-Hastings step under Gamma priors. Because theta is kept strictly positive by the MH step and variance_basis is validated to be non-negative with a leading column of ones, every per-point variance satisfies sigma_i^2 >= theta_1 > 0 by construction — no variance flooring is applied during sampling. The only floor is on the initial constant term, which is estimated from the data residuals and could otherwise be exactly zero for a perfectly fitting model.

Parameters:

Name Type Description Default
y ndarray

Centered response vector, shape (n,).

 required
X ndarray

Design matrix (principal components), shape (n, k).

 required
variance_basis ndarray

Variance design matrix phi with a leading column of ones and no negative entries, shape (n, p). See :func:pybmc.error_models.variance_basis.

 required
iterations int

Number of retained posterior samples.

 required
burn int

Burn-in iterations discarded before retention (default: 5000).

 5000
proposal_scales list[float]

Diagonal of the Gaussian random-walk proposal covariance for theta (length p). Defaults to [1e-2, 1e-3, ..., 1e-3].

 None
init_params list[float]

Initial values for the non-constant entries of theta (length p - 1, strictly positive). The constant term starts at the OLS residual variance. Defaults to 0.01 for every term.

 None
prior_spec list[tuple[float, float]]

Gamma prior (shape, scale) for each entry of theta (length p). Defaults to (2, 10) for every parameter.

 None
b_mean_prior ndarray

Prior mean for b (default zeros).

 None
b_mean_cov ndarray

Prior covariance for b (default 1e6 * I, i.e. weakly informative).

 None
adapt_proposal bool

If True (default), rescale the proposal covariance during burn-in toward target_acceptance, so the fixed defaults work across data scales. Adaptation stops at the end of burn-in, which preserves detailed balance for the retained samples.

 True
target_acceptance float

Acceptance rate targeted by the burn-in adaptation (default: 0.25).

 0.25
seed int | Generator | None

Seed for the sampler. None (default) draws from the shared package-wide generator (see :mod:pybmc.rng).

 None

Returns:

Type Description

tuple[numpy.ndarray, float]: - samples (numpy.ndarray): Posterior samples [b_1..b_k, theta_1..theta_p], shape (iterations, k + p). - acceptance_rate (float): Post-burn-in MH acceptance rate.
Source code in pybmc/inference_utils.py 
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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
def gibbs_sampler_heteroscedastic(
    y,
    X,
    variance_basis,
    iterations,
    burn=5000,
    proposal_scales=None,
    init_params=None,
    prior_spec=None,
    b_mean_prior=None,
    b_mean_cov=None,
    adapt_proposal=True,
    target_acceptance=0.25,
    seed=None,
):
    """
    Gibbs-within-Metropolis sampler for the shared error-model likelihood.

    The regression model is ``y_i ~ N(X_i . b, sigma_i^2)`` with a
    per-point variance that is linear in basis functions of the
    heteroscedasticity metrics:

        ``sigma_i^2 = variance_basis[i] . theta``

    The homoscedastic model is the special case where the basis is a
    single column of ones, so ``theta = [sigma^2]`` is the constant
    variance. All error models share this one likelihood.

    The coefficients ``b`` are updated with a conjugate (weighted
    least-squares) Gibbs step; the variance parameters ``theta`` with a
    positivity-constrained Gaussian random-walk Metropolis-Hastings step
    under Gamma priors. Because ``theta`` is kept strictly positive by
    the MH step and ``variance_basis`` is validated to be non-negative
    with a leading column of ones, every per-point variance satisfies
    ``sigma_i^2 >= theta_1 > 0`` by construction — no variance flooring
    is applied during sampling. The only floor is on the *initial*
    constant term, which is estimated from the data residuals and could
    otherwise be exactly zero for a perfectly fitting model.

    Args:
        y (numpy.ndarray): Centered response vector, shape ``(n,)``.
        X (numpy.ndarray): Design matrix (principal components), shape
            ``(n, k)``.
        variance_basis (numpy.ndarray): Variance design matrix ``phi``
            with a leading column of ones and no negative entries, shape
            ``(n, p)``. See :func:`pybmc.error_models.variance_basis`.
        iterations (int): Number of retained posterior samples.
        burn (int, optional): Burn-in iterations discarded before
            retention (default: 5000).
        proposal_scales (list[float], optional): Diagonal of the Gaussian
            random-walk proposal covariance for ``theta`` (length p).
            Defaults to ``[1e-2, 1e-3, ..., 1e-3]``.
        init_params (list[float], optional): Initial values for the
            non-constant entries of ``theta`` (length p - 1, strictly
            positive). The constant term starts at the OLS residual
            variance. Defaults to 0.01 for every term.
        prior_spec (list[tuple[float, float]], optional): Gamma prior
            ``(shape, scale)`` for each entry of ``theta`` (length p).
            Defaults to ``(2, 10)`` for every parameter.
        b_mean_prior (numpy.ndarray, optional): Prior mean for ``b``
            (default zeros).
        b_mean_cov (numpy.ndarray, optional): Prior covariance for ``b``
            (default ``1e6 * I``, i.e. weakly informative).
        adapt_proposal (bool, optional): If True (default), rescale the
            proposal covariance during burn-in toward
            ``target_acceptance``, so the fixed defaults work across
            data scales. Adaptation stops at the end of burn-in, which
            preserves detailed balance for the retained samples.
        target_acceptance (float, optional): Acceptance rate targeted by
            the burn-in adaptation (default: 0.25).
        seed (int | numpy.random.Generator | None, optional): Seed for
            the sampler. None (default) draws from the shared
            package-wide generator (see :mod:`pybmc.rng`).

    Returns:
        tuple[numpy.ndarray, float]:
            - `samples` (numpy.ndarray): Posterior samples
              ``[b_1..b_k, theta_1..theta_p]``, shape
              ``(iterations, k + p)``.
            - `acceptance_rate` (float): Post-burn-in MH acceptance rate.
    """
    if burn < 0:
        raise ValueError("Burn-in iterations must be non-negative.")

    rng = get_rng(seed)

    n_points, n_betas = X.shape
    n_params = variance_basis.shape[1]

    if not np.allclose(variance_basis[:, 0], 1.0):
        raise ValueError(
            "The first column of 'variance_basis' must be ones "
            "(the constant variance term)."
        )
    if np.any(variance_basis < 0):
        raise ValueError(
            "'variance_basis' must be non-negative so that positive "
            "variance parameters guarantee positive variances."
        )

    if proposal_scales is None:
        proposal_scales = [1e-2] + [1e-3] * (n_params - 1)
    if len(proposal_scales) != n_params:
        raise ValueError(
            f"'proposal_scales' must have length {n_params} "
            f"(got {len(proposal_scales)})."
        )
    proposal_cov = np.diag(np.asarray(proposal_scales, dtype=float))

    if prior_spec is None:
        prior_spec = [(2, 10)] * n_params
    if len(prior_spec) != n_params:
        raise ValueError(
            f"'prior_spec' must have length {n_params} "
            f"(got {len(prior_spec)})."
        )
    prior_shapes = np.array([shape for shape, _ in prior_spec], dtype=float)
    prior_scales = np.array([scale for _, scale in prior_spec], dtype=float)

    def log_likelihood(residuals, sigma2):
        return -0.5 * float(
            np.sum(np.log(2.0 * np.pi * sigma2) + residuals**2 / sigma2)
        )

    def log_prior(theta):
        # Gamma log-density without the normalization constant, which
        # cancels in the Metropolis-Hastings ratio.
        return float(
            np.sum((prior_shapes - 1.0) * np.log(theta) - theta / prior_scales)
        )

    if b_mean_prior is None:
        b_mean_prior = np.zeros(n_betas)
    if b_mean_cov is None:
        b_mean_cov = np.eye(n_betas) * 1e6
    b_mean_cov_inv = np.linalg.inv(b_mean_cov)

    b_current = np.linalg.lstsq(X, y, rcond=None)[0]
    if init_params is None:
        init_params = [0.01] * (n_params - 1)
    if len(init_params) != n_params - 1:
        raise ValueError(
            f"'init_params' must have length {n_params - 1} "
            f"(got {len(init_params)})."
        )
    if np.any(np.asarray(init_params, dtype=float) <= 0):
        raise ValueError("'init_params' must be strictly positive.")
    theta_current = np.concatenate(
        [[max(np.var(y - X @ b_current), VARIANCE_FLOOR)], init_params]
    )

    samples = []
    accept_count = 0

    # Burn-in adaptation of the proposal scale (Robbins-Monro style):
    # every `adapt_interval` iterations the covariance is rescaled toward
    # the target acceptance rate, then frozen for the sampling phase.
    proposal_scale_factor = 1.0
    adapt_interval = 100
    accept_recent = 0

    for i in range(burn + iterations):
        # --- Gibbs step: b | theta (weighted least squares) ---
        sigma2 = variance_basis @ theta_current
        Xw = X / sigma2[:, None]
        b_post_cov = np.linalg.inv(
            X.T @ Xw + b_mean_cov_inv + np.eye(n_betas) * 1e-6
        )
        b_post_mean = b_post_cov @ (
            Xw.T @ y + b_mean_cov_inv @ b_mean_prior
        )
        b_current = rng.multivariate_normal(b_post_mean, b_post_cov)

        # --- MH step: theta | b (positivity-constrained random walk) ---
        residuals = y - X @ b_current
        theta_proposed = np.atleast_1d(
            rng.multivariate_normal(
                theta_current, proposal_cov * proposal_scale_factor**2
            )
        )
        if np.all(theta_proposed > 0):
            sigma2_proposed = variance_basis @ theta_proposed
            log_ratio = (
                log_likelihood(residuals, sigma2_proposed)
                + log_prior(theta_proposed)
            ) - (
                log_likelihood(residuals, sigma2)
                + log_prior(theta_current)
            )
            if np.log(rng.uniform()) < log_ratio:
                theta_current = theta_proposed
                if i >= burn:
                    accept_count += 1
                else:
                    accept_recent += 1

        if adapt_proposal and i < burn and (i + 1) % adapt_interval == 0:
            recent_rate = accept_recent / adapt_interval
            proposal_scale_factor *= np.exp(recent_rate - target_acceptance)
            proposal_scale_factor = float(
                np.clip(proposal_scale_factor, 1e-6, 1e6)
            )
            accept_recent = 0

        if i >= burn:
            samples.append(np.concatenate([b_current, theta_current]))

    acceptance_rate = accept_count / max(iterations, 1)
    return np.array(samples), acceptance_rate

gibbs_sampler_simplex(y, X, Vt_hat, S_hat, iterations, prior_info, burn=10000, stepsize=0.001, seed=None)

Performs Gibbs sampling with simplex constraints on model weights.

The likelihood is the same homoscedastic Gaussian model used by the unconstrained sampler (y_i ~ N(X_i . b, sigma^2)); only the weight constraint differs.

Parameters:

Name Type Description Default
y ndarray

Centered response vector.

 required
X ndarray

Design matrix of principal components.

 required
Vt_hat ndarray

Normalized right singular vectors.

 required
S_hat ndarray

Singular values.

 required
iterations int

Number of sampling iterations.

 required
prior_info list[float]

[nu0, sigma20] - prior parameters for variance.

 required
burn int

Burn-in iterations (default: 10000).

 10000
stepsize float

Proposal step size (default: 0.001).

 0.001
seed int | Generator | None

Seed for the sampler. None (default) draws from the shared package-wide generator (see :mod:pybmc.rng).

 None

Returns:

Type Description

numpy.ndarray: Posterior samples [beta_1..beta_k, sigma^2].
Source code in pybmc/inference_utils.py 
 49
 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
def gibbs_sampler_simplex(
    y, X, Vt_hat, S_hat, iterations, prior_info, burn=10000, stepsize=0.001,
    seed=None,
):
    """
    Performs Gibbs sampling with simplex constraints on model weights.

    The likelihood is the same homoscedastic Gaussian model used by the
    unconstrained sampler (``y_i ~ N(X_i . b, sigma^2)``); only the
    weight constraint differs.

    Args:
        y (numpy.ndarray): Centered response vector.
        X (numpy.ndarray): Design matrix of principal components.
        Vt_hat (numpy.ndarray): Normalized right singular vectors.
        S_hat (numpy.ndarray): Singular values.
        iterations (int): Number of sampling iterations.
        prior_info (list[float]): `[nu0, sigma20]` - prior parameters for variance.
        burn (int, optional): Burn-in iterations (default: 10000).
        stepsize (float, optional): Proposal step size (default: 0.001).
        seed (int | numpy.random.Generator | None, optional): Seed for
            the sampler. None (default) draws from the shared
            package-wide generator (see :mod:`pybmc.rng`).

    Returns:
        numpy.ndarray: Posterior samples ``[beta_1..beta_k, sigma^2]``.
    """
    rng = get_rng(seed)

    bias0 = np.full(len(Vt_hat.T), 1 / len(Vt_hat.T))
    nu0, sigma20 = prior_info
    cov_matrix_step = np.diag(S_hat**2 * stepsize**2)
    n = len(y)
    b_current = np.full(len(X.T), 0)
    supermodel_current = X.dot(b_current)
    residuals_current = y - supermodel_current
    ssr_current = np.sum(residuals_current**2)
    sigma2 = ssr_current / len(residuals_current)
    samples = []
    acceptance = 0

    # Validate inputs
    if burn < 0:
        raise ValueError("Burn-in iterations must be non-negative.")
    if stepsize <= 0:
        raise ValueError("Stepsize must be positive.")

    for i in range(burn + iterations):
        b_proposed = rng.multivariate_normal(b_current, cov_matrix_step)
        omegas_proposed = np.dot(b_proposed, Vt_hat) + bias0

        # Skip proposals with negative weights
        if not np.any(omegas_proposed < 0):
            supermodel_proposed = X.dot(b_proposed)
            residuals_proposed = y - supermodel_proposed
            ssr_proposed = np.sum(residuals_proposed**2)
            # Gaussian likelihood ratio at fixed sigma^2:
            # exp(-(SSR' - SSR) / (2 sigma^2)).
            acceptance_prob = min(
                1,
                np.exp((ssr_current - ssr_proposed) / (2 * sigma2)),
            )
            if rng.uniform() < acceptance_prob:
                b_current = np.copy(b_proposed)
                ssr_current = ssr_proposed
                if i >= burn:
                    acceptance += 1

        # Sample sigma^2 from its inverse-gamma full conditional
        # (positive by construction).
        shape_post = (nu0 + n) / 2.0
        scale_post = (nu0 * sigma20 + ssr_current) / 2.0
        sigma2 = 1 / rng.gamma(shape_post, 1 / scale_post)
        if i >= burn:
            samples.append(np.append(b_current, sigma2))

    return np.array(samples)

coverage(percentiles, rndm_m, models_output, truth_column)

Calculates coverage percentages for credible intervals.

Parameters:

Name Type Description Default
percentiles list[int]

Percentiles to evaluate (e.g., [5, 10, ..., 95]).

 required
rndm_m ndarray

Posterior samples of predictions.

 required
models_output DataFrame

DataFrame containing true values.

 required
truth_column str

Name of column with true values.

 required

Returns:

Type Description

list[float]: Coverage percentages for each percentile.
Source code in pybmc/sampling_utils.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
48
49
50
def coverage(percentiles, rndm_m, models_output, truth_column):
    """
    Calculates coverage percentages for credible intervals.

    Args:
        percentiles (list[int]): Percentiles to evaluate (e.g., `[5, 10, ..., 95]`).
        rndm_m (numpy.ndarray): Posterior samples of predictions.
        models_output (pandas.DataFrame): DataFrame containing true values.
        truth_column (str): Name of column with true values.

    Returns:
        list[float]: Coverage percentages for each percentile.
    """
    #  How often the model's credible intervals actually contain the true value
    data_total = len(rndm_m.T)  # Number of data points
    M_evals = len(rndm_m)  # Number of samples
    data_true = models_output[truth_column].tolist()

    coverage_results = []

    for p in percentiles:
        count_covered = 0
        for i in range(data_total):
            # Sort model evaluations for the i-th data point
            sorted_evals = np.sort(rndm_m.T[i])
            # Find indices for lower and upper bounds of the credible interval
            lower_idx = int((0.5 - p / 200) * M_evals)
            upper_idx = int((0.5 + p / 200) * M_evals) - 1
            # Check if the true value y[i] is within this interval
            if sorted_evals[lower_idx] <= data_true[i] <= sorted_evals[upper_idx]:
                count_covered += 1
        coverage_results.append(count_covered / data_total * 100)

    return coverage_results

coverage_quality(percentiles, coverage_results)

Scalar calibration score: mean |empirical - nominal| coverage.

Parameters:

Name Type Description Default
percentiles array - like

Nominal credible-interval widths in percent (as passed to coverage).

 required
coverage_results array - like

Empirical coverage in percent.

 required

Returns:

Name Type Description
float

Mean absolute deviation in percentage points (lower is

better; 0 is perfect calibration).
Source code in pybmc/sampling_utils.py 
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
def coverage_quality(percentiles, coverage_results):
    """
    Scalar calibration score: mean |empirical - nominal| coverage.

    Args:
        percentiles (array-like): Nominal credible-interval widths in
            percent (as passed to `coverage`).
        coverage_results (array-like): Empirical coverage in percent.

    Returns:
        float: Mean absolute deviation in percentage points (lower is
        better; 0 is perfect calibration).
    """
    return float(
        np.mean(np.abs(np.asarray(coverage_results) - np.asarray(percentiles)))
    )

diagnose_coverage_shape(percentiles, coverage_results, bias_tolerance=5.0, balance_threshold=0.7)

Classifies credible intervals as under-/over-dispersed or calibrated.

Parameters:

Name Type Description Default
percentiles array - like

Nominal interval widths in percent.

 required
coverage_results array - like

Empirical coverage in percent.

 required
bias_tolerance float

Mean absolute deviation (percentage points) below which the model counts as well calibrated (default: 5.0).

 5.0
balance_threshold float

Fraction of intervals that must lie consistently on one side of nominal to call the direction (default: 0.7).

 0.7

Returns:

Name Type Description
dict

Keys 'diagnosis' ('well_calibrated' |

'underdispersed' | 'overdispersed' | 'mixed'),

'mean_bias', 'mean_abs_error', 'frac_below',

'frac_above' and 'residuals'. Underdispersed means the

intervals are too narrow (overconfident); overdispersed too wide.
Source code in pybmc/sampling_utils.py 
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
def diagnose_coverage_shape(
    percentiles, coverage_results, bias_tolerance=5.0, balance_threshold=0.7
):
    """
    Classifies credible intervals as under-/over-dispersed or calibrated.

    Args:
        percentiles (array-like): Nominal interval widths in percent.
        coverage_results (array-like): Empirical coverage in percent.
        bias_tolerance (float, optional): Mean absolute deviation
            (percentage points) below which the model counts as well
            calibrated (default: 5.0).
        balance_threshold (float, optional): Fraction of intervals that
            must lie consistently on one side of nominal to call the
            direction (default: 0.7).

    Returns:
        dict: Keys ``'diagnosis'`` (``'well_calibrated'`` |
        ``'underdispersed'`` | ``'overdispersed'`` | ``'mixed'``),
        ``'mean_bias'``, ``'mean_abs_error'``, ``'frac_below'``,
        ``'frac_above'`` and ``'residuals'``. Underdispersed means the
        intervals are too narrow (overconfident); overdispersed too wide.
    """
    percentiles = np.asarray(percentiles, dtype=float)
    coverage_results = np.asarray(coverage_results, dtype=float)

    residuals = coverage_results - percentiles
    mean_bias = float(np.mean(residuals))
    mean_abs_error = float(np.mean(np.abs(residuals)))
    frac_below = float(np.mean(residuals < 0))
    frac_above = float(np.mean(residuals > 0))

    if mean_abs_error <= bias_tolerance:
        diagnosis = "well_calibrated"
    elif mean_bias < -bias_tolerance and frac_below >= balance_threshold:
        diagnosis = "underdispersed"
    elif mean_bias > bias_tolerance and frac_above >= balance_threshold:
        diagnosis = "overdispersed"
    else:
        diagnosis = "mixed"

    return {
        "diagnosis": diagnosis,
        "mean_bias": mean_bias,
        "mean_abs_error": mean_abs_error,
        "frac_below": frac_below,
        "frac_above": frac_above,
        "residuals": residuals,
    }

mace(rndm_m, y_true, quantile_levels=None)

Mean Absolute Calibration Error (MACE): a quantile-based calibration score.

Complementary to coverage/coverage_quality's two-sided central-interval view: for each one-sided quantile level q, compares the empirical fraction of true values at or below the predictive q-th percentile against q itself.

Parameters:

Name Type Description Default
rndm_m ndarray

Posterior predictive draws, shape (n_samples, n_points).

 required
y_true array - like

True/observed values, shape (n_points,).

 required
quantile_levels array - like

Quantile levels in percent, 0-100 (default: 5, 10, ..., 95).

 None

Returns:

Name Type Description
float

Mean absolute deviation between empirical and nominal

quantile coverage, in percentage points (lower is better; 0 is

perfect calibration).
Source code in pybmc/sampling_utils.py 
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 mace(rndm_m, y_true, quantile_levels=None):
    """
    Mean Absolute Calibration Error (MACE): a quantile-based calibration score.

    Complementary to `coverage`/`coverage_quality`'s two-sided central-interval
    view: for each one-sided quantile level ``q``, compares the empirical
    fraction of true values at or below the predictive ``q``-th percentile
    against ``q`` itself.

    Args:
        rndm_m (numpy.ndarray): Posterior predictive draws, shape
            ``(n_samples, n_points)``.
        y_true (array-like): True/observed values, shape ``(n_points,)``.
        quantile_levels (array-like, optional): Quantile levels in percent,
            0-100 (default: 5, 10, ..., 95).

    Returns:
        float: Mean absolute deviation between empirical and nominal
        quantile coverage, in percentage points (lower is better; 0 is
        perfect calibration).
    """
    y_true = np.asarray(y_true, dtype=float)
    if quantile_levels is None:
        quantile_levels = np.arange(5, 100, 5)
    quantile_levels = np.asarray(quantile_levels, dtype=float)

    predicted_quantiles = np.percentile(rndm_m, quantile_levels, axis=0)
    empirical = np.mean(y_true[None, :] <= predicted_quantiles, axis=1) * 100.0
    return float(np.mean(np.abs(empirical - quantile_levels)))

reduced_chi_square(rndm_m, y_true)

Reduced chi-squared statistic of the posterior predictive distribution.

chi^2_red = mean_i[ (y_true_i - mean_i)^2 / var_i ] where mean_i/var_i are the posterior predictive mean/variance at point i. A value near 1 indicates well-calibrated predictive uncertainties; > 1 means the intervals are too narrow (overconfident, "underdispersed" in the language of diagnose_coverage_shape); < 1 means too wide ("overdispersed").

Parameters:

Name Type Description Default
rndm_m ndarray

Posterior predictive draws, shape (n_samples, n_points).

 required
y_true array - like

True/observed values, shape (n_points,).

 required

Returns:

Name Type Description
float

Reduced chi-squared statistic.

Raises:

Type Description
ValueError

If the posterior predictive variance is non-positive at any point (a degenerate/zero-noise predictive distribution).
Source code in pybmc/sampling_utils.py 
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
def reduced_chi_square(rndm_m, y_true):
    """
    Reduced chi-squared statistic of the posterior predictive distribution.

    ``chi^2_red = mean_i[ (y_true_i - mean_i)^2 / var_i ]`` where
    ``mean_i``/``var_i`` are the posterior predictive mean/variance at
    point ``i``. A value near 1 indicates well-calibrated predictive
    uncertainties; > 1 means the intervals are too narrow (overconfident,
    "underdispersed" in the language of `diagnose_coverage_shape`); < 1
    means too wide ("overdispersed").

    Args:
        rndm_m (numpy.ndarray): Posterior predictive draws, shape
            ``(n_samples, n_points)``.
        y_true (array-like): True/observed values, shape ``(n_points,)``.

    Returns:
        float: Reduced chi-squared statistic.

    Raises:
        ValueError: If the posterior predictive variance is non-positive
            at any point (a degenerate/zero-noise predictive distribution).
    """
    y_true = np.asarray(y_true, dtype=float)
    pred_mean = np.mean(rndm_m, axis=0)
    pred_var = np.var(rndm_m, axis=0)
    if np.any(pred_var <= 0):
        raise ValueError(
            "Posterior predictive variance must be positive at every point."
        )
    return float(np.mean((y_true - pred_mean) ** 2 / pred_var))

rndm_m_heteroscedastic_calculator(filtered_model_predictions, samples, Vt_hat, variance_basis, seed=DEFAULT_PREDICTIVE_SEED)

Generates posterior predictive samples for any error model.

The noise added to each prediction has a per-point variance sigma_i^2 = phi_i . theta where theta are the variance parameters of each posterior draw. The homoscedastic model is the special case of a constant-only basis (a single column of ones), handled by the rndm_m_random_calculator convenience wrapper.

Parameters:

Name Type Description Default
filtered_model_predictions ndarray

Model predictions, shape (n_points, n_models).

 required
samples ndarray

Posterior samples [beta_1..beta_k, theta_1..theta_p] from gibbs_sampler_heteroscedastic.

 required
Vt_hat ndarray

Normalized right singular vectors, shape (k, n_models).

 required
variance_basis ndarray

Variance design matrix for the prediction points, shape (n_points, p).

 required
seed int | None

Seed for the posterior predictive draws (subsampling of samples and the noise added on top). Defaults to DEFAULT_PREDICTIVE_SEED; pass None to draw from the shared package-wide stream instead.

 DEFAULT_PREDICTIVE_SEED

Returns:

Type Description

tuple[numpy.ndarray, list[numpy.ndarray]]: - rndm_m (numpy.ndarray): Posterior predictive samples. - [lower, median, upper] (list[numpy.ndarray]): 95% credible interval bounds and median.
Source code in pybmc/sampling_utils.py 
 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
def rndm_m_heteroscedastic_calculator(
    filtered_model_predictions, samples, Vt_hat, variance_basis,
    seed=DEFAULT_PREDICTIVE_SEED,
):
    """
    Generates posterior predictive samples for any error model.

    The noise added to each prediction has a per-point variance
    ``sigma_i^2 = phi_i . theta`` where ``theta`` are the variance
    parameters of each posterior draw. The homoscedastic model is the
    special case of a constant-only basis (a single column of ones),
    handled by the `rndm_m_random_calculator` convenience wrapper.

    Args:
        filtered_model_predictions (numpy.ndarray): Model predictions,
            shape ``(n_points, n_models)``.
        samples (numpy.ndarray): Posterior samples
            ``[beta_1..beta_k, theta_1..theta_p]`` from
            `gibbs_sampler_heteroscedastic`.
        Vt_hat (numpy.ndarray): Normalized right singular vectors, shape
            ``(k, n_models)``.
        variance_basis (numpy.ndarray): Variance design matrix for the
            prediction points, shape ``(n_points, p)``.
        seed (int | None, optional): Seed for the posterior predictive
            draws (subsampling of `samples` and the noise added on top).
            Defaults to `DEFAULT_PREDICTIVE_SEED`; pass None to draw
            from the shared package-wide stream instead.

    Returns:
        tuple[numpy.ndarray, list[numpy.ndarray]]:
            - `rndm_m` (numpy.ndarray): Posterior predictive samples.
            - `[lower, median, upper]` (list[numpy.ndarray]): 95% credible
              interval bounds and median.
    """
    rng = get_rng(seed)

    n_draws = min(10000, len(samples))
    replace = len(samples) < 10000
    theta_rand_selected = rng.choice(samples, n_draws, replace=replace)

    n_components = Vt_hat.shape[0]
    betas = theta_rand_selected[:, :n_components]
    variance_params = theta_rand_selected[:, n_components:]

    # Model weights and noiseless central predictions.
    default_weights = np.full(Vt_hat.shape[1], 1 / Vt_hat.shape[1])
    model_weights_random = betas @ Vt_hat + default_weights
    yvals_central = model_weights_random @ filtered_model_predictions.T

    # Per-draw, per-point variances (n_draws, n_points). The floor is
    # needed here (unlike during sampling) because the normalized metrics
    # of extrapolated points can be negative, so phi . theta can dip
    # below zero even though every theta draw is positive.
    sigma2 = variance_params @ variance_basis.T
    sigma2 = np.maximum(sigma2, VARIANCE_FLOOR)

    noise = rng.standard_normal(yvals_central.shape) * np.sqrt(sigma2)
    rndm_m = yvals_central + noise

    lower = np.percentile(rndm_m, 2.5, axis=0)
    median = np.percentile(rndm_m, 50, axis=0)
    upper = np.percentile(rndm_m, 97.5, axis=0)

    return rndm_m, [lower, median, upper]

rndm_m_random_calculator(filtered_model_predictions, samples, Vt_hat, seed=DEFAULT_PREDICTIVE_SEED)

Posterior predictive samples for the homoscedastic model.

The homoscedastic model is the constant-only special case of the heteroscedastic predictive distribution, so this delegates to rndm_m_heteroscedastic_calculator with a constant-only variance basis.

Parameters:

Name Type Description Default
filtered_model_predictions ndarray

Model predictions.

 required
samples ndarray

Posterior samples [beta, sigma^2].

 required
Vt_hat ndarray

Normalized right singular vectors.

 required
seed int | None

Seed for the posterior predictive draws (subsampling of samples and the noise added on top). Defaults to DEFAULT_PREDICTIVE_SEED; pass None to draw from the shared package-wide stream instead.

 DEFAULT_PREDICTIVE_SEED

Returns:

Type Description

tuple[numpy.ndarray, list[numpy.ndarray]]: - rndm_m (numpy.ndarray): Posterior predictive samples. - [lower, median, upper] (list[numpy.ndarray]): Credible interval arrays.
Source code in pybmc/sampling_utils.py 
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
def rndm_m_random_calculator(
    filtered_model_predictions, samples, Vt_hat, seed=DEFAULT_PREDICTIVE_SEED
):
    """
    Posterior predictive samples for the homoscedastic model.

    The homoscedastic model is the constant-only special case of the
    heteroscedastic predictive distribution, so this delegates to
    `rndm_m_heteroscedastic_calculator` with a constant-only variance
    basis.

    Args:
        filtered_model_predictions (numpy.ndarray): Model predictions.
        samples (numpy.ndarray): Posterior samples `[beta, sigma^2]`.
        Vt_hat (numpy.ndarray): Normalized right singular vectors.
        seed (int | None, optional): Seed for the posterior predictive
            draws (subsampling of `samples` and the noise added on top).
            Defaults to `DEFAULT_PREDICTIVE_SEED`; pass None to draw
            from the shared package-wide stream instead.

    Returns:
        tuple[numpy.ndarray, list[numpy.ndarray]]:
            - `rndm_m` (numpy.ndarray): Posterior predictive samples.
            - `[lower, median, upper]` (list[numpy.ndarray]): Credible interval arrays.
    """
    constant_basis = np.ones((filtered_model_predictions.shape[0], 1))
    return rndm_m_heteroscedastic_calculator(
        filtered_model_predictions, samples, Vt_hat, constant_basis, seed=seed
    )

Heteroscedastic error models for Bayesian model combination.

All error models share the mean structure of the orthogonalized BMC regression and differ only in how the noise variance depends on the data point:

============================= ==================================================== name sigma_i^2 ============================= ==================================================== homoscedastic sigma^2 (constant) hetero_pc_dist alpha + beta * d_i hetero_model_var alpha + beta * v_i hetero_pc_dist_quad alpha + beta_1 d_i + beta_2 d_i^2 hetero_model_var_quad alpha + beta_1 v_i + beta_2 v_i^2 hetero_combined_linear alpha + beta_d d_i + beta_m v_i hetero_combined_quadratic alpha + beta_d1 d_i + beta_d2 d_i^2 + beta_m1 v_i + beta_m2 v_i^2 ============================= ====================================================

with two physics-informed, per-point metrics:

  • pc_dist (d): Euclidean distance from the training-data centroid in principal-component space — grows away from the fitted region.
  • model_var (v): variance among the individual model predictions — grows where the models disagree.

Both metrics are min-max normalized using the training points only, so values on extrapolated points can exceed 1.

Every model — the homoscedastic one included — is parametrized on the variance (sigma^2) scale and shares the single likelihood y_i ~ N(X_i . b, sigma_i^2) sampled by :func:pybmc.inference_utils.gibbs_sampler_heteroscedastic; the homoscedastic model is simply the case where the variance basis is the constant column alone.

HeteroscedasticMetrics

Computes and normalizes the per-point metrics of the error models.

The object is fit on the training model predictions: it stores the training centroid in PC space and the training min/max of every metric. Metrics for any other set of points are then computed consistently and scaled with the training bounds, so extrapolated points can legitimately exceed 1.

Source code in pybmc/error_models.py 
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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
class HeteroscedasticMetrics:
    """
    Computes and normalizes the per-point metrics of the error models.

    The object is fit on the training model predictions: it stores the
    training centroid in PC space and the training min/max of every
    metric. Metrics for any other set of points are then computed
    consistently and scaled with the *training* bounds, so extrapolated
    points can legitimately exceed 1.
    """

    def __init__(self, metric_names):
        """
        :param metric_names: Metric names to compute
            (subset of ``{'pc_dist', 'model_var'}``).
        """
        unknown = set(metric_names) - {"pc_dist", "model_var"}
        if unknown:
            raise ValueError(f"Unknown metric names: {sorted(unknown)}")
        self.metric_names = list(metric_names)
        self.Vt_hat = None
        self.pc_centroid = None
        self.bounds = {}

    def fit(self, train_model_predictions, Vt_hat):
        """
        Fit the normalization on the training predictions.

        Args:
            train_model_predictions (numpy.ndarray): Training model
                outputs, shape ``(n_train, n_models)``.
            Vt_hat (numpy.ndarray): Scaled right singular vectors from
                the orthogonalization step.

        Returns:
            HeteroscedasticMetrics: ``self``, for chaining.
        """
        self.Vt_hat = np.asarray(Vt_hat)
        self.pc_centroid = np.mean(
            train_model_predictions @ self.Vt_hat.T, axis=0
        )
        train_metrics = self._raw_metrics(train_model_predictions)
        self.bounds = {}
        for name, values in train_metrics.items():
            lo, hi = float(np.min(values)), float(np.max(values))
            self.bounds[name] = (lo, hi)
        return self

    def compute(self, model_predictions):
        """
        Computes normalized metrics for a set of points.

        Args:
            model_predictions (numpy.ndarray): Model outputs, shape
                ``(n_points, n_models)``.

        Returns:
            dict[str, numpy.ndarray]: Normalized metric arrays keyed by name.
        """
        if self.Vt_hat is None:
            raise ValueError("Call `fit()` before `compute()`.")
        metrics = self._raw_metrics(model_predictions)
        for name, values in metrics.items():
            lo, hi = self.bounds[name]
            if hi > lo:
                metrics[name] = (values - lo) / (hi - lo)
            # A metric that is constant on the training set is left
            # unscaled; the sampler treats it like an extra constant.
        return metrics

    def _raw_metrics(self, model_predictions):
        metrics = {}
        if "pc_dist" in self.metric_names:
            metrics["pc_dist"] = pc_distance_metric(
                model_predictions, self.Vt_hat, self.pc_centroid
            )
        if "model_var" in self.metric_names:
            metrics["model_var"] = model_variance_metric(model_predictions)
        return metrics

__init__(metric_names)

:param metric_names: Metric names to compute (subset of {'pc_dist', 'model_var'}).

Source code in pybmc/error_models.py 
226
227
228
229
230
231
232
233
234
235
236
237
def __init__(self, metric_names):
    """
    :param metric_names: Metric names to compute
        (subset of ``{'pc_dist', 'model_var'}``).
    """
    unknown = set(metric_names) - {"pc_dist", "model_var"}
    if unknown:
        raise ValueError(f"Unknown metric names: {sorted(unknown)}")
    self.metric_names = list(metric_names)
    self.Vt_hat = None
    self.pc_centroid = None
    self.bounds = {}

compute(model_predictions)

Computes normalized metrics for a set of points.

Parameters:

Name Type Description Default
model_predictions ndarray

Model outputs, shape (n_points, n_models).

 required

Returns:

Type Description

dict[str, numpy.ndarray]: Normalized metric arrays keyed by name.
Source code in pybmc/error_models.py 
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
def compute(self, model_predictions):
    """
    Computes normalized metrics for a set of points.

    Args:
        model_predictions (numpy.ndarray): Model outputs, shape
            ``(n_points, n_models)``.

    Returns:
        dict[str, numpy.ndarray]: Normalized metric arrays keyed by name.
    """
    if self.Vt_hat is None:
        raise ValueError("Call `fit()` before `compute()`.")
    metrics = self._raw_metrics(model_predictions)
    for name, values in metrics.items():
        lo, hi = self.bounds[name]
        if hi > lo:
            metrics[name] = (values - lo) / (hi - lo)
        # A metric that is constant on the training set is left
        # unscaled; the sampler treats it like an extra constant.
    return metrics

fit(train_model_predictions, Vt_hat)

Fit the normalization on the training predictions.

Parameters:

Name Type Description Default
train_model_predictions ndarray

Training model outputs, shape (n_train, n_models).

 required
Vt_hat ndarray

Scaled right singular vectors from the orthogonalization step.

 required

Returns:

Name Type Description
HeteroscedasticMetrics

self, for chaining.
Source code in pybmc/error_models.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def fit(self, train_model_predictions, Vt_hat):
    """
    Fit the normalization on the training predictions.

    Args:
        train_model_predictions (numpy.ndarray): Training model
            outputs, shape ``(n_train, n_models)``.
        Vt_hat (numpy.ndarray): Scaled right singular vectors from
            the orthogonalization step.

    Returns:
        HeteroscedasticMetrics: ``self``, for chaining.
    """
    self.Vt_hat = np.asarray(Vt_hat)
    self.pc_centroid = np.mean(
        train_model_predictions @ self.Vt_hat.T, axis=0
    )
    train_metrics = self._raw_metrics(train_model_predictions)
    self.bounds = {}
    for name, values in train_metrics.items():
        lo, hi = float(np.min(values)), float(np.max(values))
        self.bounds[name] = (lo, hi)
    return self

model_variance_metric(model_predictions)

Variance among the model predictions for each point (NaN-aware).

Parameters:

Name Type Description Default
model_predictions ndarray

Model outputs, shape (n_points, n_models).

 required

Returns:

Type Description

numpy.ndarray: Per-point variance across models, shape (n_points,).
Source code in pybmc/error_models.py 
180
181
182
183
184
185
186
187
188
189
190
191
def model_variance_metric(model_predictions):
    """
    Variance among the model predictions for each point (NaN-aware).

    Args:
        model_predictions (numpy.ndarray): Model outputs, shape
            ``(n_points, n_models)``.

    Returns:
        numpy.ndarray: Per-point variance across models, shape ``(n_points,)``.
    """
    return np.nanvar(model_predictions, axis=1)

pc_distance_metric(model_predictions, Vt_hat, pc_centroid)

Distance of each point from a reference centroid in PC space.

The principal-component coordinates of a point are obtained by projecting its raw model predictions with Vt_hat (for a row-centered SVD this is identical to projecting the centered predictions, since the right singular vectors are orthogonal to the all-ones direction).

Parameters:

Name Type Description Default
model_predictions ndarray

Model outputs, shape (n_points, n_models).

 required
Vt_hat ndarray

Scaled right singular vectors, shape (components_kept, n_models).

 required
pc_centroid ndarray

Training centroid in PC space, shape (components_kept,).

 required

Returns:

Type Description

numpy.ndarray: Euclidean distances, shape (n_points,).
Source code in pybmc/error_models.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def pc_distance_metric(model_predictions, Vt_hat, pc_centroid):
    """
    Distance of each point from a reference centroid in PC space.

    The principal-component coordinates of a point are obtained by
    projecting its raw model predictions with ``Vt_hat`` (for a
    row-centered SVD this is identical to projecting the centered
    predictions, since the right singular vectors are orthogonal to the
    all-ones direction).

    Args:
        model_predictions (numpy.ndarray): Model outputs, shape
            ``(n_points, n_models)``.
        Vt_hat (numpy.ndarray): Scaled right singular vectors, shape
            ``(components_kept, n_models)``.
        pc_centroid (numpy.ndarray): Training centroid in PC space,
            shape ``(components_kept,)``.

    Returns:
        numpy.ndarray: Euclidean distances, shape ``(n_points,)``.
    """
    pc_coords = model_predictions @ Vt_hat.T
    return np.linalg.norm(pc_coords - pc_centroid, axis=1)

required_metrics(error_model)

Returns the metric names an error model needs.

Parameters:

Name Type Description Default
error_model str

One of the keys of VARIANCE_MODELS.

 required

Returns:

Type Description

list[str]: Unique metric names (empty for homoscedastic).
Source code in pybmc/error_models.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def required_metrics(error_model):
    """
    Returns the metric names an error model needs.

    Args:
        error_model (str): One of the keys of `VARIANCE_MODELS`.

    Returns:
        list[str]: Unique metric names (empty for ``homoscedastic``).
    """
    if error_model not in VARIANCE_MODELS:
        raise ValueError(
            f"Unknown error model '{error_model}'. "
            f"Must be one of {tuple(VARIANCE_MODELS)}."
        )
    return sorted({metric for metric, _ in VARIANCE_MODELS[error_model]})

variance_basis(metrics, terms)

Builds the variance design matrix phi for a heteroscedastic model.

The per-point variance is sigma_i^2 = phi[i] . theta where the first column of phi is ones (the constant term alpha).

Parameters:

Name Type Description Default
metrics dict[str, ndarray]

Metric arrays keyed by name.

 required
terms list[tuple[str, int]]

(metric_name, power) pairs.

 required

Returns:

Type Description

numpy.ndarray: Basis matrix, shape (n_points, 1 + len(terms)).
Source code in pybmc/error_models.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def variance_basis(metrics, terms):
    """
    Builds the variance design matrix ``phi`` for a heteroscedastic model.

    The per-point variance is ``sigma_i^2 = phi[i] . theta`` where the
    first column of ``phi`` is ones (the constant term ``alpha``).

    Args:
        metrics (dict[str, numpy.ndarray]): Metric arrays keyed by name.
        terms (list[tuple[str, int]]): ``(metric_name, power)`` pairs.

    Returns:
        numpy.ndarray: Basis matrix, shape ``(n_points, 1 + len(terms))``.
    """
    n_points = len(next(iter(metrics.values())))
    columns = [np.ones(n_points)]
    for metric, power in terms:
        columns.append(np.asarray(metrics[metric], dtype=float) ** power)
    return np.column_stack(columns)

variance_parameter_names(error_model)

Returns human-readable names of the variance parameters of a model.

Parameters:

Name Type Description Default
error_model str

One of the keys of VARIANCE_MODELS.

 required

Returns:

Type Description

list[str]: Names such as ['alpha', 'beta_pc_dist^1', ...];

['sigma^2'] for the homoscedastic model (whose constant term

is the variance itself).
Source code in pybmc/error_models.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def variance_parameter_names(error_model):
    """
    Returns human-readable names of the variance parameters of a model.

    Args:
        error_model (str): One of the keys of `VARIANCE_MODELS`.

    Returns:
        list[str]: Names such as ``['alpha', 'beta_pc_dist^1', ...]``;
        ``['sigma^2']`` for the homoscedastic model (whose constant term
        is the variance itself).
    """
    terms = VARIANCE_MODELS[error_model] if error_model in VARIANCE_MODELS else None
    if terms is None:
        raise ValueError(
            f"Unknown error model '{error_model}'. "
            f"Must be one of {tuple(VARIANCE_MODELS)}."
        )
    if not terms:
        return ["sigma^2"]
    return ["alpha"] + [f"beta_{metric}^{power}" for metric, power in terms]

Package-wide random-number generation.

All stochastic pybmc routines (the MCMC samplers in :mod:pybmc.inference_utils and the posterior-predictive draws in :mod:pybmc.sampling_utils) obtain their generator through :func:get_rng, so the whole pipeline is driven by one seeded state:

  • With seed=None a function draws from the shared package-wide generator, which is seeded with DEFAULT_SEED at import time. A fresh session that performs the same sequence of calls is therefore fully reproducible, training included.
  • With an explicit seed a function uses an independent generator seeded with that value, so a single call is reproducible in isolation regardless of what ran before it.

Use set_seed to re-seed the shared generator mid-session (e.g. at the top of a script or between repetitions of an experiment).

get_rng(seed=None)

Returns the generator to use for a stochastic routine.

Parameters:

Name Type Description Default
seed int | Generator | None

If None, the shared package-wide generator (seeded with DEFAULT_SEED at import, or the last set_seed call). If an integer, a fresh independent generator seeded with it. A ready-made numpy.random.Generator is returned unchanged.

 None

Returns:

Type Description

numpy.random.Generator: The generator to draw from.
Source code in pybmc/rng.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def get_rng(seed=None):
    """
    Returns the generator to use for a stochastic routine.

    Args:
        seed (int | numpy.random.Generator | None): If None, the shared
            package-wide generator (seeded with `DEFAULT_SEED` at import,
            or the last `set_seed` call). If an integer, a fresh
            independent generator seeded with it. A ready-made
            `numpy.random.Generator` is returned unchanged.

    Returns:
        numpy.random.Generator: The generator to draw from.
    """
    if seed is None:
        return _global_rng
    if isinstance(seed, np.random.Generator):
        return seed
    return np.random.default_rng(seed)

set_seed(seed=DEFAULT_SEED)

Re-seeds the shared package-wide generator.

Parameters:

Name Type Description Default
seed int

New seed (default: DEFAULT_SEED).

 DEFAULT_SEED

Returns:

Type Description

numpy.random.Generator: The freshly seeded shared generator.
Source code in pybmc/rng.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
def set_seed(seed=DEFAULT_SEED):
    """
    Re-seeds the shared package-wide generator.

    Args:
        seed (int, optional): New seed (default: `DEFAULT_SEED`).

    Returns:
        numpy.random.Generator: The freshly seeded shared generator.
    """
    global _global_rng
    _global_rng = np.random.default_rng(seed)
    return _global_rng