API Reference

This is the code documentation for the humlpy package, which implements the methodology from:

Göbler, K., Drton, M., Mukherjee, S. and Miloschewski, A. (2024). High-dimensional undirected graphical models for arbitrary mixed data. Electronic Journal of Statistics, 18(1), 2339–2404. https://doi.org/10.1214/24-EJS2254

Main Classes

Estimate a sample correlation matrix for mixed continuous/ordinal data.

Pairwise correlations are estimated as follows:

continuous - continuous: Spearman sin-transform :math:\hat{\sigma} = 2 \sin(\pi/6 \cdot \hat{\rho}_S).
continuous - ordinal: ad-hoc polyserial correlation via :class:~humlpy.correlation.PolyserialCorrelation.
ordinal - ordinal: maximum-likelihood polychoric correlation via :class:~humlpy.correlation.PolychoricCorrelation.

After calling :meth:fit, the estimated matrix is available as :attr:correlation_matrix_.

Parameters:

Name	Type	Description	Default
`n_levels_threshold`	`int`	Columns with strictly fewer unique values are treated as ordinal. Defaults to 20.	`20`

Example::

sc = SampleCorrelation().fit(X)
print(sc.correlation_matrix_)

Source code in humlpy/estimation.py

class SampleCorrelation:
    r"""Estimate a sample correlation matrix for mixed continuous/ordinal data.

    Pairwise correlations are estimated as follows:

    * **continuous - continuous**: Spearman sin-transform
      :math:`\hat{\sigma} = 2 \sin(\pi/6 \cdot \hat{\rho}_S)`.
    * **continuous - ordinal**: ad-hoc polyserial correlation via
      :class:`~humlpy.correlation.PolyserialCorrelation`.
    * **ordinal - ordinal**: maximum-likelihood polychoric correlation via
      :class:`~humlpy.correlation.PolychoricCorrelation`.

    After calling :meth:`fit`, the estimated matrix is available as
    :attr:`correlation_matrix_`.

    Args:
        n_levels_threshold: Columns with strictly fewer unique values are
            treated as ordinal.  Defaults to 20.

    Example::

        sc = SampleCorrelation().fit(X)
        print(sc.correlation_matrix_)
    """

    def __init__(
        self,
        n_levels_threshold: int = 20,
    ) -> None:
        """Construct a SampleCorrelation estimator.

        Args:
            n_levels_threshold: Unique-value threshold for ordinal detection.
                Defaults to 20.
        """
        self.n_levels_threshold = n_levels_threshold

        # set after fit
        self.correlation_matrix_: pd.DataFrame | None = None
        self.feature_names_: list[str] | None = None

    # ------------------------------------------------------------------
    # Public API
    # ------------------------------------------------------------------

    def fit(self, x: pd.DataFrame | NDArray[Any]) -> SampleCorrelation:
        """Estimate the correlation matrix from *x*.

        Args:
            x: Data matrix of shape ``(n_samples, n_features)``.
                A :class:`pandas.DataFrame` is preferred so that column names
                are preserved; a 2-D :class:`numpy.ndarray` is also accepted.

        Returns:
            *self* -- enables method chaining.

        Raises:
            ValueError: If *x* has fewer than 2 columns.
        """
        df = x if isinstance(x, pd.DataFrame) else pd.DataFrame(x)
        if df.shape[1] < 2:
            raise ValueError("X must have at least 2 columns.")

        self.feature_names_ = list(df.columns.astype(str))

        rho = self._estimate(df)

        _check_correlation_bounds(rho)
        self.correlation_matrix_ = pd.DataFrame(rho, index=self.feature_names_, columns=self.feature_names_)
        return self

    # ------------------------------------------------------------------
    # Internals
    # ------------------------------------------------------------------

    def _is_discrete(self, df: pd.DataFrame) -> NDArray[np.bool_]:
        return np.array([df.iloc[:, j].nunique() < self.n_levels_threshold for j in range(df.shape[1])])

    def _estimate(self, df: pd.DataFrame) -> NDArray[Any]:
        r"""Estimate pairwise correlations under the latent Gaussian copula.

        Rules per pair type:

        * **cont - cont**: :math:`2 \sin(\pi/6 \cdot \hat{\rho}_S)` where
            :math:`\hat{\rho}_S` is Spearman's rank correlation.
        * **cont - ord**: ad-hoc polyserial via :class:`PolyserialCorrelation`.
        * **ord - ord**: MLE polychoric via :class:`PolychoricCorrelation`.
        """
        d = df.shape[1]
        rho = np.eye(d)
        disc = self._is_discrete(df)

        for i in range(d - 1):
            for j in range(i + 1, d):
                xi = df.iloc[:, i].to_numpy(dtype=float)
                xj = df.iloc[:, j].to_numpy(dtype=float)

                if not disc[i] and not disc[j]:
                    rho_s = float(stats.spearmanr(xi, xj)[0])  # pyright: ignore[reportArgumentType]
                    r = 2.0 * np.sin(np.pi / 6.0 * rho_s)
                elif disc[i] and disc[j]:
                    r = PolychoricCorrelation().fit(xi, xj).correlation
                else:
                    r = PolyserialCorrelation(n_levels_threshold=self.n_levels_threshold).fit(xi, xj).correlation

                rho[i, j] = rho[j, i] = r

        return rho

`init(n_levels_threshold=20)`

Construct a SampleCorrelation estimator.

Parameters:

Name	Type	Description	Default
`n_levels_threshold`	`int`	Unique-value threshold for ordinal detection. Defaults to 20.	`20`

Source code in humlpy/estimation.py

def __init__(
    self,
    n_levels_threshold: int = 20,
) -> None:
    """Construct a SampleCorrelation estimator.

    Args:
        n_levels_threshold: Unique-value threshold for ordinal detection.
            Defaults to 20.
    """
    self.n_levels_threshold = n_levels_threshold

    # set after fit
    self.correlation_matrix_: pd.DataFrame | None = None
    self.feature_names_: list[str] | None = None

`fit(x)`

Estimate the correlation matrix from x.

Parameters:

Name	Type	Description	Default
`x`	`DataFrame \| NDArray[Any]`	Data matrix of shape `(n_samples, n_features)`. A :class:`pandas.DataFrame` is preferred so that column names are preserved; a 2-D :class:`numpy.ndarray` is also accepted.	required

Returns:

Type	Description
`SampleCorrelation`	self -- enables method chaining.

Raises:

Type	Description
`ValueError`	If x has fewer than 2 columns.

Source code in humlpy/estimation.py

def fit(self, x: pd.DataFrame | NDArray[Any]) -> SampleCorrelation:
    """Estimate the correlation matrix from *x*.

    Args:
        x: Data matrix of shape ``(n_samples, n_features)``.
            A :class:`pandas.DataFrame` is preferred so that column names
            are preserved; a 2-D :class:`numpy.ndarray` is also accepted.

    Returns:
        *self* -- enables method chaining.

    Raises:
        ValueError: If *x* has fewer than 2 columns.
    """
    df = x if isinstance(x, pd.DataFrame) else pd.DataFrame(x)
    if df.shape[1] < 2:
        raise ValueError("X must have at least 2 columns.")

    self.feature_names_ = list(df.columns.astype(str))

    rho = self._estimate(df)

    _check_correlation_bounds(rho)
    self.correlation_matrix_ = pd.DataFrame(rho, index=self.feature_names_, columns=self.feature_names_)
    return self

Sparse Gaussian graphical model for mixed continuous/ordinal data.

Estimates a precision matrix and the associated :class:~humlpy.graphs.UGRAPH by:

Computing the sample correlation matrix via :class:SampleCorrelation.
Projecting to the positive-definite cone if necessary.
Computing the full graphical lasso regularisation path.
Selecting the model with the smallest eBIC via :func:omega_select.
Re-fitting graphical lasso at the selected alpha for a clean final solution.

The overall design mirrors sklearn.covariance.GraphicalLassoCV but replaces cross-validation with eBIC-based selection.

Parameters:

Name	Type	Description	Default
`n_lambdas`	`int`	Number of regularisation steps in the path. Defaults to 50.	`50`
`ebic_gamma`	`float`	eBIC hyper-parameter :math:`\gamma`. `0` gives standard BIC; larger values impose a stronger dimensionality penalty. Defaults to 0.1.	`0.1`
`n_levels_threshold`	`int`	Unique-value threshold for declaring a column ordinal. Defaults to 20.	`20`

Attributes set after :meth:fit: precision_matrix_ (pd.DataFrame): Partial correlation matrix derived from the estimated precision (diagonal = 1). correlation_matrix_ (pd.DataFrame): Sample correlation matrix used as input to graphical lasso. graph_ (UGRAPH): Undirected graph whose edges correspond to non-zero off-diagonal entries in precision_matrix_. alpha_ (float): Selected regularisation parameter. ebic_scores_ (np.ndarray): eBIC values for each point on the path. singular_ (bool): Whether positive-definite projection was needed. feature_names_ (list[str]): Variable names inferred from the input.

Example::

mgl = MixedGraphicalLasso(ebic_gamma=0.1)
mgl.fit(X)
print(mgl.graph_.edges)
print(mgl.precision_matrix_)

