### Install Kulprit via pip

Source: https://kulprit.readthedocs.io/en/latest/index.html

Standard installation command for the Kulprit package using pip.

```bash
pip install kulprit
```

--------------------------------

### Install Kulprit from GitHub

Source: https://kulprit.readthedocs.io/en/latest/index.html

Installs the latest development version of Kulprit directly from the GitHub repository.

```bash
pip install git+https://github.com/bambinos/kulprit.git
```

--------------------------------

### Importing libraries and setting global configuration

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Initial setup including importing necessary libraries and setting the random seed for reproducibility.

```python
from arviz_plots import style
import bambi as bmb
import kulprit as kpt
import numpy as np
import pandas as pd

style.use("arviz-variat")
SEED = 9352
np.random.seed(SEED)
```

--------------------------------

### Install Kulprit with Lasso support

Source: https://kulprit.readthedocs.io/en/latest/index.html

Installs Kulprit along with the scikit-learn dependency required for L1 (Lasso) search.

```bash
pip install kulprit[lasso]
```

--------------------------------

### Install Kulprit via conda

Source: https://kulprit.readthedocs.io/en/latest/index.html

Installation command for Kulprit using the conda package manager from the conda-forge channel.

```bash
conda install -c conda-forge kulprit
```

--------------------------------

### Plot forest of submodels with kpt.plot_forest()

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Plots the posterior marginals of submodels using the plot_forest function. The example includes the reference model and submodel 13, with a figure size of (10, 5). The output shows FutureWarnings from arviz_plots about xarray join defaults.

```python
kpt.plot_forest(ppi,
                   submodels=[13],
                   include_reference=True,
                   figure_kwargs={"figsize": (10, 5)},
                  );

```

```output
/home/osvaldo/anaconda3/envs/arviz_1/lib/python3.14/site-packages/arviz_plots/plot_collection.py:56: FutureWarning: In a future version of xarray the default value for join will change from join='outer' to join='exact'. This change will result in the following ValueError: cannot be aligned with join='exact' because index/labels/sizes are not equal along these coordinates (dimensions): 'chain' ('chain',) The recommendation is to set join explicitly for this case.
  data = xr.concat(ds_list, dim="model").assign_coords(model=list(data))
/home/osvaldo/anaconda3/envs/arviz_1/lib/python3.14/site-packages/arviz_plots/plot_collection.py:56: FutureWarning: In a future version of xarray the default value for join will change from join='outer' to join='exact'. This change will result in the following ValueError: cannot be aligned with join='exact' because index/labels/sizes are not equal along these coordinates (dimensions): 'draw' ('draw',) The recommendation is to set join explicitly for this case.
  data = xr.concat(ds_list, dim="model").assign_coords(model=list(data))

```

--------------------------------

### Bambi model initialization output

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

The console output generated during the initialization of the NUTS sampler for the body fat model.

```text
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, Intercept, age, weight, height, neck, chest, abdomen, hip, thigh, knee, ankle, biceps, forearm, wrist]
```

--------------------------------

### Model preparation helper function _get_models_to_plot

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/plots/plots.html

Prepares a dictionary of models and filters variable names for plotting based on the provided submodels and reference model.