Source code in humlpy/estimation.py

class MixedGraphicalLasso:
    r"""Sparse Gaussian graphical model for mixed continuous/ordinal data.

    Estimates a precision matrix and the associated :class:`~humlpy.graphs.UGRAPH`
    by:

    1. Computing the sample correlation matrix via :class:`SampleCorrelation`.
    2. Projecting to the positive-definite cone if necessary.
    3. Computing the full graphical lasso regularisation path.
    4. Selecting the model with the smallest eBIC via :func:`omega_select`.
    5. Re-fitting graphical lasso at the selected alpha for a clean final solution.

    The overall design mirrors ``sklearn.covariance.GraphicalLassoCV`` but
    replaces cross-validation with eBIC-based selection.

    Args:
        n_lambdas: Number of regularisation steps in the path.  Defaults to 50.
        ebic_gamma: eBIC hyper-parameter :math:`\gamma`.  ``0`` gives standard
            BIC; larger values impose a stronger dimensionality penalty.
            Defaults to 0.1.
        n_levels_threshold: Unique-value threshold for declaring a column
            ordinal.  Defaults to 20.

    Attributes set after :meth:`fit`:
        precision_matrix_ (pd.DataFrame): Partial correlation matrix derived
            from the estimated precision (diagonal = 1).
        correlation_matrix_ (pd.DataFrame): Sample correlation matrix used as
            input to graphical lasso.
        graph_ (UGRAPH): Undirected graph whose edges correspond to non-zero
            off-diagonal entries in *precision_matrix_*.
        alpha_ (float): Selected regularisation parameter.
        ebic_scores_ (np.ndarray): eBIC values for each point on the path.
        singular_ (bool): Whether positive-definite projection was needed.
        feature_names_ (list[str]): Variable names inferred from the input.

    Example::

        mgl = MixedGraphicalLasso(ebic_gamma=0.1)
        mgl.fit(X)
        print(mgl.graph_.edges)
        print(mgl.precision_matrix_)
    """

    def __init__(
        self,
        n_lambdas: int = 50,
        ebic_gamma: float = 0.1,
        n_levels_threshold: int = 20,
    ) -> None:
        """Construct a MixedGraphicalLasso estimator.

        Args:
            n_lambdas: Length of the regularisation path.  Defaults to 50.
            ebic_gamma: eBIC dimensionality penalty weight in [0, 1].
                Defaults to 0.1.
            n_levels_threshold: Ordinal detection threshold.  Defaults to 20.

        Raises:
            ValueError: If *ebic_gamma* is not in [0, 1].
        """
        if not (0.0 <= ebic_gamma <= 1.0):
            raise ValueError(f"`ebic_gamma` must be in [0, 1], got {ebic_gamma}.")

        self.n_lambdas = n_lambdas
        self.ebic_gamma = ebic_gamma
        self.n_levels_threshold = n_levels_threshold

        # set after fit
        self.precision_matrix_: pd.DataFrame | None = None
        self.correlation_matrix_: pd.DataFrame | None = None
        self.graph_: UGRAPH | None = None
        self.alpha_: float | None = None
        self.ebic_scores_: NDArray[Any] | None = None
        self.singular_: bool | None = None
        self.feature_names_: list[str] | None = None

    # ------------------------------------------------------------------
    # Public API
    # ------------------------------------------------------------------

    def fit(self, x: pd.DataFrame | NDArray[Any]) -> MixedGraphicalLasso:
        """Fit the model to *x*.

        Args:
            x: Data matrix of shape ``(n_samples, n_features)``.
                A :class:`pandas.DataFrame` with named columns is preferred.

        Returns:
            *self* -- enables method chaining.
        """
        df = x if isinstance(x, pd.DataFrame) else pd.DataFrame(x)
        n, _ = df.shape

        # ----------------------------------------------------------------
        # Step 1 -- sample correlation
        # ----------------------------------------------------------------
        sc = SampleCorrelation(
            n_levels_threshold=self.n_levels_threshold,
        ).fit(df)

        self.feature_names_ = sc.feature_names_
        rho_raw: NDArray[Any] = sc.correlation_matrix_.to_numpy()  # type: ignore[union-attr]

        # ----------------------------------------------------------------
        # Step 2 -- ensure positive definiteness
        # ----------------------------------------------------------------
        rho, singular = _make_positive_definite(rho_raw, keep_diag=True)
        self.singular_ = singular
        self.correlation_matrix_ = pd.DataFrame(rho, index=self.feature_names_, columns=self.feature_names_)

        # ----------------------------------------------------------------
        # Step 3 -- regularisation path
        # ----------------------------------------------------------------
        precision_path, lambda_path = _glasso_path(rho, n_lambdas=self.n_lambdas)

        # ----------------------------------------------------------------
        # Step 4 -- eBIC model selection
        # ----------------------------------------------------------------
        _, selected_alpha, ebic_scores = omega_select(
            precision_path,
            lambda_path,
            n=n,
            s=rho,
            gamma=self.ebic_gamma,
        )
        self.alpha_ = selected_alpha
        self.ebic_scores_ = ebic_scores

        # ----------------------------------------------------------------
        # Step 5 -- final fit at selected alpha
        # ----------------------------------------------------------------
        try:
            _, omega_final = graphical_lasso(rho, alpha=selected_alpha, mode="cd", max_iter=200)
        except (FloatingPointError, ValueError):
            logger.warning(
                "Final graphical_lasso fit failed for alpha=%.4f; falling back to path estimate.",
                selected_alpha,
            )
            best_idx = int(np.argmin(ebic_scores))
            omega_final = precision_path[best_idx]

        omega_pcor = _precision_to_partial(omega_final)

        # ----------------------------------------------------------------
        # Store results
        # ----------------------------------------------------------------
        names = self.feature_names_
        assert names is not None  # guaranteed by SampleCorrelation.fit
        self.precision_matrix_ = pd.DataFrame(omega_pcor, index=names, columns=names)
        self.graph_ = self._build_ugraph(omega_pcor, names)

        return self

    @property
    def n_edges_(self) -> int:
        """Number of edges in the estimated graph.

        Raises:
            RuntimeError: If :meth:`fit` has not been called.
        """
        if self.graph_ is None:
            raise RuntimeError("Call .fit() before accessing n_edges_.")
        return self.graph_.num_edges

    # ------------------------------------------------------------------
    # Internals
    # ------------------------------------------------------------------

    @staticmethod
    def _build_ugraph(pcor: NDArray[Any], names: list[str]) -> UGRAPH:
        """Construct a UGRAPH from a partial correlation matrix.

        An edge ``(i, j)`` is included whenever ``pcor[i, j] != 0``.

        Args:
            pcor: Partial correlation matrix (diagonal = 1).
            names: Variable names corresponding to rows/columns.

        Returns:
            UGRAPH with one node per variable and one undirected edge per
            non-zero off-diagonal entry.
        """
        d = pcor.shape[0]
        edges: list[tuple[str, str]] = [
            (names[i], names[j]) for i in range(d - 1) for j in range(i + 1, d) if pcor[i, j] != 0.0
        ]
        return UGRAPH(nodes=names, edges=edges)

`n_edges_` `property`

Number of edges in the estimated graph.

Raises:

Type	Description
`RuntimeError`	If :meth:`fit` has not been called.

`init(n_lambdas=50, ebic_gamma=0.1, n_levels_threshold=20)`

Construct a MixedGraphicalLasso estimator.

Parameters:

Name	Type	Description	Default
`n_lambdas`	`int`	Length of the regularisation path. Defaults to 50.	`50`
`ebic_gamma`	`float`	eBIC dimensionality penalty weight in [0, 1]. Defaults to 0.1.	`0.1`
`n_levels_threshold`	`int`	Ordinal detection threshold. Defaults to 20.	`20`

Raises:

Type	Description
`ValueError`	If ebic_gamma is not in [0, 1].

Source code in humlpy/estimation.py

def __init__(
    self,
    n_lambdas: int = 50,
    ebic_gamma: float = 0.1,
    n_levels_threshold: int = 20,
) -> None:
    """Construct a MixedGraphicalLasso estimator.

    Args:
        n_lambdas: Length of the regularisation path.  Defaults to 50.
        ebic_gamma: eBIC dimensionality penalty weight in [0, 1].
            Defaults to 0.1.
        n_levels_threshold: Ordinal detection threshold.  Defaults to 20.

    Raises:
        ValueError: If *ebic_gamma* is not in [0, 1].
    """
    if not (0.0 <= ebic_gamma <= 1.0):
        raise ValueError(f"`ebic_gamma` must be in [0, 1], got {ebic_gamma}.")

    self.n_lambdas = n_lambdas
    self.ebic_gamma = ebic_gamma
    self.n_levels_threshold = n_levels_threshold

    # set after fit
    self.precision_matrix_: pd.DataFrame | None = None
    self.correlation_matrix_: pd.DataFrame | None = None
    self.graph_: UGRAPH | None = None
    self.alpha_: float | None = None
    self.ebic_scores_: NDArray[Any] | None = None
    self.singular_: bool | None = None
    self.feature_names_: list[str] | None = None

`fit(x)`

Fit the model to x.

Parameters:

Name	Type	Description	Default
`x`	`DataFrame \| NDArray[Any]`	Data matrix of shape `(n_samples, n_features)`. A :class:`pandas.DataFrame` with named columns is preferred.	required

Returns:

Type	Description
`MixedGraphicalLasso`	self -- enables method chaining.

Source code in humlpy/estimation.py

def fit(self, x: pd.DataFrame | NDArray[Any]) -> MixedGraphicalLasso:
    """Fit the model to *x*.

    Args:
        x: Data matrix of shape ``(n_samples, n_features)``.
            A :class:`pandas.DataFrame` with named columns is preferred.

    Returns:
        *self* -- enables method chaining.
    """
    df = x if isinstance(x, pd.DataFrame) else pd.DataFrame(x)
    n, _ = df.shape

    # ----------------------------------------------------------------
    # Step 1 -- sample correlation
    # ----------------------------------------------------------------
    sc = SampleCorrelation(
        n_levels_threshold=self.n_levels_threshold,
    ).fit(df)

    self.feature_names_ = sc.feature_names_
    rho_raw: NDArray[Any] = sc.correlation_matrix_.to_numpy()  # type: ignore[union-attr]

    # ----------------------------------------------------------------
    # Step 2 -- ensure positive definiteness
    # ----------------------------------------------------------------
    rho, singular = _make_positive_definite(rho_raw, keep_diag=True)
    self.singular_ = singular
    self.correlation_matrix_ = pd.DataFrame(rho, index=self.feature_names_, columns=self.feature_names_)

    # ----------------------------------------------------------------
    # Step 3 -- regularisation path
    # ----------------------------------------------------------------
    precision_path, lambda_path = _glasso_path(rho, n_lambdas=self.n_lambdas)

    # ----------------------------------------------------------------
    # Step 4 -- eBIC model selection
    # ----------------------------------------------------------------
    _, selected_alpha, ebic_scores = omega_select(
        precision_path,
        lambda_path,
        n=n,
        s=rho,
        gamma=self.ebic_gamma,
    )
    self.alpha_ = selected_alpha
    self.ebic_scores_ = ebic_scores

    # ----------------------------------------------------------------
    # Step 5 -- final fit at selected alpha
    # ----------------------------------------------------------------
    try:
        _, omega_final = graphical_lasso(rho, alpha=selected_alpha, mode="cd", max_iter=200)
    except (FloatingPointError, ValueError):
        logger.warning(
            "Final graphical_lasso fit failed for alpha=%.4f; falling back to path estimate.",
            selected_alpha,
        )
        best_idx = int(np.argmin(ebic_scores))
        omega_final = precision_path[best_idx]

    omega_pcor = _precision_to_partial(omega_final)

    # ----------------------------------------------------------------
    # Store results
    # ----------------------------------------------------------------
    names = self.feature_names_
    assert names is not None  # guaranteed by SampleCorrelation.fit
    self.precision_matrix_ = pd.DataFrame(omega_pcor, index=names, columns=names)
    self.graph_ = self._build_ugraph(omega_pcor, names)

    return self

Model Selection

Select the best precision matrix from a glasso path using eBIC.

The extended BIC criterion is:

.. math::

\mathrm{eBIC}(\Omega) =
    -2 \ell(\Omega \mid E)
    + |E| \log n
    + 4 \gamma |E| \log d

where :math:|E| is the number of edges, :math:n the sample size, :math:d the number of variables, and :math:\gamma \in [0,1] an additional penalty for high-dimensional settings.

Parameters:

Name	Type	Description	Default
`precision_path`	`list[NDArray[Any]]`	List of precision matrices along the regularisation path.	required
`lambda_path`	`NDArray[Any]`	Corresponding array of `alpha` values (same length).	required
`n`	`int`	Sample size.	required
`s`	`NDArray[Any]`	Sample correlation/covariance matrix used for estimation.	required
`gamma`	`float`	eBIC hyper-parameter. `0` recovers standard BIC. Defaults to 0.1.	`0.1`

Returns:

Type	Description
`NDArray[Any]`	Tuple `(selected_precision, selected_alpha, ebic_scores)` where
`float`	selected_precision is the raw precision matrix for the chosen model,
`NDArray[Any]`	selected_alpha is the corresponding regularisation value, and
`tuple[NDArray[Any], float, NDArray[Any]]`	ebic_scores is the full eBIC array along the path.

References

Foygel, Rina and Drton, Mathias. (2010). Extended Bayesian Information Criteria for Gaussian Graphical Models. Advances in Neural Information Processing Systems, Volume 23, pp. 604-612.

Source code in humlpy/estimation.py

def omega_select(
    precision_path: list[NDArray[Any]],
    lambda_path: NDArray[Any],
    n: int,
    s: NDArray[Any],
    *,
    gamma: float = 0.1,
) -> tuple[NDArray[Any], float, NDArray[Any]]:
    r"""Select the best precision matrix from a glasso path using eBIC.

    The extended BIC criterion is:

    .. math::

        \mathrm{eBIC}(\Omega) =
            -2 \ell(\Omega \mid E)
            + |E| \log n
            + 4 \gamma |E| \log d

    where :math:`|E|` is the number of edges, :math:`n` the sample size,
    :math:`d` the number of variables, and :math:`\gamma \in [0,1]` an
    additional penalty for high-dimensional settings.

    Args:
        precision_path: List of precision matrices along the regularisation path.
        lambda_path: Corresponding array of ``alpha`` values (same length).
        n: Sample size.
        s: Sample correlation/covariance matrix used for estimation.
        gamma: eBIC hyper-parameter.  ``0`` recovers standard BIC.
            Defaults to 0.1.

    Returns:
        Tuple ``(selected_precision, selected_alpha, ebic_scores)`` where
        *selected_precision* is the raw precision matrix for the chosen model,
        *selected_alpha* is the corresponding regularisation value, and
        *ebic_scores* is the full eBIC array along the path.

    References:
        Foygel, Rina and Drton, Mathias. (2010).
        Extended Bayesian Information Criteria for Gaussian Graphical Models.
        Advances in Neural Information Processing Systems, Volume 23, pp. 604-612.
    """
    d = s.shape[0]
    ebic = np.zeros(len(lambda_path))

    for k, omega in enumerate(precision_path):
        n_edges = _edgenumber(omega)
        sign, logdet = np.linalg.slogdet(omega)
        logdet = logdet if sign > 0 else -np.inf
        loglik = 0.5 * n * (logdet - np.trace(s @ omega))
        ebic[k] = -2 * loglik + n_edges * np.log(n) + 4 * gamma * n_edges * np.log(d)

    best = int(np.argmin(ebic))
    return precision_path[best].copy(), float(lambda_path[best]), ebic

Graph Classes

Bases: GRAPH

Class for dealing with undirected graph i.e. graphs that only contain undirected edges.

Source code in humlpy/graphs.py