```python
def _get_models_to_plot(ppi, var_names, submodels, include_reference):
    """Prepare a dictionary of models to be plotted."""

    add_intercept = []

    if submodels is None:
        submodels = ppi
    else:
        if 0 in submodels:
            add_intercept = ["Intercept"]
        submodels = ppi[submodels]

    if not var_names:
        if include_reference:
            var_names = add_intercept + ppi._ref_terms  # pylint: disable=protected-access
        else:
            var_names = add_intercept + sorted(submodel.term_names for submodel in submodels)[-1]

    # Expand Bambi term names (e.g. "1|Subject") to the actual PyMC variable names
    # stored in the posterior (e.g. "1|Subject_offset", "1|Subject_sigma").
    expanded_names = []
    term_to_vars = ppi._term_to_vars  # pylint: disable=protected-access
    for name in var_names:
        if name in term_to_vars:
            expanded_names.extend(sorted(term_to_vars[name]))
            # For non-centered group-specific terms the term name itself corresponds to
            # the centered coefficient, which is stored in the projected idata for plotting.
            if name not in term_to_vars[name]:
                expanded_names.append(name)
        else:
            expanded_names.append(name)
    var_names = expanded_names

    if include_reference:
        models_to_plot = {"Reference": ppi.reference_model.idata.posterior}
    else:
        models_to_plot = {}

    models_to_plot.update({submodel.size: submodel.idata.posterior for submodel in submodels})

    # Keep only variable names that actually exist in the plotted models.
    available_vars = [set(dt.data_vars) for dt in models_to_plot.values()]
    if include_reference:
        keep_vars = set.intersection(*available_vars)
    else:
        keep_vars = set.union(*available_vars)
    var_names = [name for name in var_names if name in keep_vars]

    return models_to_plot, var_names
```

--------------------------------

### Perform model projection with project()

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

The project method executes the projection procedure using either a forward search or L1-regularization, or by projecting user-defined terms.

```python
def project(
    self,
    method="forward",
    user_terms=None,
    num_samples=400,
    num_clusters=20,
    early_stop=None,
    require_lower_terms=True,
    tolerance=1,
):
    """Perform model projection.

    Parameters:
    -----------
    method : str
        The search method to employ, either "forward" for a forward search, or "l1" for
        an L1-regularized search. Ignored if "user_terms" is provided.
    user_terms : list of list of str
        If a nested list of terms is provided, model with those terms will be projected
        directly.
    num_samples : int
        The number of samples to draw from the posterior predictive distribution for the
        projection procedure and ELPD computation. Defaults to 400.
    num_clusters : int
        The number of clusters to use during the forward search. Defaults to 20.
    If None, the number of clusters is set to the number of samples.
    If num_clusters is larger than num_samples, it is set to num_samples.
    early_stop : str or int, optional
        Whether to stop the search earlier. If an integer is provided, the search stops
        when the submodel size is equal to the integer. If a string is provided, the search
        stops when the difference in ELPD between the reference and submodel is small.
        There are two criteria to define what small is, "mean" and "se".
        The "mean" criterion stops the search when the difference between a the ELPD is smaller
        than 4. The "se" criterion stops the search when the ELPD of the submodel is within
        one standard error of the reference model. Defaults to None.
    require_lower_terms : bool
        Include higher-order interactions only if all lower-order interactions and main effects
        are already in the subset. Defaults to True. Ignored if user_terms is provided or
        if the method is not "forward".
    tolerance : float
        The tolerance for the optimization procedure. Defaults to 1. Decreasing this value
        will increase the accuracy of the projection at the cost of speed.
    """
    self.num_samples = num_samples
    self.num_clusters = num_clusters
    self.early_stop = early_stop
    self.tolerance = tolerance
    self._pps, self._ppc, self._weights = get_pps(
        self.reference_model.idata,
        self._response_name,
        self.num_samples,
        self.num_clusters,
        self._rng,
    )

    _check_interactions(self._ref_terms, method, require_lower_terms)

    # if user provided the terms we used them directly, no search is performed
    if user_terms is not None:
        # check if the terms are a list of lists
        if not isinstance(user_terms, list) or not all(
            isinstance(term, list) for term in user_terms
        ):
            raise ValueError("Please provide a list of lists of terms.")
        # check if the length of the submodels always increase
        if not all(
            len(user_terms[idx]) < len(user_terms[idx + 1])
            for idx in range(len(user_terms) - 1)
        ):
            raise ValueError("Please provide a list of terms in increasing order")
        # check if the listed terms are valid
        for idx, term_names in enumerate(user_terms):
            if not set(term_names).issubset(self.reference_model.term_names):
                raise ValueError(f"Term {idx} is not a valid term in the reference")

        self._list_of_submodels = user_path(self._project, user_terms, self.reference_model)
    else:
        if method not in ["forward", "l1"]:
            raise ValueError("Please select either forward search or L1 search.")

        max_terms = len(self._ref_terms)

        if isinstance(self.early_stop, int):
            if self.early_stop < 0:
                raise ValueError("The early stopping value must be a positive integer.")
```