class UGRAPH(GRAPH):
    """Class for dealing with undirected graph i.e. graphs that only contain undirected edges."""

    def __init__(
        self,
        nodes: list[str] | None = None,
        edges: list[tuple[str, str]] | None = None,
    ) -> None:
        """UGRAPH constructor.

        Args:
            nodes (list[str] | None, optional): Nodes. Defaults to None.
            edges (list[tuple[str,str]] | None, optional): Edges. Defaults to None.
        """
        if nodes is None:
            nodes = []
        if edges is None:
            edges = []

        self._nodes: set[str] = set(nodes)
        self._edges: set[tuple[str, str]] = set()
        self._neighbors: defaultdict[str, set[str]] = defaultdict(set)

        for edge in edges:
            self._add_edge(*edge)

    def _add_edge(self, i: str, j: str) -> None:
        self._nodes.add(i)
        self._nodes.add(j)
        self._edges.add((i, j))

        self._neighbors[i].add(j)
        self._neighbors[j].add(i)

    def neighbors(self, node: str) -> set[str]:
        """Gives all neighbors of node `node`.

        Args:
            node (str): node in current UGRAPH.

        Returns:
            set: set of neighbors.
        """
        if node in self._neighbors:
            return self._neighbors[node]
        else:
            return set()

    def is_adjacent(self, i: str, j: str) -> bool:
        """Return True if the graph contains an undirected edge between i and j.

        Args:
            i (str): node i.
            j (str): node j.

        Returns:
            bool: True if i - j
        """
        return (i, j) in self._edges or (j, i) in self._edges

    def is_clique(self, potential_clique: set[str]) -> bool:
        """Check every pair of nodes in potential_clique is adjacent."""
        return all(self.is_adjacent(i, j) for i, j in combinations(potential_clique, 2))

    @classmethod
    def from_pandas_adjacency(cls, pd_amat: pd.DataFrame) -> UGRAPH:
        """Build UGRAPH from a Pandas adjacency matrix.

        Args:
            pd_amat (pd.DataFrame): input adjacency matrix.

        Returns:
            UGRAPH
        """
        assert pd_amat.shape[0] == pd_amat.shape[1]
        nodes = list(pd_amat.columns)

        all_connections = []
        start, end = np.where(pd_amat != 0)
        for idx, _ in enumerate(start):
            all_connections.append((pd_amat.columns[start[idx]], pd_amat.columns[end[idx]]))

        edges = [tuple(item) for item in set(frozenset(item) for item in all_connections)]

        return UGRAPH(nodes=nodes, edges=edges)

    def remove_edge(self, i: str, j: str) -> None:
        """Removes edge in question.

        Args:
            i (str): first node
            j (str): second node

        Raises:
            AssertionError: if edge does not exist
        """
        if not self.is_adjacent(i, j):
            raise AssertionError("Edge does not exist in current UGRAPH")

        self._edges.discard((i, j))
        self._edges.discard((j, i))
        self._neighbors[i].discard(j)
        self._neighbors[j].discard(i)

    def remove_node(self, node: str) -> None:
        """Remove a node from the graph.

        Args:
            node (str): node to remove
        """
        self._nodes.remove(node)

        self._edges = {(i, j) for i, j in self._edges if node not in {i, j}}

        for nbr in self._neighbors[node]:
            self._neighbors[nbr].discard(node)

        self._neighbors.pop(node, "I was never here")

    @property
    def adjacency_matrix(self) -> pd.DataFrame:
        """Returns adjacency matrix.

        The i,jth entry being one indicates that there is an undirected edge
        between i and j. A zero indicates that there is no edge. The matrix
        is symmetric.

        Returns:
            pd.DataFrame: adjacency matrix
        """
        amat = pd.DataFrame(
            np.zeros([self.num_nodes, self.num_nodes]),
            index=self.nodes,
            columns=self.nodes,
        )
        for edge in self.edges:
            amat.loc[edge] = amat.loc[edge[::-1]] = 1
        return amat

    @property
    def causal_order(self) -> None:
        """Causal order is None.

        This is because undirected graphs do not imply a causal order.

        Returns:
            None: None
        """
        return None

    def copy(self) -> UGRAPH:
        """Return a copy of the graph."""
        return UGRAPH(nodes=list(self._nodes), edges=list(self._edges))

    def show(self) -> None:
        """Plot UGRAPH."""
        graph = self.to_networkx()
        pos = nx.circular_layout(graph)
        nx.draw(graph, pos=pos, with_labels=True)

    def to_networkx(self) -> nx.Graph:
        """Convert to networkx graph.

        Returns:
            nx.Graph: Undirected networkx graph.
        """
        nx_ugraph = nx.Graph()
        nx_ugraph.add_nodes_from(self.nodes)
        nx_ugraph.add_edges_from(self.edges)
        return nx_ugraph

    @property
    def nodes(self) -> list[str]:
        """Get all nodes in current UGRAPH.

        Returns:
            list: list of nodes.
        """
        return sorted(list(self._nodes))

    @property
    def num_nodes(self) -> int:
        """Number of nodes in current UGRAPH.

        Returns:
            int: Number of nodes
        """
        return len(self._nodes)

    @property
    def num_edges(self) -> int:
        """Number of edges in current UGRAPH.

        Returns:
            int: Number of edges
        """
        return len(self._edges)

    @property
    def edges(self) -> list[tuple[str, str]]:
        """Gives all edges in current UGRAPH.

        Returns:
            list[tuple[str,str]]: List of edges.
        """
        return list(self._edges)

`adjacency_matrix` `property`

Returns adjacency matrix.

The i,jth entry being one indicates that there is an undirected edge between i and j. A zero indicates that there is no edge. The matrix is symmetric.

Returns:

Type	Description
`DataFrame`	pd.DataFrame: adjacency matrix

`causal_order` `property`

Causal order is None.

This is because undirected graphs do not imply a causal order.

Returns:

Name	Type	Description
`None`	`None`	None

`edges` `property`

Gives all edges in current UGRAPH.

Returns:

Type	Description
`list[tuple[str, str]]`	list[tuple[str,str]]: List of edges.

`nodes` `property`

Get all nodes in current UGRAPH.

Returns:

Name	Type	Description
`list`	`list[str]`	list of nodes.

`num_edges` `property`

Number of edges in current UGRAPH.

Returns:

Name	Type	Description
`int`	`int`	Number of edges

`num_nodes` `property`

Number of nodes in current UGRAPH.

Returns:

Name	Type	Description
`int`	`int`	Number of nodes

`init(nodes=None, edges=None)`

UGRAPH constructor.

Parameters:

Name	Type	Description	Default
`nodes`	`list[str] \| None`	Nodes. Defaults to None.	`None`
`edges`	`list[tuple[str, str]] \| None`	Edges. Defaults to None.	`None`

Source code in humlpy/graphs.py

def __init__(
    self,
    nodes: list[str] | None = None,
    edges: list[tuple[str, str]] | None = None,
) -> None:
    """UGRAPH constructor.

    Args:
        nodes (list[str] | None, optional): Nodes. Defaults to None.
        edges (list[tuple[str,str]] | None, optional): Edges. Defaults to None.
    """
    if nodes is None:
        nodes = []
    if edges is None:
        edges = []

    self._nodes: set[str] = set(nodes)
    self._edges: set[tuple[str, str]] = set()
    self._neighbors: defaultdict[str, set[str]] = defaultdict(set)

    for edge in edges:
        self._add_edge(*edge)

`copy()`

Return a copy of the graph.

Source code in humlpy/graphs.py

def copy(self) -> UGRAPH:
    """Return a copy of the graph."""
    return UGRAPH(nodes=list(self._nodes), edges=list(self._edges))

`from_pandas_adjacency(pd_amat)` `classmethod`

Build UGRAPH from a Pandas adjacency matrix.

Parameters:

Name	Type	Description	Default
`pd_amat`	`DataFrame`	input adjacency matrix.	required

Returns:

Type	Description
`UGRAPH`	UGRAPH

Source code in humlpy/graphs.py

@classmethod
def from_pandas_adjacency(cls, pd_amat: pd.DataFrame) -> UGRAPH:
    """Build UGRAPH from a Pandas adjacency matrix.

    Args:
        pd_amat (pd.DataFrame): input adjacency matrix.

    Returns:
        UGRAPH
    """
    assert pd_amat.shape[0] == pd_amat.shape[1]
    nodes = list(pd_amat.columns)

    all_connections = []
    start, end = np.where(pd_amat != 0)
    for idx, _ in enumerate(start):
        all_connections.append((pd_amat.columns[start[idx]], pd_amat.columns[end[idx]]))

    edges = [tuple(item) for item in set(frozenset(item) for item in all_connections)]

    return UGRAPH(nodes=nodes, edges=edges)

`is_adjacent(i, j)`

Return True if the graph contains an undirected edge between i and j.

Parameters:

Name	Type	Description	Default
`i`	`str`	node i.	required
`j`	`str`	node j.	required

Returns:

Name	Type	Description
`bool`	`bool`	True if i - j

Source code in humlpy/graphs.py

def is_adjacent(self, i: str, j: str) -> bool:
    """Return True if the graph contains an undirected edge between i and j.

    Args:
        i (str): node i.
        j (str): node j.

    Returns:
        bool: True if i - j
    """
    return (i, j) in self._edges or (j, i) in self._edges

`is_clique(potential_clique)`

Check every pair of nodes in potential_clique is adjacent.

Source code in humlpy/graphs.py

def is_clique(self, potential_clique: set[str]) -> bool:
    """Check every pair of nodes in potential_clique is adjacent."""
    return all(self.is_adjacent(i, j) for i, j in combinations(potential_clique, 2))

`neighbors(node)`

Gives all neighbors of node node.

Parameters:

Name	Type	Description	Default
`node`	`str`	node in current UGRAPH.	required

Returns:

Name	Type	Description
`set`	`set[str]`	set of neighbors.

Source code in humlpy/graphs.py

def neighbors(self, node: str) -> set[str]:
    """Gives all neighbors of node `node`.

    Args:
        node (str): node in current UGRAPH.

    Returns:
        set: set of neighbors.
    """
    if node in self._neighbors:
        return self._neighbors[node]
    else:
        return set()

`remove_edge(i, j)`

Removes edge in question.

Parameters:

Name	Type	Description	Default
`i`	`str`	first node	required
`j`	`str`	second node	required

Raises:

Type	Description
`AssertionError`	if edge does not exist

Source code in humlpy/graphs.py

def remove_edge(self, i: str, j: str) -> None:
    """Removes edge in question.

    Args:
        i (str): first node
        j (str): second node

    Raises:
        AssertionError: if edge does not exist
    """
    if not self.is_adjacent(i, j):
        raise AssertionError("Edge does not exist in current UGRAPH")

    self._edges.discard((i, j))
    self._edges.discard((j, i))
    self._neighbors[i].discard(j)
    self._neighbors[j].discard(i)

`remove_node(node)`

Remove a node from the graph.

Parameters:

Name	Type	Description	Default
`node`	`str`	node to remove	required

Source code in humlpy/graphs.py

def remove_node(self, node: str) -> None:
    """Remove a node from the graph.

    Args:
        node (str): node to remove
    """
    self._nodes.remove(node)

    self._edges = {(i, j) for i, j in self._edges if node not in {i, j}}

    for nbr in self._neighbors[node]:
        self._neighbors[nbr].discard(node)

    self._neighbors.pop(node, "I was never here")

`show()`

Plot UGRAPH.

Source code in humlpy/graphs.py

def show(self) -> None:
    """Plot UGRAPH."""
    graph = self.to_networkx()
    pos = nx.circular_layout(graph)
    nx.draw(graph, pos=pos, with_labels=True)

`to_networkx()`

Convert to networkx graph.

Returns:

Type	Description
`Graph`	nx.Graph: Undirected networkx graph.

Source code in humlpy/graphs.py

def to_networkx(self) -> nx.Graph:
    """Convert to networkx graph.

    Returns:
        nx.Graph: Undirected networkx graph.
    """
    nx_ugraph = nx.Graph()
    nx_ugraph.add_nodes_from(self.nodes)
    nx_ugraph.add_edges_from(self.edges)
    return nx_ugraph

Correlation Functions

Bases: CorrelationMeasure

Maximum-likelihood polychoric correlation between two ordinal variables.

Both variables are treated as discretised versions of latent normal variates. The MLE of their latent correlation is found by one of two solvers:

"brent" (default) — Direct maximisation of $\ell(\rho) = \sum_{ij} n_{ij}\log\pi_{ij}$ via :func:scipy.optimize.minimize_scalar with Brent's bounded method. Each function evaluation requires only :func:_pi_rs (no derivatives). Converges super-linearly and is consistently faster once the number of ordinal categories exceeds ~3, because the per-evaluation cost grows as $O(k_x \times k_y)$ CDF calls rather than CDF + PDF calls.
"newton" — Fisher scoring. Uses the score $S(\rho) = \sum_{ij} n_{ij}\,(\pi_{ij}'/\pi_{ij})$ and the empirical Fisher information $\hat{I}(\rho) = \sum_{ij} n_{ij}\,(\pi_{ij}'/\pi_{ij})^2$ to iterate $\rho \leftarrow \rho + S/\hat{I}$ until convergence. Converges quadratically near the MLE and can be ~2× faster than "brent" for binary variables ($k_x = k_y = 2$), but becomes progressively slower as the category count grows because each iteration requires both :func:_pi_rs and :func:_pi_rs_derivative per cell. Consider this solver when both variables are binary or ternary.

Solver selection guide (rule of thumb from benchmarks):

Parameters:

Name	Type	Description	Default
`max_cor`	`float`	Absolute upper bound for the estimated correlation. Defaults to 0.9999.	`0.9999`
`solver`	`Literal['newton', 'brent']`	Optimisation strategy, `"brent"` or `"newton"`. Defaults to `"brent"`.	`'brent'`
`max_iter`	`int`	Maximum Fisher-scoring iterations (ignored for `"brent"`). Defaults to 100.	`100`
`tol`	`float`	Convergence tolerance on the absolute score value and step size (ignored for `"brent"`). Defaults to 1e-10.	`1e-10`

Example::

rho = PolychoricCorrelation().fit(x_ord, y_ord).correlation
rho = PolychoricCorrelation(solver="newton").fit(x_ord, y_ord).correlation

Source code in humlpy/correlation.py

class PolychoricCorrelation(CorrelationMeasure):
    r"""Maximum-likelihood polychoric correlation between two ordinal variables.

    Both variables are treated as discretised versions of latent normal variates.
    The MLE of their latent correlation is found by one of two solvers:

    * ``"brent"`` (**default**) — Direct maximisation of
      $\\ell(\\rho) = \\sum_{ij} n_{ij}\\log\\pi_{ij}$ via
      :func:`scipy.optimize.minimize_scalar` with Brent's bounded method.
      Each function evaluation requires only :func:`_pi_rs` (no derivatives).
      Converges super-linearly and is consistently faster once the number of
      ordinal categories exceeds ~3, because the per-evaluation cost grows as
      $O(k_x \\times k_y)$ CDF calls rather than CDF + PDF calls.

    * ``"newton"`` — Fisher scoring.  Uses the score
      $S(\\rho) = \\sum_{ij} n_{ij}\\,(\\pi_{ij}'/\\pi_{ij})$ and the empirical
      Fisher information $\\hat{I}(\\rho) = \\sum_{ij} n_{ij}\\,(\\pi_{ij}'/\\pi_{ij})^2$
      to iterate $\\rho \\leftarrow \\rho + S/\\hat{I}$ until convergence.
      Converges quadratically near the MLE and can be ~2× faster than
      ``"brent"`` for binary variables ($k_x = k_y = 2$), but becomes
      progressively slower as the category count grows because each iteration
      requires both :func:`_pi_rs` and :func:`_pi_rs_derivative` per cell.
      Consider this solver when both variables are binary or ternary.

    **Solver selection guide** (rule of thumb from benchmarks):

    +-------------------+------------------+
    | Category count    | Recommended      |
    +===================+==================+
    | binary (2 × 2)    | ``"newton"``     |
    +-------------------+------------------+
    | ternary (3 × 3)   | either           |
    +-------------------+------------------+
    | 4+ categories     | ``"brent"``      |
    +-------------------+------------------+

    Args:
        max_cor: Absolute upper bound for the estimated correlation.
            Defaults to 0.9999.
        solver: Optimisation strategy, ``"brent"`` or ``"newton"``.
            Defaults to ``"brent"``.
        max_iter: Maximum Fisher-scoring iterations (ignored for ``"brent"``).
            Defaults to 100.
        tol: Convergence tolerance on the absolute score value and step size
            (ignored for ``"brent"``).
            Defaults to 1e-10.

    Example::

        rho = PolychoricCorrelation().fit(x_ord, y_ord).correlation
        rho = PolychoricCorrelation(solver="newton").fit(x_ord, y_ord).correlation
    """

    def __init__(
        self,
        max_cor: float = 0.9999,
        solver: Literal["newton", "brent"] = "brent",
        max_iter: int = 100,
        tol: float = 1e-10,
    ) -> None:
        """Initialise the polychoric correlation estimator.

        Args:
            max_cor: Absolute upper bound for the estimated correlation.
                Defaults to 0.9999.
            solver: Optimisation strategy.  ``"brent"`` (default) is faster
                for variables with 4 or more ordinal categories.
                ``"newton"`` (Fisher scoring) can be ~2× faster for binary
                variables but slows down as category count grows.
            max_iter: Maximum Fisher-scoring iterations; ignored when
                ``solver="brent"``.  Defaults to 100.
            tol: Convergence tolerance on the absolute score value and step
                size; ignored when ``solver="brent"``.  Defaults to 1e-10.

        Raises:
            ValueError: If *solver* is not ``"newton"`` or ``"brent"``.
        """
        super().__init__(max_cor=max_cor)
        if solver not in {"newton", "brent"}:
            raise ValueError(f"`solver` must be 'newton' or 'brent', got '{solver}'.")
        self._solver = solver
        self._max_iter = max_iter
        self._tol = tol

    def fit(self, x: np.ndarray, y: np.ndarray) -> PolychoricCorrelation:
        """Fit the polychoric correlation model.

        Args:
            x: First ordinal variable (integer-valued or small number of
                distinct float levels).
            y: Second ordinal variable.

        Returns:
            *self*

        Raises:
            ValueError: On invalid or incompatible inputs.
            RuntimeError: If the solver cannot converge.
        """
        x_arr, y_arr = self._prepare(x, y)
        _validate_ordinal(x_arr, "x")
        _validate_ordinal(y_arr, "y")

        self._correlation = self._clip(self._polychoric(x_arr, y_arr))
        return self

    # ------------------------------------------------------------------
    # Dispatch
    # ------------------------------------------------------------------

    def _polychoric(self, x: np.ndarray, y: np.ndarray) -> float:
        """Build shared data structures, then dispatch to the chosen solver."""
        n = x.size
        ux = np.unique(x)
        uy = np.unique(y)

        n_rs = np.array(
            [[np.sum((x == xi) & (y == yj)) for yj in uy] for xi in ux],
            dtype=float,
        )
        assert n_rs.sum() == n, "Joint frequency table does not sum to n."

        tx = _thresholds(x)
        ty = _thresholds(y)

        if self._solver == "newton":
            return self._polychoric_newton(n_rs, tx, ty, ux, uy)
        return self._polychoric_brent(n_rs, tx, ty, ux, uy)

    # ------------------------------------------------------------------
    # Solver: Fisher scoring (Newton-type)
    # ------------------------------------------------------------------

    def _polychoric_newton(
        self,
        n_rs: np.ndarray,
        tx: np.ndarray,
        ty: np.ndarray,
        ux: np.ndarray,
        uy: np.ndarray,
    ) -> float:
        r"""Fisher scoring iteration on the score equation.

        Update rule: $\\rho \\leftarrow \\rho + S(\\rho)/\\hat{I}(\\rho)$

        where $S = \\sum n_{ij}\\,(dp/p)$ is the score and
        $\\hat{I} = \\sum n_{ij}\\,(dp/p)^2$ is the empirical Fisher information.
        Both $\\pi_{ij}$ and $\\pi_{ij}'$ are needed per cell per iteration;
        the step converges quadratically in the neighbourhood of the MLE.
        """
        rho = 0.0
        bound = self._max_cor
        score_val = 0.0

        for iteration in range(self._max_iter):
            score_val = 0.0
            info_val = 0.0

            for i in range(len(ux)):
                for j in range(len(uy)):
                    lower = (float(tx[i]), float(ty[j]))
                    upper = (float(tx[i + 1]), float(ty[j + 1]))
                    p = max(_pi_rs(lower=lower, upper=upper, corr=rho), _CELL_FLOOR)
                    dp = _pi_rs_derivative(lower=np.array(lower), upper=np.array(upper), corr=rho)
                    ratio = dp / p
                    score_val += n_rs[i, j] * ratio
                    info_val += n_rs[i, j] * ratio**2

            if abs(score_val) < self._tol:
                break
            if info_val < 1e-14:  # nearly flat — cannot determine direction
                logger.warning(
                    "Fisher information ≈ 0 at ρ=%.4f after %d iteration(s); returning current estimate.",
                    rho,
                    iteration,
                )
                break

            step = score_val / info_val
            rho_new = float(np.clip(rho + step, -bound, bound))

            if abs(rho_new - rho) < self._tol:
                rho = rho_new
                break
            rho = rho_new
        else:
            logger.warning(
                "Fisher scoring did not converge in %d iterations (|score|=%.2e). "
                "Consider increasing max_iter or switching to solver='brent'.",
                self._max_iter,
                abs(score_val),
            )

        return rho

    # ------------------------------------------------------------------
    # Solver: direct log-likelihood via Brent's bounded method
    # ------------------------------------------------------------------

    def _polychoric_brent(
        self,
        n_rs: np.ndarray,
        tx: np.ndarray,
        ty: np.ndarray,
        ux: np.ndarray,
        uy: np.ndarray,
    ) -> float:
        r"""Maximise $\\ell(\\rho)$ directly via ``minimize_scalar`` (Brent bounded).

        Each function evaluation costs one :func:`_pi_rs` call per cell.
        No derivative information is used; convergence is super-linear.
        Near-zero cell probabilities are floored at *_CELL_FLOOR* so that
        empty cells never produce ``-inf`` contributions.
        """
        bound = self._max_cor

        def neg_log_likelihood(rho: float) -> float:
            total = 0.0
            for i in range(len(ux)):
                for j in range(len(uy)):
                    if n_rs[i, j] == 0:
                        continue
                    lower = (float(tx[i]), float(ty[j]))
                    upper = (float(tx[i + 1]), float(ty[j + 1]))
                    p = max(_pi_rs(lower=lower, upper=upper, corr=rho), _CELL_FLOOR)
                    total += n_rs[i, j] * np.log(p)
            return -total

        result = minimize_scalar(neg_log_likelihood, bounds=(-bound, bound), method="bounded")
        return float(result.x)