--------------------------------

### Implement __len__ and __getitem__ for submodel access

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

These methods allow the Projector to report the number of submodels and retrieve them via index or slice.

```python
def __len__(self):
    """Return the number of submodels."""
    return len(self._list_of_submodels)

def __getitem__(self, index):
    """Return submodels by index or slice."""

    if isinstance(index, int):
        return self._list_of_submodels[index]
    else:
        n_submodels = len(self._list_of_submodels)
        if not all(-n_submodels <= i < n_submodels for i in index):
            warnings.warn(
                "At least one index is out of bounds. Ignoring out of bounds indices."
            )
            index = [i for i in index if -n_submodels <= i < n_submodels]
        return [self._list_of_submodels[i] for i in index]
```

--------------------------------

### Project parameters to submodel

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Maps active variable names to flat indices, solves the projection optimization, and constructs a SubModel object.

```python
# map active variable names to flat indices in the optimization vector
        keys = list(self._ref_var_info.keys())
        sizes = [self._ref_var_info[k][1] for k in keys]
        flat_keys = []
        for k, s in zip(keys, sizes):
            flat_keys += [k] * s
        flat_keys = np.array(flat_keys)
        active_idx = np.where(np.isin(flat_keys, list(active_var_names)))[0]
        if active_idx.size == 0:
            raise ValueError(
                "No active parameters selected for projection. "
                f"Submodel terms: {term_names}; base terms: {self._base_terms}."
            )

        new_idata, loss = solve(
            self._neg_log_likelihood,
            samples,
            self._initial_guess,
            self._ref_var_info,
            weights,
            self.tolerance,
            active_idx=active_idx,
        )

        # Add observed data and log-likelihood to the projected InferenceData object
        # We only do this for the selected projected model, not the intermediate ones
        if new_idata is not None:
            new_idata["observed_data"] = self._observed_dataset
            new_idata["log_likelihood"] = compute_llk(new_idata, self._pymc_model)

            # For non-centered group-specific terms the projected idata contains the
            # offset variables (e.g. "1|Subject_offset") while the reference idata
            # contains the deterministic coefficients (e.g. "1|Subject"). Compute and
            # store those coefficients, and keep the sigmas, so plotting can align
            # reference and projected models on the same variable names.
            vars_to_keep = set(active_var_names)
            for term_name in term_names:
                if "|" not in term_name:
                    continue
                offset_name = next(
                    (v for v in self._term_to_vars[term_name] if v.endswith("_offset")),
                    None,
                )
                if offset_name is None:
                    continue
                sigma_name = f"{term_name}_sigma"
                if sigma_name not in new_idata["posterior"].data_vars:
                    continue
                coef_name = term_name
                offset = new_idata["posterior"][offset_name]
                sigma = new_idata["posterior"][sigma_name]
                coef = offset * sigma
                # Rename the group dimension and copy its coordinates from the reference
                # model idata, so arviz-plots can align reference and projected models.
                ref_var = self.reference_model.idata.posterior[coef_name]
                old_dim = coef.dims[-1]
                new_dim = ref_var.dims[-1]
                if old_dim != new_dim:
                    coef = coef.rename({old_dim: new_dim})
                coef = coef.assign_coords({new_dim: ref_var.coords[new_dim]})
                new_idata["posterior"][coef_name] = coef
                vars_to_keep.add(sigma_name)
                vars_to_keep.add(coef_name)

            # remove the variables that are not in the submodel
            vars_to_drop = [
                var for var in new_idata["posterior"].ds.data_vars if var not in vars_to_keep
            ]
            new_idata["posterior"] = new_idata["posterior"].ds.drop_vars(vars_to_drop)

        # build SubModel object and return
        sub_model = SubModel(
            model=self._pymc_model_sw,
            idata=new_idata,
            loss=loss,
            elpd=None,
            elpd_se=None,
            elpd_dse=None,
            elpd_i=None,
            size=len(term_names),
            term_names=term_names,
            has_intercept=self._has_intercept,
        )
        return sub_model
```