`init(max_cor=0.9999, solver='brent', max_iter=100, tol=1e-10)`

Initialise the polychoric correlation estimator.

Parameters:

Name	Type	Description	Default
`max_cor`	`float`	Absolute upper bound for the estimated correlation. Defaults to 0.9999.	`0.9999`
`solver`	`Literal['newton', 'brent']`	Optimisation strategy. `"brent"` (default) is faster for variables with 4 or more ordinal categories. `"newton"` (Fisher scoring) can be ~2× faster for binary variables but slows down as category count grows.	`'brent'`
`max_iter`	`int`	Maximum Fisher-scoring iterations; ignored when `solver="brent"`. Defaults to 100.	`100`
`tol`	`float`	Convergence tolerance on the absolute score value and step size; ignored when `solver="brent"`. Defaults to 1e-10.	`1e-10`

Raises:

Type	Description
`ValueError`	If solver is not `"newton"` or `"brent"`.

Source code in humlpy/correlation.py

def __init__(
    self,
    max_cor: float = 0.9999,
    solver: Literal["newton", "brent"] = "brent",
    max_iter: int = 100,
    tol: float = 1e-10,
) -> None:
    """Initialise the polychoric correlation estimator.

    Args:
        max_cor: Absolute upper bound for the estimated correlation.
            Defaults to 0.9999.
        solver: Optimisation strategy.  ``"brent"`` (default) is faster
            for variables with 4 or more ordinal categories.
            ``"newton"`` (Fisher scoring) can be ~2× faster for binary
            variables but slows down as category count grows.
        max_iter: Maximum Fisher-scoring iterations; ignored when
            ``solver="brent"``.  Defaults to 100.
        tol: Convergence tolerance on the absolute score value and step
            size; ignored when ``solver="brent"``.  Defaults to 1e-10.

    Raises:
        ValueError: If *solver* is not ``"newton"`` or ``"brent"``.
    """
    super().__init__(max_cor=max_cor)
    if solver not in {"newton", "brent"}:
        raise ValueError(f"`solver` must be 'newton' or 'brent', got '{solver}'.")
    self._solver = solver
    self._max_iter = max_iter
    self._tol = tol

`fit(x, y)`

Fit the polychoric correlation model.

Parameters:

Name	Type	Description	Default
`x`	`ndarray`	First ordinal variable (integer-valued or small number of distinct float levels).	required
`y`	`ndarray`	Second ordinal variable.	required

Returns:

Type	Description
`PolychoricCorrelation`	self

Raises:

Type	Description
`ValueError`	On invalid or incompatible inputs.
`RuntimeError`	If the solver cannot converge.

Source code in humlpy/correlation.py

def fit(self, x: np.ndarray, y: np.ndarray) -> PolychoricCorrelation:
    """Fit the polychoric correlation model.

    Args:
        x: First ordinal variable (integer-valued or small number of
            distinct float levels).
        y: Second ordinal variable.

    Returns:
        *self*

    Raises:
        ValueError: On invalid or incompatible inputs.
        RuntimeError: If the solver cannot converge.
    """
    x_arr, y_arr = self._prepare(x, y)
    _validate_ordinal(x_arr, "x")
    _validate_ordinal(y_arr, "y")

    self._correlation = self._clip(self._polychoric(x_arr, y_arr))
    return self

Bases: CorrelationMeasure

Ad-hoc polyserial correlation between one continuous and one ordinal variable.

The estimator applies the nonparanormal (rank-based) transformation to the continuous variable and uses a closed-form correction factor derived from the ordinal thresholds to approximate the underlying latent correlation.

Parameters:

Name	Type	Description	Default
`max_cor`	`float`	Absolute upper bound for the estimated correlation. Defaults to 0.9999.	`0.9999`
`n_levels_threshold`	`int`	Variables with fewer unique values than this threshold are treated as ordinal. Defaults to 20.	`20`

Example::

rho = PolyserialCorrelation().fit(x_cont, y_ord).correlation

Source code in humlpy/correlation.py

class PolyserialCorrelation(CorrelationMeasure):
    """Ad-hoc polyserial correlation between one continuous and one ordinal variable.

    The estimator applies the nonparanormal (rank-based) transformation to the
    continuous variable and uses a closed-form correction factor derived from
    the ordinal thresholds to approximate the underlying latent correlation.

    Args:
        max_cor: Absolute upper bound for the estimated correlation.
            Defaults to 0.9999.
        n_levels_threshold: Variables with fewer unique values than this
            threshold are treated as ordinal.  Defaults to 20.

    Example::

        rho = PolyserialCorrelation().fit(x_cont, y_ord).correlation
    """

    def __init__(
        self,
        max_cor: float = 0.9999,
        n_levels_threshold: int = 20,
    ) -> None:
        """Inits polyserial correlation class.

        Args:
            max_cor (float, optional): Defaults to 0.9999.
            n_levels_threshold (int, optional): Until what level set to treat a variable as ordinal. Defaults to 20.

        Raises:
            ValueError: _description_
        """
        super().__init__(max_cor=max_cor)
        if n_levels_threshold < 2:
            raise ValueError("`n_levels_threshold` must be ≥ 2.")
        self._n_levels_threshold = n_levels_threshold

    def fit(self, x: np.ndarray, y: np.ndarray) -> PolyserialCorrelation:
        """Fit the polyserial correlation model.

        The method automatically identifies which argument is the continuous
        variable and which is the ordinal variable based on the number of
        unique values.

        Args:
            x: One variable (continuous or ordinal).
            y: The other variable (ordinal or continuous).

        Returns:
            *self*

        Raises:
            ValueError: On invalid inputs or if both variables appear to be
                continuous (use a Spearman/NPP estimator instead) or if both
                appear to be ordinal (use :class:`PolychoricCorrelation`
                instead).
        """
        x_arr, y_arr = self._prepare(x, y)

        x_is_ord = len(np.unique(x_arr)) < self._n_levels_threshold
        y_is_ord = len(np.unique(y_arr)) < self._n_levels_threshold

        if not x_is_ord and not y_is_ord:
            raise ValueError(
                "Both variables appear continuous (≥ n_levels_threshold unique values). "
                "Use a Spearman or nonparanormal estimator instead."
            )
        if x_is_ord and y_is_ord:
            raise ValueError(
                "Both variables appear ordinal (< n_levels_threshold unique values). Use PolychoricCorrelation instead."
            )

        if x_is_ord:
            cont_arr, ord_arr = y_arr, x_arr
            cont_name, ord_name = "y", "x"
        else:
            cont_arr, ord_arr = x_arr, y_arr
            cont_name, ord_name = "x", "y"

        _validate_continuous(cont_arr, cont_name)
        _validate_ordinal(ord_arr, ord_name)

        self._correlation = self._clip(self._polyserial(cont_arr, ord_arr))
        return self

    # ------------------------------------------------------------------
    # Implementation
    # ------------------------------------------------------------------

    def _polyserial(self, cont: np.ndarray, disc: np.ndarray) -> float:
        unique_vals, _ = np.unique(disc, return_counts=True)
        threshold_estimate = _thresholds(disc)

        values = np.sort(unique_vals.astype(float))
        interior_thresholds = threshold_estimate[1:-1]
        value_diff = np.diff(values)
        lambda_val: float = float(np.sum(stats.norm.pdf(interior_thresholds) * value_diff))

        if lambda_val == 0:
            raise ValueError(
                "Denominator λ in ad-hoc polyserial estimator is zero. "
                "This can happen when ordinal categories are far apart in "
                "probability space.  Consider recoding the ordinal variable."
            )

        s_disc = float(np.std(disc.astype(float), ddof=1))
        r = _npn_pearson(cont, disc)
        return r * s_disc / lambda_val

`init(max_cor=0.9999, n_levels_threshold=20)`

Inits polyserial correlation class.

Parameters:

Name	Type	Description	Default
`max_cor`	`float`	Defaults to 0.9999.	`0.9999`
`n_levels_threshold`	`int`	Until what level set to treat a variable as ordinal. Defaults to 20.	`20`

Raises:

Type	Description
`ValueError`	description

Source code in humlpy/correlation.py

def __init__(
    self,
    max_cor: float = 0.9999,
    n_levels_threshold: int = 20,
) -> None:
    """Inits polyserial correlation class.

    Args:
        max_cor (float, optional): Defaults to 0.9999.
        n_levels_threshold (int, optional): Until what level set to treat a variable as ordinal. Defaults to 20.

    Raises:
        ValueError: _description_
    """
    super().__init__(max_cor=max_cor)
    if n_levels_threshold < 2:
        raise ValueError("`n_levels_threshold` must be ≥ 2.")
    self._n_levels_threshold = n_levels_threshold

`fit(x, y)`

Fit the polyserial correlation model.

The method automatically identifies which argument is the continuous variable and which is the ordinal variable based on the number of unique values.

Parameters:

Name	Type	Description	Default
`x`	`ndarray`	One variable (continuous or ordinal).	required
`y`	`ndarray`	The other variable (ordinal or continuous).	required

Returns:

Type	Description
`PolyserialCorrelation`	self

Raises:

Type	Description
`ValueError`	On invalid inputs or if both variables appear to be continuous (use a Spearman/NPP estimator instead) or if both appear to be ordinal (use :class:`PolychoricCorrelation` instead).

Source code in humlpy/correlation.py

def fit(self, x: np.ndarray, y: np.ndarray) -> PolyserialCorrelation:
    """Fit the polyserial correlation model.

    The method automatically identifies which argument is the continuous
    variable and which is the ordinal variable based on the number of
    unique values.

    Args:
        x: One variable (continuous or ordinal).
        y: The other variable (ordinal or continuous).

    Returns:
        *self*

    Raises:
        ValueError: On invalid inputs or if both variables appear to be
            continuous (use a Spearman/NPP estimator instead) or if both
            appear to be ordinal (use :class:`PolychoricCorrelation`
            instead).
    """
    x_arr, y_arr = self._prepare(x, y)

    x_is_ord = len(np.unique(x_arr)) < self._n_levels_threshold
    y_is_ord = len(np.unique(y_arr)) < self._n_levels_threshold

    if not x_is_ord and not y_is_ord:
        raise ValueError(
            "Both variables appear continuous (≥ n_levels_threshold unique values). "
            "Use a Spearman or nonparanormal estimator instead."
        )
    if x_is_ord and y_is_ord:
        raise ValueError(
            "Both variables appear ordinal (< n_levels_threshold unique values). Use PolychoricCorrelation instead."
        )

    if x_is_ord:
        cont_arr, ord_arr = y_arr, x_arr
        cont_name, ord_name = "y", "x"
    else:
        cont_arr, ord_arr = x_arr, y_arr
        cont_name, ord_name = "x", "y"

    _validate_continuous(cont_arr, cont_name)
    _validate_ordinal(ord_arr, ord_name)

    self._correlation = self._clip(self._polyserial(cont_arr, ord_arr))
    return self

Spearman's rank correlation coefficient.

Parameters:

Name	Type	Description	Default
`x`	`ndarray`	First numeric array.	required
`y`	`ndarray`	Second numeric array (same length as x).	required

Returns:

Type	Description
`float`	Spearman's ρ in [−1, 1].

Source code in humlpy/correlation.py

def spearman(x: np.ndarray, y: np.ndarray) -> float:
    """Spearman's rank correlation coefficient.

    Args:
        x: First numeric array.
        y: Second numeric array (same length as *x*).

    Returns:
        Spearman's ρ in [−1, 1].
    """
    x_arr, y_arr = _to_array(x), _to_array(y)
    _validate_pair(x_arr, y_arr)
    rho, _ = stats.spearmanr(x_arr, y_arr)  # pyright: ignore[reportArgumentType]
    return float(rho)

Winsorized nonparanormal transformation.

Implements the estimator from Liu, Lafferty & Wasserman (2009): each observation is mapped to its normal score and the result is scaled to unit variance. Extreme ranks are Winsorized to avoid ±∞.

Parameters:

Name	Type	Description	Default
`x`	`ndarray`	A 1-D numeric array.	required

Returns:

Type	Description
`ndarray`	Transformed array of the same length, scaled to unit variance.

Source code in humlpy/correlation.py

def f_hat(x: np.ndarray) -> np.ndarray:
    """Winsorized nonparanormal transformation.

    Implements the estimator from Liu, Lafferty & Wasserman (2009):  each
    observation is mapped to its normal score and the result is scaled to
    unit variance.  Extreme ranks are Winsorized to avoid ±∞.

    Args:
        x: A 1-D numeric array.

    Returns:
        Transformed array of the same length, scaled to unit variance.
    """
    x_arr = _to_array(x)
    if x_arr.size < 2:
        raise ValueError("f_hat requires at least 2 observations.")
    return _f_hat(x_arr)

Nonparanormal product-moment correlation.

Parameters:

Name	Type	Description	Default
`cont`	`ndarray`	Continuous numeric array.	required
`disc`	`ndarray`	Discrete/ordinal array.	required

Returns:

Type	Description
`float`	Correlation coefficient in [−1, 1].

Source code in humlpy/correlation.py

def npn_pearson(cont: np.ndarray, disc: np.ndarray) -> float:
    """Nonparanormal product-moment correlation.

    Args:
        cont: Continuous numeric array.
        disc: Discrete/ordinal array.

    Returns:
        Correlation coefficient in [−1, 1].
    """
    cont_arr, disc_arr = _to_array(cont), _to_array(disc)
    _validate_pair(cont_arr, disc_arr)
    return _npn_pearson(cont_arr, disc_arr)

Convenience dispatcher for mixed-type correlation estimation.

Selects the appropriate estimator based on the number of unique values:

Both continuous → nonparanormal Spearman sin-transformation.
Both ordinal → :class:PolychoricCorrelation.
Mixed → :class:PolyserialCorrelation.

Parameters:

Name	Type	Description	Default
`x`	`ndarray`	First array (discrete or continuous).	required
`y`	`ndarray`	Second array (discrete or continuous).	required
`max_cor`	`float`	Maximum allowed absolute correlation; clips to ±max_cor.	`0.9999`
`n_levels_threshold`	`int`	Unique-value count below which a variable is treated as ordinal.	`20`
`verbose`	`bool`	Emit an :mod:`logging` info message describing the chosen estimator.	`False`

Returns:

Type	Description
`float`	Correlation estimate in [−1, 1].

Source code in humlpy/correlation.py

def adhoc_polyserial(
    x: np.ndarray,
    y: np.ndarray,
    *,
    max_cor: float = 0.9999,
    n_levels_threshold: int = 20,
    verbose: bool = False,
) -> float:
    """Convenience dispatcher for mixed-type correlation estimation.

    Selects the appropriate estimator based on the number of unique values:

    * Both continuous → nonparanormal Spearman sin-transformation.
    * Both ordinal → :class:`PolychoricCorrelation`.
    * Mixed → :class:`PolyserialCorrelation`.

    Args:
        x: First array (discrete or continuous).
        y: Second array (discrete or continuous).
        max_cor: Maximum allowed absolute correlation; clips to ±*max_cor*.
        n_levels_threshold: Unique-value count below which a variable is
            treated as ordinal.
        verbose: Emit an :mod:`logging` info message describing the chosen
            estimator.

    Returns:
        Correlation estimate in [−1, 1].
    """
    x_arr, y_arr = _to_array(x), _to_array(y)
    _validate_pair(x_arr, y_arr)

    n_unique_x = len(np.unique(x_arr))
    n_unique_y = len(np.unique(y_arr))
    x_is_discrete = n_unique_x < n_levels_threshold
    y_is_discrete = n_unique_y < n_levels_threshold

    if not x_is_discrete and not y_is_discrete:
        if verbose:
            logger.info("Both variables continuous — using nonparanormal Spearman.")
        rho = spearman(x_arr, y_arr)
        return float(np.clip(2 * np.sin(np.pi / 6 * rho), -max_cor, max_cor))

    if x_is_discrete and y_is_discrete:
        if verbose:
            logger.info("Both variables ordinal — using polychoric correlation.")
        return PolychoricCorrelation(max_cor=max_cor).fit(x_arr, y_arr).correlation

    if verbose:
        logger.info("Mixed pair — using ad-hoc polyserial correlation.")
    return PolyserialCorrelation(max_cor=max_cor, n_levels_threshold=n_levels_threshold).fit(x_arr, y_arr).correlation