--------------------------------

### Imports and module docstring for kulprit.projector

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Imports required libraries and internal modules for the projector. Includes pylint disables and the module docstring.

```python
# pylint: disable=undefined-loop-variable
# pylint: disable=too-many-instance-attributes
"""Core reference model class."""

import warnings
from copy import copy

import numpy as np
from pandas import DataFrame

from bambi import formula

from kulprit.projection.arviz_io import check_idata, compute_loo, get_observed_data, get_pps
from kulprit.projection.pymc_io import (
    add_switches,
    compile_mllk,
    compute_llk,
    turn_off_terms,
    get_model_information,
    get_term_variables,
)
from kulprit.projection.search_strategies import (
    user_path,
    forward_search,
    l1_search,
    _missing_lower_order_terms,
    _missing_fixed_effect_for_group_term,
    _get_fixed_part,
)
from kulprit.projection.solver import solve
```

--------------------------------

### Load bike rental data

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Reads the bikes.csv file into a pandas DataFrame named 'bikes'.

```python
bikes = pd.read_csv("bikes.csv")

```

--------------------------------

### Project with early stopping

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Performs projection with early stopping set to 'mean', which stops when the mean criterion is met. Then plots the comparison. The early_stop argument defaults to None and can be 'mean', 'se', or an integer for maximum model size.

```python
ppi = kpt.ProjectionPredictive(model, idata)
ppi.project(early_stop="mean")
kpt.plot_compare(ppi.compare(min_model_size=1));

```

--------------------------------

### Compare model performance

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Returns a DataFrame with performance statistics for the reference and submodels.

```python
def compare(self, stats="elpd", min_model_size=0, relative_to="reference", round_to=None):
        """Return a DataFrame with the performance stats of the reference and submodels.

        Parameters:
        -----------
        stats : str
            The statistics to compute. Defaults to "elpd".
            * "elpd": expected log (pointwise) predictive density (ELPD).
            * "mlpd": mean log predictive density (MLPD), that is, the ELPD divided by the
            number of observations.
            * "gmpd": geometric mean predictive density (GMPD), that is, exp(MLPD).
            For discrete response families the GMPD is bounded by zero and one.
        min_model_size : int
            The minimum size of the submodels to compare. Defaults to 0, which means the
            intercept-only model is included in the comparison.
        relative_to : str
            If "reference", the ELPD differences are computed relative to the reference model.
            If "full", the ELPD differences are computed relative to the full submodel.
            Defaults to "reference".
        round_to : int
            Number of decimals used to round results. Defaults to None
```

--------------------------------

### Loading data and defining the Bambi model

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Load the body fat dataset, define a linear regression model using Bambi, and compute the log-likelihood required for ELPD calculations.

```python
body = pd.read_csv("body_fat.csv")
model = bmb.Model("siri ~" +  " + ".join([c for c in body.columns if c!="siri"]), data=body)
idata = model.fit(random_seed=SEED)
model.compute_log_likelihood(idata)
```

--------------------------------

### __len__

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Returns the number of submodels in the projector.

```APIDOC
## __len__

### Description
Returns the number of submodels.

### Method
__len__

### Parameters
None

### Returns
- **int** - The number of submodels.
```

--------------------------------

### ProjectionPredictive class definition and __init__ method

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Defines the ProjectionPredictive class with its constructor. The constructor initializes attributes, sets the random number generator, validates the model is built, extracts reference model information from Bambi and PyMC, adds switches, compiles the negative log-likelihood, and computes the reference model's ELPD.