Legacy Functions

Estimate a mixed graphical model (backward-compatible wrapper).

.. deprecated:: Prefer :class:MixedGraphicalLasso directly.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame \| NDArray[Any]`	Data matrix `(n_samples, n_features)`.	required
`verbose`	`bool`	Emit logging warnings.	`True`
`n_lambdas`	`int`	Length of the regularisation path.	`50`
`param`	`float`	eBIC dimensionality penalty weight (gamma).	`0.1`
`feature_names`	`list[str] \| None`	Override column names.	`None`
`n_levels_threshold`	`int`	Ordinal detection threshold.	`20`

Returns:

Type	Description
`MixedGraphResult`	class:`MixedGraphResult`

Source code in humlpy/estimation.py

def mixed_graph_nonpara(
    data: pd.DataFrame | NDArray[Any],
    *,
    verbose: bool = True,
    n_lambdas: int = 50,
    param: float = 0.1,
    feature_names: list[str] | None = None,
    n_levels_threshold: int = 20,
) -> MixedGraphResult:
    """Estimate a mixed graphical model (backward-compatible wrapper).

    .. deprecated::
        Prefer :class:`MixedGraphicalLasso` directly.

    Args:
        data: Data matrix ``(n_samples, n_features)``.
        verbose: Emit logging warnings.
        n_lambdas: Length of the regularisation path.
        param: eBIC dimensionality penalty weight (gamma).
        feature_names: Override column names.
        n_levels_threshold: Ordinal detection threshold.

    Returns:
        :class:`MixedGraphResult`
    """
    return _legacy_fit(
        data,
        verbose=verbose,
        n_lambdas=n_lambdas,
        param=param,
        feature_names=feature_names,
        n_levels_threshold=n_levels_threshold,
    )

Estimate a mixed graphical model (backward-compatible wrapper).

.. deprecated:: Prefer :class:MixedGraphicalLasso directly.

Parameters:

Name	Type	Description	Default
`data`	`DataFrame \| NDArray[Any]`	Data matrix `(n_samples, n_features)`.	required
`verbose`	`bool`	Emit logging warnings.	`True`
`n_lambdas`	`int`	Length of the regularisation path.	`50`
`param`	`float`	eBIC dimensionality penalty weight (gamma).	`0.1`
`feature_names`	`list[str] \| None`	Override column names.	`None`
`n_levels_threshold`	`int`	Ordinal detection threshold.	`20`

Returns:

Type	Description
`MixedGraphResult`	class:`MixedGraphResult`

Source code in humlpy/estimation.py

def mixed_graph_gauss(
    data: pd.DataFrame | NDArray[Any],
    *,
    verbose: bool = True,
    n_lambdas: int = 50,
    param: float = 0.1,
    feature_names: list[str] | None = None,
    n_levels_threshold: int = 20,
) -> MixedGraphResult:
    """Estimate a mixed graphical model (backward-compatible wrapper).

    .. deprecated::
        Prefer :class:`MixedGraphicalLasso` directly.

    Args:
        data: Data matrix ``(n_samples, n_features)``.
        verbose: Emit logging warnings.
        n_lambdas: Length of the regularisation path.
        param: eBIC dimensionality penalty weight (gamma).
        feature_names: Override column names.
        n_levels_threshold: Ordinal detection threshold.

    Returns:
        :class:`MixedGraphResult`
    """
    return _legacy_fit(
        data,
        verbose=verbose,
        n_lambdas=n_lambdas,
        param=param,
        feature_names=feature_names,
        n_levels_threshold=n_levels_threshold,
    )

Legacy result container kept for backward compatibility.

New code should use :class:MixedGraphicalLasso directly and access results via its attributes.

Attributes:

Name	Type	Description
`precision_matrix`	`NDArray[Any]`	Estimated partial correlation matrix as a numpy array.
`adjacency_matrix`	`NDArray[bool_]`	Boolean adjacency matrix (True = edge present).
`correlation_matrix`	`NDArray[Any]`	Sample correlation matrix fed to graphical lasso.
`n_edges`	`int`	Number of edges in the estimated graph.
`max_degree`	`int`	Largest node degree.
`initial_mat_singular`	`bool`	Whether the correlation matrix needed PD projection.
`feature_names`	`list[str] \| None`	Variable names (None if not provided).

Source code in humlpy/estimation.py

@dataclass
class MixedGraphResult:
    """Legacy result container kept for backward compatibility.

    New code should use :class:`MixedGraphicalLasso` directly and access
    results via its attributes.

    Attributes:
        precision_matrix: Estimated partial correlation matrix as a numpy array.
        adjacency_matrix: Boolean adjacency matrix (True = edge present).
        correlation_matrix: Sample correlation matrix fed to graphical lasso.
        n_edges: Number of edges in the estimated graph.
        max_degree: Largest node degree.
        initial_mat_singular: Whether the correlation matrix needed PD projection.
        feature_names: Variable names (None if not provided).
    """

    precision_matrix: NDArray[Any]
    adjacency_matrix: NDArray[np.bool_]
    correlation_matrix: NDArray[Any]
    n_edges: int
    max_degree: int
    initial_mat_singular: bool
    feature_names: list[str] | None = None

    def __repr__(self) -> str:
        """Return string representation."""
        n_nodes = self.precision_matrix.shape[0]
        return f"MixedGraphResult(n_nodes={n_nodes}, n_edges={self.n_edges}, max_degree={self.max_degree})"

`repr()`

Return string representation.

Source code in humlpy/estimation.py

def __repr__(self) -> str:
    """Return string representation."""
    n_nodes = self.precision_matrix.shape[0]
    return f"MixedGraphResult(n_nodes={n_nodes}, n_edges={self.n_edges}, max_degree={self.max_degree})"

Calculate number of edges given a precision matrix.

Thin wrapper around :func:_edgenumber kept for backward compatibility.

Parameters:

Name	Type	Description	Default
`precision`	`NDArray[Any]`	Square precision or partial-correlation matrix.	required
`cut`	`float`	Entries with \|value\| <= cut are treated as zero.	`0.0`

Returns:

Type	Description
`int`	Number of edges (lower-triangular non-zeros).

Source code in humlpy/estimation.py

def edgenumber(precision: NDArray[Any], *, cut: float = 0.0) -> int:
    """Calculate number of edges given a precision matrix.

    Thin wrapper around :func:`_edgenumber` kept for backward compatibility.

    Args:
        precision: Square precision or partial-correlation matrix.
        cut: Entries with |value| <= *cut* are treated as zero.

    Returns:
        Number of edges (lower-triangular non-zeros).
    """
    return _edgenumber(precision, cut=cut)

API Reference

Main Classes

__init__(n_levels_threshold=20)

fit(x)

n_edges_ property

__init__(n_lambdas=50, ebic_gamma=0.1, n_levels_threshold=20)

fit(x)

Model Selection

Graph Classes

adjacency_matrix property

causal_order property

edges property

nodes property

num_edges property

num_nodes property

__init__(nodes=None, edges=None)

copy()

from_pandas_adjacency(pd_amat) classmethod

is_adjacent(i, j)

is_clique(potential_clique)

neighbors(node)

remove_edge(i, j)

remove_node(node)

show()

to_networkx()

Correlation Functions

__init__(max_cor=0.9999, solver='brent', max_iter=100, tol=1e-10)

fit(x, y)

__init__(max_cor=0.9999, n_levels_threshold=20)

fit(x, y)

Legacy Functions

__repr__()

`init(n_levels_threshold=20)`

`fit(x)`

`n_edges_` `property`

`init(n_lambdas=50, ebic_gamma=0.1, n_levels_threshold=20)`

`fit(x)`

`adjacency_matrix` `property`

`causal_order` `property`

`edges` `property`

`nodes` `property`

`num_edges` `property`

`num_nodes` `property`

`init(nodes=None, edges=None)`

`copy()`

`from_pandas_adjacency(pd_amat)` `classmethod`

`is_adjacent(i, j)`

`is_clique(potential_clique)`

`neighbors(node)`

`remove_edge(i, j)`

`remove_node(node)`

`show()`

`to_networkx()`

`init(max_cor=0.9999, solver='brent', max_iter=100, tol=1e-10)`

`fit(x, y)`

`init(max_cor=0.9999, n_levels_threshold=20)`

`fit(x, y)`

`repr()`