```python
[docs]
class ProjectionPredictive:
    """
    Projection Predictive class from which we perform the model selection procedure.

    Parameters:
    ----------
    model : Bambi model
        The reference model to project
    idata : InferenceData or DataTree
        The result of fitting reference model
    rng : RandomState
        Random number generator used for sampling from the posterior predictive if
        the group is not present in idata.
    """

    def __init__(self, model, idata, rng=456):
        """Builder for projection predictive model selection."""
        # initialize attributes
        self.num_samples = None
        self.num_clusters = None
        self.early_stop = None
        self.tolerance = None
        self._pps = None
        self._ppc = None
        self._weights = None
        self._list_of_submodels = []

        # set random number generator
        if rng is None:
            self._rng = np.random.default_rng()
        else:
            self._rng = np.random.default_rng(rng)

        # check we have the model fitted
        if not model.built:
            raise ValueError(
                "Before projecting, please fit the model, using the `fit` method.\n"
                "Additionally, make sure that the sampling converged, "
                "and that the model fits well the data."
            )

        # get information from Bambi's reference model
        self._has_intercept = formula.formula_has_intercept(model.formula.main)
        self._response_name = model.response_component.term.name
        component = model.components[model.family.likelihood.parent]
        self._ref_terms = [
            v.alias if v.alias is not None else k for k, v in component.common_terms.items()
        ] + [
            v.alias if v.alias is not None else k for k, v in component.group_specific_terms.items()
        ]
        self._categorical_terms = sum(term.categorical for term in component.common_terms.values())

        self._base_terms = _get_base_terms(self._has_intercept, model.constant_components)

        # get information from PyMC's reference model
        self._pymc_model = copy(model.backend.model)

        # add switches to the model to turn on/off terms in the model
        # without having to rebuild the model.
        self._pymc_model_sw, self._switches = add_switches(self._pymc_model, self._ref_terms)

        initial_point = self._pymc_model_sw.initial_point()
        self._ref_var_info = get_model_information(self._pymc_model_sw, initial_point)
        self._initial_guess = np.concatenate([np.ravel(value) for value in initial_point.values()])

        # map each switchable term to the free RVs that belong to it
        self._term_to_vars = get_term_variables(self._pymc_model_sw, self._ref_terms)

        self._neg_log_likelihood = compile_mllk(self._pymc_model_sw, initial_point)

        # get information from ArviZ's InferenceData object
        idata = check_idata(idata, model, self._rng)
        self._observed_dataset, self._observed_array = get_observed_data(idata, self._response_name)

        elpd_ref = compute_loo(idata=idata)

        self.reference_model = RefModel(
            model=model,
            idata=idata,
            elpd=elpd_ref.elpd,
            elpd_se=elpd_ref.se,
            elpd_i=elpd_ref.elpd_i,
            term_names=[fvar.name for fvar in self._pymc_model.free_RVs],
        )
```

--------------------------------

### project

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Performs model projection using either forward search or L1-regularized search, or directly projects user-specified terms.

```APIDOC
## project

### Description
Performs model projection. The method can use a forward search or an L1-regularized search, or directly project user-specified terms. The projection procedure draws samples from the posterior predictive distribution and computes ELPD for submodel selection.

### Method
project

### Parameters
- **method** (str) - Optional - The search method to employ, either "forward" for a forward search, or "l1" for an L1-regularized search. Ignored if "user_terms" is provided. Defaults to "forward".
- **user_terms** (list of list of str) - Optional - If a nested list of terms is provided, model with those terms will be projected directly.
- **num_samples** (int) - Optional - The number of samples to draw from the posterior predictive distribution for the projection procedure and ELPD computation. Defaults to 400.
- **num_clusters** (int) - Optional - The number of clusters to use during the forward search. Defaults to 20. If None, the number of clusters is set to the number of samples. If num_clusters is larger than num_samples, it is set to num_samples.
- **early_stop** (str or int, optional) - Optional - Whether to stop the search earlier. If an integer is provided, the search stops when the submodel size is equal to the integer. If a string is provided, the search stops when the difference in ELPD between the reference and submodel is small. There are two criteria to define what small is, "mean" and "se". The "mean" criterion stops the search when the difference between a the ELPD is smaller than 4. The "se" criterion stops the search when the ELPD of the submodel is within one standard error of the reference model. Defaults to None.
- **require_lower_terms** (bool) - Optional - Include higher-order interactions only if all lower-order interactions and main effects are already in the subset. Defaults to True. Ignored if user_terms is provided or if the method is not "forward".
- **tolerance** (float) - Optional - The tolerance for the optimization procedure. Defaults to 1. Decreasing this value will increase the accuracy of the projection at the cost of speed.

### Returns
None

### Raises
- **ValueError** - If user_terms is not a list of lists, if the submodel lengths do not increase, if a term is not valid in the reference, or if method is not "forward" or "l1".
- **ValueError** - If early_stop is an integer and is negative.
```

--------------------------------

### Bambi model sampling performance output

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

The console output showing the duration and iteration count for the MCMC sampling process.

```text
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.
```

--------------------------------

### Inspect ppi object

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Prints the ppi object to view submodels ordered by ELPD.

```python
ppi
```

--------------------------------

### Main plotting function implementation

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/plots/plots.html

This function handles the configuration of visual elements and statistics before calling the ArviZ PlotCollection to generate the plot.

```python
    if visuals is None:
        visuals = {}

    if stats is None:
        stats = {}

    visuals.setdefault("point_estimate_text", False)
    visuals.setdefault("credible_interval", False)
    stats.setdefault("point_estimate", {"skipna": True})

    models_to_plot, var_names = _get_models_to_plot(ppi, var_names, submodels, include_reference)

    pc = azp_plot_dist(
        models_to_plot,
        var_names=var_names,
        filter_vars=filter_vars,
        coords=coords,
        sample_dims=sample_dims,
        kind=kind,
        point_estimate=point_estimate,
        ci_kind=ci_kind,
        ci_prob=ci_prob,
        plot_collection=plot_collection,
        backend=backend,
        labeller=labeller,
        aes_by_visuals=aes_by_visuals,
        visuals=visuals,
        stats=stats,
        **pc_kwargs,
    )

    pc.add_legend("model")
    return pc
```

--------------------------------

### __getitem__

Source: https://kulprit.readthedocs.io/en/latest/_modules/kulprit/projector.html

Returns submodels by index or slice. Supports integer indexing and list of indices with bounds checking.

```APIDOC
## __getitem__

### Description
Returns submodels by index or slice. If an integer is provided, returns the submodel at that index. If a list of indices is provided, returns a list of submodels, ignoring out-of-bounds indices with a warning.

### Method
__getitem__

### Parameters
- **index** (int or list of int) - Required - Index or list of indices of submodels to retrieve.

### Returns
- **Submodel** or **list of Submodel** - The requested submodel(s).
```

--------------------------------

### Fit Bambi model with splines

Source: https://kulprit.readthedocs.io/en/latest/examples/Examples.html

Defines knots and priors, then fits a negative binomial Bambi model using a spline on 'hour' and all other predictors. Sets an alias for the spline term and fits the model, producing MCMC sampling output.

```python
knots = range(1, 22, 5)
priors = {
    "alpha": bmb.Prior("Exponential", lam=1)
}

model = bmb.Model("count ~ bs(hour, knots=knots, intercept=True) + " + " + ".join(bikes.columns.drop(["count", "hour"])),
                  bikes,
                  family="negativebinomial",
                  priors=priors)

model.set_alias({"bs(hour, knots=knots, intercept=True)": "bs_hour"})
idata = model.fit()

```

```text
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [alpha, Intercept, bs_hour, season, month, holiday, weekday, workingday, weather, temperature, humidity, windspeed]

```

```text


```

```text
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 262 seconds.

```