quanttoolbox.stats

stats.distributions

Python alternatives

Simple wrappers (normal/t/chi2/F/MVN/beta/skew-normal) are one-liners around scipy.stats — switch/keep as call-site convenience, no reason to avoid calling scipy.stats directly either. Skew-normal in particular fully switches to scipy.stats.skewnorm (verified numerically equivalent — no hand-rolled Newton iteration needed for its quantile function). GQF1/GQF2, skew-t (Azzalini's, not scipy's Jones-Faddy jf_skew_t), the Bates distribution, and the normal-ratio (Hinkley) distribution are genuinely niche — nothing in scipy, statsmodels, or elsewhere implements these. Keep. Poisson-binomial matches scipy.stats.poisson_binom exactly — switch. Order-statistic CDF is a thin scipy.stats.binom.sf-backed convenience; its quantile counterpart (grid search) has no scipy equivalent — hybrid.

quanttoolbox.stats.distributions

Probability distribution functions: CDF/PDF/quantile wrappers and the GQF (generalized quadratic form) distribution family.

Ported from QuantToolBox/stats/{cdfn,cdfni,cdft,cdfti,cdftc,cdfchi2, cdfchi2c,cdff,cdffc,cdfmvn,pdfmvn,pdfn,rndmvn,gqf1_,gqf2_}.m, extended with HSF toolbox stats/{cdfSN,cdfSNi,pdfSN,momSN,rndSN,cdfST,cdfSTi,pdfST, momST,rndST,cdfBates,pdfBates,cdfbeta,pdfbeta,cdfig,pdfig,cdfln,pdfln, cdfNormalRatio,pdfNormalRatio,pdfPoissonBinomial,cdfchi2i,pdft, compute_cdf_order_statistics,compute_inv_cdf_order_statistics, constant_correlation_matrix}.m (HSF toolbox port -- see docs/migration_map.md).

Translation notes:

• The simple distribution wrappers (normal/Student-t/chi-square/F/MVN) are thin passthroughs to scipy.stats -- MATLAB's Statistics Toolbox functions (tcdf, chi2cdf, mvncdf, ...) map directly onto scipy.stats.{t,chi2,f,multivariate_normal}, so there's no need to hand-roll any of this. student_t_pdf/chi2_ppf (from the HSF toolbox's pdft.m/cdfchi2i.m) fill out this same family.
• GQF #1 is the distribution of sum(a_i * (Z_i + b_i)^2) for independent standard normals Z_i (a "weighted noncentral chi-square mixture"); GQF #2 is the distribution of the general quadratic form (X - mu)' Q (X - mu) (or similar) for X ~ N(mu, Sigma). Both use cumulant-based approximations (a Laguerre-type series for GQF1, a three-cumulant noncentral-chi-square matching for GQF2 -- see Solomon & Stephens (1978) / Buckley & Eagleson (1988) for the underlying theory). These have no scipy equivalent and are ported algorithm-for-algorithm.
• Beta/lognormal/inverse-Gaussian (cdfbeta/pdfbeta, cdfln/ pdfln, cdfig/pdfig): thin wrappers with the original's call signature, same "call-site compatibility" reasoning as the simple wrappers above -- beta_cdf/beta_pdf go straight to scipy.stats.beta; lognormal_cdf/lognormal_pdf and inverse_gaussian_cdf/inverse_gaussian_pdf are evaluated directly from the closed-form formulas the originals use (verified algebraically equivalent to scipy.stats.lognorm/invgauss under those distributions' own shape/scale reparameterization, but ported directly rather than requiring callers to compute that reparameterization themselves).
• Bates distribution (bates_cdf/bates_pdf): no scipy equivalent (scipy has no named Bates distribution) -- ported algorithm-for-algorithm (scipy.special.comb replaces MATLAB's nchoosek).
• Poisson-binomial PMF (poisson_binomial_pmf): matches scipy.stats.poisson_binom exactly (verified numerically against both the original's FFT and direct-recursion branches) -- switched to the scipy-backed distribution rather than hand-rolling either branch.
• Order statistics (order_statistic_cdf/order_statistic_ppf): the order-statistic CDF formula F_{i:n}(x) = P(Binom(n, F_x) >= i) is evaluated via scipy.stats.binom.sf rather than hand-computing binomial coefficients in a loop -- same formula, more numerically stable for large n. The quantile version is a grid search (given sample points x and their CDF values F_x, find the first point where the order-statistic CDF crosses alpha) with no scipy equivalent -- ported directly.
• Normal-ratio distribution (normal_ratio_cdf/normal_ratio_pdf, Hinkley's distribution for the ratio of two independent normals): no scipy equivalent -- ported directly, using this module's own quanttoolbox.stats.multivariate.bvn_cdf in place of the original's cdfbvn calls, matching the original's (p, a_z, b_z, c, rho_z) multi-output signature.
• Skew-normal (skew_normal_{cdf,ppf,pdf,moments,rvs}): switched entirely to scipy.stats.skewnorm -- verified numerically that pdfSN.m's formula, momSN.m's moment formulas, and scipy.stats.skewnorm's pdf/cdf/stats(moments="mvsk") all agree exactly under the direct parameter mapping (eta -> a, xi -> loc, omega -> scale). The original's two alternate numerical CDF branches (mtd 0 vs 1, both bivariate-normal-orthant identities for the same value) and its hand-rolled Newton-iteration quantile function (cdfSNi.m) are therefore unnecessary -- scipy computes the exact analytic answer directly via .cdf/.ppf.
• Skew-t (skew_t_{cdf,ppf,pdf,moments,rvs}): not the same family as scipy.stats.jf_skew_t (Jones & Faddy's skew-t, a different two-shape- parameter family) -- this is Azzalini's skew-t (one skewness parameter eta plus degrees of freedom nu), which has no scipy equivalent. Ported algorithm-for-algorithm, including the original's Newton-iteration quantile function (cdfSTi.m, via quanttoolbox.config.NewtonConfig) and its bivariate-Student-t-based CDF (via quanttoolbox.stats.multivariate.bvt_cdf).
• max_size.m (a MATLAB shape-broadcasting helper) is not ported as a function -- superseded entirely by numpy's native broadcasting (np.broadcast_arrays), used throughout this module and stats.multivariate instead.
• MATLAB's missex(cdf, cond) (replace-with-missing where cond is true) is inlined as np.where(cond, np.nan, cdf).
• MATLAB's eig is replaced with numpy.linalg.eigh in gqf2_to_gqf1 since the matrix involved is symmetric by construction; the eigenvalue order may differ from MATLAB's, but the resulting (a, b) pairs describe the same distribution regardless of order.

NormalRatioResult(p, a_z, b_z, c, rho_z) dataclass

CDF or PDF of Z = X / Y (X, Y independent normals), plus the intermediate quantities (a_z, b_z, c, rho_z) the original returns alongside the probability.

bates_cdf(x, n)

CDF of the Bates distribution: the mean of n iid Uniform(0, 1) variables. No scipy equivalent.

Original: stats/cdfBates.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def bates_cdf(x: np.ndarray, n: int) -> np.ndarray:
    """CDF of the Bates distribution: the mean of `n` iid Uniform(0, 1)
    variables. No scipy equivalent.

    Original: stats/cdfBates.m (HSF toolbox)
    """
    x = np.asarray(x, dtype=float)
    cdf = np.zeros_like(x)
    for k in range(n + 1):
        s = (n * x > k).astype(float)
        cdf = cdf + ((-1.0) ** k) * comb(n, k) * (n * x - k) ** n * s
    cdf = cdf / factorial(n)
    cdf = np.where(x == 0, 0.0, cdf)
    cdf = np.where(x == 1, 1.0, cdf)
    return cdf

bates_pdf(x, n)

PDF of the Bates distribution: the mean of n iid Uniform(0, 1) variables. No scipy equivalent.

Original: stats/pdfBates.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def bates_pdf(x: np.ndarray, n: int) -> np.ndarray:
    """PDF of the Bates distribution: the mean of `n` iid Uniform(0, 1)
    variables. No scipy equivalent.

    Original: stats/pdfBates.m (HSF toolbox)
    """
    x = np.asarray(x, dtype=float)
    pdf = np.zeros_like(x)
    for k in range(n + 1):
        pdf = pdf + ((-1.0) ** k) * comb(n, k) * (n * x - k) ** (n - 1) * np.sign(n * x - k)
    pdf = 0.5 * pdf * n / factorial(n - 1)
    if n > 1:
        pdf = np.where(x == 0, 0.0, pdf)
        pdf = np.where(x == 1, 0.0, pdf)
    else:
        pdf = np.where(x == 0, 1.0, pdf)
        pdf = np.where(x == 1, 1.0, pdf)
    return pdf

beta_cdf(x, a, b)

Original: stats/cdfbeta.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def beta_cdf(x: np.ndarray, a: float, b: float) -> np.ndarray:
    """Original: stats/cdfbeta.m (HSF toolbox)"""
    return beta_dist.cdf(x, a, b)

beta_pdf(x, a, b)

Original: stats/pdfbeta.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def beta_pdf(x: np.ndarray, a: float, b: float) -> np.ndarray:
    """Original: stats/pdfbeta.m (HSF toolbox)"""
    return beta_dist.pdf(x, a, b)

chi2_cdf(x, nu)

Original: stats/cdfchi2.m

Source code in src/quanttoolbox/stats/distributions.py

def chi2_cdf(x: np.ndarray, nu: float) -> np.ndarray:
    """Original: stats/cdfchi2.m"""
    return chi2.cdf(x, df=nu)

chi2_ppf(p, nu)

Original: stats/cdfchi2i.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def chi2_ppf(p: np.ndarray, nu: float) -> np.ndarray:
    """Original: stats/cdfchi2i.m (HSF toolbox)"""
    return chi2.ppf(p, df=nu)

chi2_sf(x, nu)

Upper-tail (survival) chi-square. Original: stats/cdfchi2c.m

Source code in src/quanttoolbox/stats/distributions.py

def chi2_sf(x: np.ndarray, nu: float) -> np.ndarray:
    """Upper-tail (survival) chi-square. Original: stats/cdfchi2c.m"""
    return chi2.sf(x, df=nu)

constant_correlation_matrix(n, rho)

An n x n correlation matrix with constant off-diagonal correlation rho and unit diagonal.

Original: stats/constant_correlation_matrix.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def constant_correlation_matrix(n: int, rho: float) -> np.ndarray:
    """An n x n correlation matrix with constant off-diagonal correlation
    `rho` and unit diagonal.

    Original: stats/constant_correlation_matrix.m (HSF toolbox)
    """
    c = np.full((n, n), rho, dtype=float)
    np.fill_diagonal(c, 1.0)
    return c

f_cdf(x, nu1, nu2)

Original: stats/cdff.m

Source code in src/quanttoolbox/stats/distributions.py

def f_cdf(x: np.ndarray, nu1: float, nu2: float) -> np.ndarray:
    """Original: stats/cdff.m"""
    return f.cdf(x, dfn=nu1, dfd=nu2)

f_sf(x, nu1, nu2)

Upper-tail (survival) F distribution. Original: stats/cdffc.m

Source code in src/quanttoolbox/stats/distributions.py

def f_sf(x: np.ndarray, nu1: float, nu2: float) -> np.ndarray:
    """Upper-tail (survival) F distribution. Original: stats/cdffc.m"""
    return f.sf(x, dfn=nu1, dfd=nu2)

gqf1_cdf(x, a, b, beta=None, order=200)

CDF of the GQF #1 distribution. Original: stats/gqf1_cdf.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf1_cdf(
    x: np.ndarray, a: np.ndarray, b: np.ndarray, beta: float | None = None, order: int = 200
) -> np.ndarray:
    """CDF of the GQF #1 distribution. Original: stats/gqf1_cdf.m"""
    a = np.asarray(a, dtype=float)
    n = a.shape[0]
    if beta is None or beta == 0:
        beta = 0.8 * a.min()

    coeffs = gqf1_coeffs(a, b, beta, order)
    x_star = np.asarray(x, dtype=float) / beta

    cdf = np.zeros_like(x_star, dtype=float)
    for j in range(order + 1):
        cdf = cdf + coeffs[j] * chi2.cdf(x_star, df=n + 2 * j)

    return np.where((cdf < -1e-4) | (cdf > 1.0001), np.nan, cdf)

gqf1_coeffs(a, b, beta, order)

Laguerre-series mixing coefficients for the GQF #1 CDF/PDF.

Original: stats/gqf1_coeffs.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf1_coeffs(a: np.ndarray, b: np.ndarray, beta: float, order: int) -> np.ndarray:
    """Laguerre-series mixing coefficients for the GQF #1 CDF/PDF.

    Original: stats/gqf1_coeffs.m
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    n = max(a.shape[0], b.shape[0])
    a = np.broadcast_to(a, (n,))
    b = np.broadcast_to(b, (n,))

    beta_a = beta / a
    beta_a_c = 1 - beta_a
    b2_a = b**2 / a
    zeta = np.sum(b**2)

    g = np.zeros(order)
    for m in range(1, order + 1):
        g[m - 1] = np.sum(beta_a_c**m) + m * beta * np.sum(b2_a * beta_a_c ** (m - 1))

    coeffs = np.zeros(order + 1)
    coeffs[0] = np.exp(-zeta / 2) * np.prod(np.sqrt(beta_a))
    for j in range(1, order + 1):
        coeffs[j] = np.sum(g[0:j][::-1] * coeffs[0:j]) / (2 * j)

    return coeffs

gqf1_moments(a, b)

Mean, std dev, skewness, and excess kurtosis of the GQF #1 distribution.

Original: stats/gqf1_moments.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf1_moments(a: np.ndarray, b: np.ndarray) -> tuple[float, float, float, float]:
    """Mean, std dev, skewness, and excess kurtosis of the GQF #1 distribution.

    Original: stats/gqf1_moments.m
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    b2, a2, a3, a4 = b**2, a**2, a**3, a**4

    mean = np.sum(a * (1 + b2))
    denom = np.sum(a2 * (1 + 2 * b2))
    sigma = np.sqrt(2 * denom)
    gamma1 = 2 * np.sqrt(2) * np.sum(a3 * (1 + 3 * b2)) / denom**1.5
    gamma2 = 12 * np.sum(a4 * (1 + 4 * b2)) / denom**2
    return mean, sigma, gamma1, gamma2

gqf1_pdf(x, a, b, beta=None, order=100)

PDF of the GQF #1 distribution. Original: stats/gqf1_pdf.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf1_pdf(
    x: np.ndarray, a: np.ndarray, b: np.ndarray, beta: float | None = None, order: int = 100
) -> np.ndarray:
    """PDF of the GQF #1 distribution. Original: stats/gqf1_pdf.m"""
    a = np.asarray(a, dtype=float)
    n = a.shape[0]
    if beta is None or beta == 0:
        beta = 0.8 * a.min()

    coeffs = gqf1_coeffs(a, b, beta, order)
    x_star = np.asarray(x, dtype=float) / beta

    pdf = np.zeros_like(x_star, dtype=float)
    for j in range(order + 1):
        pdf = pdf + coeffs[j] * chi2.pdf(x_star, df=n + 2 * j) / beta

    return pdf

gqf1_to_gqf2(a, b)

Convert GQF #1 parameters (a, b) to GQF #2 parameters (mu, Sigma, Q).

Original: stats/gqf1_to_gqf2.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf1_to_gqf2(a: np.ndarray, b: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Convert GQF #1 parameters (a, b) to GQF #2 parameters (mu, Sigma, Q).

    Original: stats/gqf1_to_gqf2.m
    """
    a = np.asarray(a, dtype=float)
    b = np.asarray(b, dtype=float)
    n = max(a.shape[0], b.shape[0])
    a = np.broadcast_to(a, (n,))
    b = np.broadcast_to(b, (n,))

    mu = b
    sigma = np.eye(n)
    q = np.diag(a)
    return mu, sigma, q

gqf2_cdf(x, mu, sigma, q)

CDF of the GQF #2 distribution, via noncentral-chi-square matching.

Original: stats/gqf2_cdf.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf2_cdf(x: np.ndarray, mu: np.ndarray, sigma: np.ndarray, q: np.ndarray) -> np.ndarray:
    """CDF of the GQF #2 distribution, via noncentral-chi-square matching.

    Original: stats/gqf2_cdf.m
    """
    m, sig, _, _, s1, s2 = gqf2_moments(mu, sigma, q)
    nu, zeta = _gqf2_noncentral_chi2_params(s1, s2)

    mu_star = nu + zeta
    sigma_star = np.sqrt(2 * nu + 4 * zeta)

    x_star = (np.asarray(x, dtype=float) - m) / sig
    x_star = mu_star + sigma_star * x_star
    return ncx2.cdf(x_star, df=nu, nc=zeta)

gqf2_moments(mu, sigma, q)

Mean, std dev, skewness, excess kurtosis, and the two Pearson-type shape statistics (s1, s2) used by gqf2_cdf/gqf2_pdf.

Original: stats/gqf2_moments.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf2_moments(
    mu: np.ndarray, sigma: np.ndarray, q: np.ndarray
) -> tuple[float, float, float, float, float, float]:
    """Mean, std dev, skewness, excess kurtosis, and the two Pearson-type
    shape statistics (s1, s2) used by gqf2_cdf/gqf2_pdf.

    Original: stats/gqf2_moments.m
    """
    mu = np.asarray(mu, dtype=float).flatten()
    sigma = np.asarray(sigma, dtype=float)
    q = np.asarray(q, dtype=float)
    n = sigma.shape[0]
    n_kappa = 4

    kappa = np.zeros(n_kappa)
    p = q @ sigma
    p_old = np.eye(n)
    for k in range(1, n_kappa + 1):
        p_new = p_old @ p
        kappa[k - 1] = (
            (2 ** (k - 1)) * factorial(k - 1) * (np.trace(p_new) + k * mu.T @ p_old @ q @ mu)
        )
        p_old = p_new

    mean = kappa[0]
    sigma_out = np.sqrt(kappa[1])
    gamma1 = kappa[2] / sigma_out**3
    gamma2 = kappa[3] / sigma_out**4

    s1 = gamma1 / np.sqrt(8)
    s2 = gamma2 / 12
    return mean, sigma_out, gamma1, gamma2, s1, s2

gqf2_pdf(x, mu, sigma, q)

PDF of the GQF #2 distribution, via noncentral-chi-square matching.

Original: stats/gqf2_pdf.m

Source code in src/quanttoolbox/stats/distributions.py

def gqf2_pdf(x: np.ndarray, mu: np.ndarray, sigma: np.ndarray, q: np.ndarray) -> np.ndarray:
    """PDF of the GQF #2 distribution, via noncentral-chi-square matching.

    Original: stats/gqf2_pdf.m
    """
    m, sig, _, _, s1, s2 = gqf2_moments(mu, sigma, q)
    nu, zeta = _gqf2_noncentral_chi2_params(s1, s2)

    mu_star = nu + zeta
    sigma_star = np.sqrt(2 * nu + 4 * zeta)

    x_star = (np.asarray(x, dtype=float) - m) / sig
    x_star = mu_star + sigma_star * x_star
    return ncx2.pdf(x_star, df=nu, nc=zeta) * (sigma_star / sig)

gqf2_to_gqf1(mu, sigma, q)

Convert GQF #2 parameters (mu, Sigma, Q) to GQF #1 parameters (a, b).

Original: stats/gqf2_to_gqf1.m

Note: eigenvalue ordering may differ from MATLAB's eig (which does not sort), but the resulting (a, b) pairs describe the same distribution regardless of order.

Source code in src/quanttoolbox/stats/distributions.py

def gqf2_to_gqf1(mu: np.ndarray, sigma: np.ndarray, q: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Convert GQF #2 parameters (mu, Sigma, Q) to GQF #1 parameters (a, b).

    Original: stats/gqf2_to_gqf1.m

    Note: eigenvalue ordering may differ from MATLAB's ``eig`` (which does
    not sort), but the resulting (a, b) pairs describe the same
    distribution regardless of order.
    """
    mu = np.asarray(mu, dtype=float).flatten()
    sigma = np.asarray(sigma, dtype=float)
    q = np.asarray(q, dtype=float)

    sigma_sqrt = linalg.sqrtm(sigma).real
    b_mat = sigma_sqrt @ q @ sigma_sqrt
    eigenvalues, eigenvectors = np.linalg.eigh(b_mat)

    m_vec = eigenvectors.T @ np.linalg.inv(sigma_sqrt) @ mu

    a = eigenvalues
    b = m_vec
    return a, b

inverse_gaussian_cdf(x, mu, lam)

CDF of the inverse Gaussian (Wald) distribution, mean mu, shape lam.

Original: stats/cdfig.m (HSF toolbox; algebraically equivalent to scipy.stats.invgauss.cdf(x, mu / lam, scale=lam), ported directly to keep the original's (mu, lam) call signature)

Source code in src/quanttoolbox/stats/distributions.py

def inverse_gaussian_cdf(x: np.ndarray, mu: float, lam: float) -> np.ndarray:
    """CDF of the inverse Gaussian (Wald) distribution, mean `mu`, shape
    `lam`.

    Original: stats/cdfig.m (HSF toolbox; algebraically equivalent to
    ``scipy.stats.invgauss.cdf(x, mu / lam, scale=lam)``, ported directly
    to keep the original's ``(mu, lam)`` call signature)
    """
    x = np.asarray(x, dtype=float)
    alpha = (x - mu) / mu * np.sqrt(lam / x)
    alpha_bar = -(x + mu) / mu * np.sqrt(lam / x)
    return normal_cdf(alpha) + np.exp(2 * lam / mu) * normal_cdf(alpha_bar)

inverse_gaussian_pdf(x, mu, lam)

PDF of the inverse Gaussian (Wald) distribution, mean mu, shape lam.

Original: stats/pdfig.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def inverse_gaussian_pdf(x: np.ndarray, mu: float, lam: float) -> np.ndarray:
    """PDF of the inverse Gaussian (Wald) distribution, mean `mu`, shape
    `lam`.

    Original: stats/pdfig.m (HSF toolbox)
    """
    x = np.asarray(x, dtype=float)
    return np.sqrt(lam / (2 * np.pi * x**3)) * np.exp(-0.5 * lam / mu**2 * (x - mu) ** 2 / x)

lognormal_cdf(x, mu, sigma)

CDF of a lognormal variable, i.e. log(X) ~ N(mu, sigma^2).

Original: stats/cdfln.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def lognormal_cdf(x: np.ndarray, mu: float, sigma: float) -> np.ndarray:
    """CDF of a lognormal variable, i.e. log(X) ~ N(mu, sigma^2).

    Original: stats/cdfln.m (HSF toolbox)
    """
    x = np.asarray(x, dtype=float)
    return np.asarray(normal_cdf((np.log(x) - mu) / sigma))

lognormal_pdf(x, mu, sigma)

PDF of a lognormal variable, i.e. log(X) ~ N(mu, sigma^2).

Original: stats/pdfln.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def lognormal_pdf(x: np.ndarray, mu: float, sigma: float) -> np.ndarray:
    """PDF of a lognormal variable, i.e. log(X) ~ N(mu, sigma^2).

    Original: stats/pdfln.m (HSF toolbox)
    """
    x = np.asarray(x, dtype=float)
    y = (np.log(x) - mu) / sigma
    return 1.0 / (x * np.sqrt(2 * np.pi) * sigma) * np.exp(-0.5 * y**2)

mvn_cdf(x, mu, sigma)

CDF of the multivariate normal N(mu, Sigma). Original: stats/cdfmvn.m

Source code in src/quanttoolbox/stats/distributions.py

def mvn_cdf(x: np.ndarray, mu: np.ndarray, sigma: np.ndarray) -> np.ndarray:
    """CDF of the multivariate normal N(mu, Sigma). Original: stats/cdfmvn.m"""
    return multivariate_normal.cdf(x, mean=np.asarray(mu).flatten(), cov=sigma)

mvn_pdf(x, mu, sigma)

PDF of the multivariate normal N(mu, Sigma). Original: stats/pdfmvn.m

Source code in src/quanttoolbox/stats/distributions.py

def mvn_pdf(x: np.ndarray, mu: np.ndarray, sigma: np.ndarray) -> np.ndarray:
    """PDF of the multivariate normal N(mu, Sigma). Original: stats/pdfmvn.m"""
    return multivariate_normal.pdf(x, mean=np.asarray(mu).flatten(), cov=sigma)

mvn_rvs(mu, sigma, n_samples, random_state=None)

Draw samples from the multivariate normal N(mu, Sigma).

Original: stats/rndmvn.m

Source code in src/quanttoolbox/stats/distributions.py

def mvn_rvs(mu: np.ndarray, sigma: np.ndarray, n_samples: int, random_state=None) -> np.ndarray:
    """Draw samples from the multivariate normal N(mu, Sigma).

    Original: stats/rndmvn.m
    """
    return multivariate_normal.rvs(
        mean=np.asarray(mu).flatten(), cov=sigma, size=n_samples, random_state=random_state
    )

normal_cdf(x, mu=0.0, sigma=1.0)

Original: stats/cdfn.m

Source code in src/quanttoolbox/stats/distributions.py

def normal_cdf(x: float | np.ndarray, mu: float = 0.0, sigma: float = 1.0) -> float | np.ndarray:
    """Original: stats/cdfn.m"""
    return norm.cdf(x, loc=mu, scale=sigma)

normal_pdf(x, mu=0.0, sigma=1.0)

Original: stats/pdfn.m

Source code in src/quanttoolbox/stats/distributions.py

def normal_pdf(x: np.ndarray, mu: float = 0.0, sigma: float = 1.0) -> np.ndarray:
    """Original: stats/pdfn.m"""
    return norm.pdf(x, loc=mu, scale=sigma)

normal_ppf(p)

Quantile function of N(0,1). Original: stats/cdfni.m

Source code in src/quanttoolbox/stats/distributions.py

def normal_ppf(p: float | np.ndarray) -> float | np.ndarray:
    """Quantile function of N(0,1). Original: stats/cdfni.m"""
    return norm.ppf(p)

normal_ratio_cdf(z, mu_x, sigma_x, mu_y, sigma_y)

CDF of Z = X / Y, for independent X ~ N(mu_x, sigma_x^2) and Y ~ N(mu_y, sigma_y^2) (Hinkley's ratio distribution). No scipy equivalent.

Original: stats/cdfNormalRatio.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def normal_ratio_cdf(
    z: np.ndarray, mu_x: float, sigma_x: float, mu_y: float, sigma_y: float
) -> NormalRatioResult:
    """CDF of Z = X / Y, for independent X ~ N(mu_x, sigma_x^2) and
    Y ~ N(mu_y, sigma_y^2) (Hinkley's ratio distribution). No scipy
    equivalent.

    Original: stats/cdfNormalRatio.m (HSF toolbox)
    """
    z, a_z, b_z, c, rho_z = _normal_ratio_terms(z, mu_x, sigma_x, mu_y, sigma_y)

    x1 = (mu_x - mu_y * z) / (sigma_x * sigma_y * a_z)
    y1 = -mu_y / sigma_y
    x2, y2 = -x1, -y1

    p = bvn_cdf(x1, y1, rho_z) + bvn_cdf(x2, y2, rho_z)
    return NormalRatioResult(p=p, a_z=a_z, b_z=b_z, c=np.full_like(a_z, c), rho_z=rho_z)

normal_ratio_pdf(z, mu_x, sigma_x, mu_y, sigma_y)

PDF of Z = X / Y, for independent X ~ N(mu_x, sigma_x^2) and Y ~ N(mu_y, sigma_y^2) (Hinkley's ratio distribution). No scipy equivalent.

Original: stats/pdfNormalRatio.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def normal_ratio_pdf(
    z: np.ndarray, mu_x: float, sigma_x: float, mu_y: float, sigma_y: float
) -> NormalRatioResult:
    """PDF of Z = X / Y, for independent X ~ N(mu_x, sigma_x^2) and
    Y ~ N(mu_y, sigma_y^2) (Hinkley's ratio distribution). No scipy
    equivalent.

    Original: stats/pdfNormalRatio.m (HSF toolbox)
    """
    z, a_z, b_z, c, rho_z = _normal_ratio_terms(z, mu_x, sigma_x, mu_y, sigma_y)
    a2_z, a3_z, b2_z = a_z**2, a_z**3, b_z**2

    p1 = b_z / (sigma_x * sigma_y * np.sqrt(2 * np.pi) * a3_z)
    p2 = normal_cdf(b_z / a_z) - normal_cdf(-b_z / a_z)
    p3 = np.exp((b2_z - c * a2_z) / (2 * a2_z))
    p4 = np.exp(-c / 2) / (sigma_x * sigma_y * a2_z * np.pi)
    p = p1 * p2 * p3 + p4
    return NormalRatioResult(p=p, a_z=a_z, b_z=b_z, c=np.full_like(a_z, c), rho_z=rho_z)

order_statistic_cdf(f_x, n, i_select=None)

CDF of the i-th order statistic (i = 1..n) of n iid draws from a distribution with CDF value(s) f_x: F_{i:n} = P(Binom(n, F_x) >= i), evaluated via scipy.stats.binom.sf rather than hand-computing binomial coefficients (same formula, more numerically stable for large n).

f_x may be an array (one row per evaluation point); returns an array of shape (len(f_x), len(i_select)) with one column per selected order statistic i (default: all of 1..n).

Original: stats/compute_cdf_order_statistics.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def order_statistic_cdf(f_x: np.ndarray, n: int, i_select: np.ndarray | None = None) -> np.ndarray:
    """CDF of the i-th order statistic (i = 1..n) of n iid draws from a
    distribution with CDF value(s) `f_x`: ``F_{i:n} = P(Binom(n, F_x) >=
    i)``, evaluated via ``scipy.stats.binom.sf`` rather than hand-computing
    binomial coefficients (same formula, more numerically stable for large
    n).

    `f_x` may be an array (one row per evaluation point); returns an array
    of shape (len(f_x), len(i_select)) with one column per selected order
    statistic i (default: all of 1..n).

    Original: stats/compute_cdf_order_statistics.m (HSF toolbox)
    """
    f_x = np.atleast_1d(np.asarray(f_x, dtype=float))
    if i_select is None:
        i_select = np.arange(1, n + 1)
    i_select = np.asarray(i_select, dtype=int)

    # F_{i:n}(x) = P(Binom(n, F_x) >= i) = sf(i - 1, n, F_x)
    out = np.stack([binom.sf(i - 1, n, f_x) for i in i_select], axis=-1)
    return out

order_statistic_ppf(alpha, x, f_x, n, i_select=None)

Quantile of the i-th order statistic, found by grid search: given sample points x with CDF values f_x, returns (for each alpha and each selected order statistic i) the first x value where order_statistic_cdf crosses alpha. No scipy equivalent -- ported directly.

Original: stats/compute_inv_cdf_order_statistics.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def order_statistic_ppf(
    alpha: np.ndarray,
    x: np.ndarray,
    f_x: np.ndarray,
    n: int,
    i_select: np.ndarray | None = None,
) -> np.ndarray:
    """Quantile of the i-th order statistic, found by grid search: given
    sample points `x` with CDF values `f_x`, returns (for each `alpha` and
    each selected order statistic i) the first `x` value where
    `order_statistic_cdf` crosses `alpha`. No scipy equivalent -- ported
    directly.

    Original: stats/compute_inv_cdf_order_statistics.m (HSF toolbox)
    """
    alpha = np.atleast_1d(np.asarray(alpha, dtype=float))
    x = np.asarray(x, dtype=float)
    if i_select is None:
        i_select = np.arange(1, n + 1)
    i_select = np.asarray(i_select, dtype=int)

    f_i_n = order_statistic_cdf(f_x, n, i_select)  # shape (len(x), len(i_select))

    q = np.full((alpha.shape[0], i_select.shape[0]), np.nan)
    for a_idx, a in enumerate(alpha):
        for k in range(i_select.shape[0]):
            candidates = np.flatnonzero(f_i_n[:, k] >= a)
            if candidates.size > 0:
                q[a_idx, k] = x[candidates[0]]
    return q

poisson_binomial_pmf(p)

PMF of the sum of n independent Bernoulli(p_i) variables, for support k = 0..n. Matches scipy.stats.poisson_binom exactly (verified numerically against the original's own FFT and direct-recursion branches) -- neither branch is hand-rolled here.

Original: stats/pdfPoissonBinomial.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def poisson_binomial_pmf(p: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """PMF of the sum of `n` independent Bernoulli(p_i) variables, for
    support k = 0..n. Matches ``scipy.stats.poisson_binom`` exactly
    (verified numerically against the original's own FFT and
    direct-recursion branches) -- neither branch is hand-rolled here.

    Original: stats/pdfPoissonBinomial.m (HSF toolbox)
    """
    p = np.asarray(p, dtype=float).flatten()
    n = p.shape[0]
    k = np.arange(n + 1)
    return k, poisson_binom.pmf(k, p)

poisson_binomial_pmf_brute_force(p)

PMF of the sum of n independent Bernoulli(p_i) variables, via exact brute-force enumeration over all 2^n outcomes -- feasible only for small n; intended purely as an independent cross-check of poisson_binomial_pmf_dp/poisson_binomial_pmf_dft, not for general use.

Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive implements this. Promoted from HSF-Notebooks chapter 16d as an independent reference implementation, kept deliberately separate from poisson_binomial_pmf -- see that notebook's three-way cross-check of DP/DFT/brute-force methods against the packaged function.

Source code in src/quanttoolbox/stats/distributions.py

def poisson_binomial_pmf_brute_force(p: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """PMF of the sum of `n` independent Bernoulli(p_i) variables, via
    exact brute-force enumeration over all 2^n outcomes -- feasible only
    for small `n`; intended purely as an independent cross-check of
    `poisson_binomial_pmf_dp`/`poisson_binomial_pmf_dft`, not for general
    use.

    Not ported from the MATLAB HSF toolbox -- no `.m` file in `hfs-archive`
    implements this. Promoted from HSF-Notebooks chapter 16d as an
    independent reference implementation, kept deliberately separate from
    `poisson_binomial_pmf` -- see that notebook's three-way cross-check of
    DP/DFT/brute-force methods against the packaged function.
    """
    from itertools import product

    p = np.asarray(p, dtype=float)
    n = len(p)
    pmf = np.zeros(n + 1)
    for outcome_tuple in product([0, 1], repeat=n):
        outcome = np.array(outcome_tuple)
        prob = np.prod(np.where(outcome == 1, p, 1 - p))
        pmf[outcome.sum()] += prob
    return np.arange(n + 1), pmf

poisson_binomial_pmf_dft(p)

PMF of the sum of n independent Bernoulli(p_i) variables, via a DFT of the characteristic function (Fernandez & Williams (2010) / the poibin package's "DFT-CF" method).

Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive implements this. Promoted from HSF-Notebooks chapter 16d as an independent reference implementation, kept deliberately separate from poisson_binomial_pmf -- see that notebook's three-way cross-check of DP/DFT/brute-force methods against the packaged function.

Source code in src/quanttoolbox/stats/distributions.py

def poisson_binomial_pmf_dft(p: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """PMF of the sum of `n` independent Bernoulli(p_i) variables, via a
    DFT of the characteristic function (Fernandez & Williams (2010) /
    the `poibin` package's "DFT-CF" method).

    Not ported from the MATLAB HSF toolbox -- no `.m` file in `hfs-archive`
    implements this. Promoted from HSF-Notebooks chapter 16d as an
    independent reference implementation, kept deliberately separate from
    `poisson_binomial_pmf` -- see that notebook's three-way cross-check of
    DP/DFT/brute-force methods against the packaged function.
    """
    p = np.asarray(p, dtype=float)
    n = len(p)
    length = n + 1
    ell = np.arange(length)
    omega = 2 * np.pi * ell / length
    exp_iw = np.exp(1j * omega)
    chi = np.prod(1 - p[:, None] + p[:, None] * exp_iw[None, :], axis=0)
    # np.fft.fft(chi)/length matches the DFT-CF inversion formula's sign
    # convention (pmf[k] = (1/length) * sum_l chi(omega_l) * exp(-i*omega_l*k));
    # np.fft.ifft uses the opposite sign and gives the *reversed* PMF instead --
    # verified against brute-force enumeration.
    pmf = np.real(np.fft.fft(chi)) / length
    k = np.arange(length)
    return k, pmf

poisson_binomial_pmf_dp(p)

PMF of the sum of n independent Bernoulli(p_i) variables, via the direct O(n^2) DP/convolution recursion (each Bernoulli in turn either shifts the running PMF by one and scales by p_i, or leaves it in place scaled by 1 - p_i).

Not ported from the MATLAB HSF toolbox -- no .m file in hfs-archive implements this. Promoted from HSF-Notebooks chapter 16d as an independent reference implementation, kept deliberately separate from poisson_binomial_pmf -- see that notebook's three-way cross-check of DP/DFT/brute-force methods against the packaged function.

Source code in src/quanttoolbox/stats/distributions.py

def poisson_binomial_pmf_dp(p: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """PMF of the sum of `n` independent Bernoulli(p_i) variables, via the
    direct O(n^2) DP/convolution recursion (each Bernoulli in turn either
    shifts the running PMF by one and scales by `p_i`, or leaves it in
    place scaled by `1 - p_i`).

    Not ported from the MATLAB HSF toolbox -- no `.m` file in `hfs-archive`
    implements this. Promoted from HSF-Notebooks chapter 16d as an
    independent reference implementation, kept deliberately separate from
    `poisson_binomial_pmf` -- see that notebook's three-way cross-check of
    DP/DFT/brute-force methods against the packaged function.
    """
    p = np.asarray(p, dtype=float)
    pmf = np.array([1.0])
    for pi in p:
        pmf = np.concatenate([pmf, [0.0]]) * (1 - pi) + np.concatenate([[0.0], pmf]) * pi
    k = np.arange(len(pmf))
    return k, pmf

skew_normal_cdf(x, xi, omega, eta)

CDF of Azzalini's skew-normal distribution (location xi, scale omega, shape eta).

Original: stats/cdfSN.m (HSF toolbox; both of the original's numerical branches are exact identities for this same value -- scipy computes it directly, see module docstring)

Source code in src/quanttoolbox/stats/distributions.py

def skew_normal_cdf(x: np.ndarray, xi: float, omega: float, eta: float) -> np.ndarray:
    """CDF of Azzalini's skew-normal distribution (location `xi`, scale
    `omega`, shape `eta`).

    Original: stats/cdfSN.m (HSF toolbox; both of the original's numerical
    branches are exact identities for this same value -- ``scipy`` computes
    it directly, see module docstring)
    """
    return skewnorm.cdf(x, eta, loc=xi, scale=omega)

skew_normal_moments(xi, omega, eta)

Mean, std dev, skewness, and excess kurtosis of Azzalini's skew-normal distribution.

Original: stats/momSN.m (HSF toolbox; matches scipy.stats.skewnorm.stats(eta, loc=xi, scale=omega, moments="mvsk") exactly, see module docstring)

Source code in src/quanttoolbox/stats/distributions.py

def skew_normal_moments(xi: float, omega: float, eta: float) -> tuple[float, float, float, float]:
    """Mean, std dev, skewness, and excess kurtosis of Azzalini's
    skew-normal distribution.

    Original: stats/momSN.m (HSF toolbox; matches
    ``scipy.stats.skewnorm.stats(eta, loc=xi, scale=omega,
    moments="mvsk")`` exactly, see module docstring)
    """
    mean, var, skew, kurt = skewnorm.stats(eta, loc=xi, scale=omega, moments="mvsk")
    return float(mean), float(np.sqrt(var)), float(skew), float(kurt)

skew_normal_pdf(x, xi, omega, eta)

PDF of Azzalini's skew-normal distribution.

Original: stats/pdfSN.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_normal_pdf(x: np.ndarray, xi: float, omega: float, eta: float) -> np.ndarray:
    """PDF of Azzalini's skew-normal distribution.

    Original: stats/pdfSN.m (HSF toolbox)
    """
    return skewnorm.pdf(x, eta, loc=xi, scale=omega)

skew_normal_ppf(p, xi, omega, eta)

Quantile function of Azzalini's skew-normal distribution.

Original: stats/cdfSNi.m (HSF toolbox; the original's Newton iteration is unnecessary -- scipy.stats.skewnorm.ppf is exact, see module docstring)

Source code in src/quanttoolbox/stats/distributions.py

def skew_normal_ppf(p: np.ndarray, xi: float, omega: float, eta: float) -> np.ndarray:
    """Quantile function of Azzalini's skew-normal distribution.

    Original: stats/cdfSNi.m (HSF toolbox; the original's Newton iteration
    is unnecessary -- ``scipy.stats.skewnorm.ppf`` is exact, see module
    docstring)
    """
    return skewnorm.ppf(p, eta, loc=xi, scale=omega)

skew_normal_rvs(xi, omega, eta, size=1, random_state=None)

Draw samples from Azzalini's skew-normal distribution.

Original: stats/rndSN.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_normal_rvs(
    xi: float, omega: float, eta: float, size: int = 1, random_state=None
) -> np.ndarray:
    """Draw samples from Azzalini's skew-normal distribution.

    Original: stats/rndSN.m (HSF toolbox)
    """
    return skewnorm.rvs(eta, loc=xi, scale=omega, size=size, random_state=random_state)

skew_t_cdf(x, xi, omega, eta, nu, method=0)

CDF of Azzalini's skew-t distribution (location xi, scale omega, shape eta, degrees of freedom nu).

Original: stats/cdfST.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_t_cdf(
    x: np.ndarray, xi: float, omega: float, eta: float, nu: float, method: int = 0
) -> np.ndarray:
    """CDF of Azzalini's skew-t distribution (location `xi`, scale `omega`,
    shape `eta`, degrees of freedom `nu`).

    Original: stats/cdfST.m (HSF toolbox)
    """
    xc = (np.asarray(x, dtype=float) - xi) / omega
    if method == 1:
        delta = eta / np.sqrt(1 + eta**2)
        return 2.0 * bvt_cdf(xc, 0.0, -delta, nu)

    e = float(eta >= 0)
    eta_abs = max(abs(eta), 1e-8)
    delta = (1 - eta_abs**2) / (1 + eta_abs**2)
    cdf = bvt_cdf(xc, xc, delta, nu)
    return cdf * e + (2 * student_t_cdf(xc, nu) - cdf) * (1 - e)

skew_t_moments(xi, omega, eta, nu)

Mean, std dev, skewness, and excess kurtosis of Azzalini's skew-t distribution (requires nu > 4 for the kurtosis to be finite).

Original: stats/momST.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_t_moments(
    xi: float, omega: float, eta: float, nu: float
) -> tuple[float, float, float, float]:
    """Mean, std dev, skewness, and excess kurtosis of Azzalini's skew-t
    distribution (requires nu > 4 for the kurtosis to be finite).

    Original: stats/momST.m (HSF toolbox)
    """
    delta = eta / np.sqrt(1 + eta**2)
    m0 = delta * np.sqrt(nu / np.pi) * np.exp(gammaln(0.5 * (nu - 1)) - gammaln(0.5 * nu))
    mean = xi + omega * m0
    sigma = omega * np.sqrt(nu / (nu - 2) - m0**2)

    gamma1 = (
        m0
        * (nu * (3 - delta**2) / (nu - 3) - 3 * nu / (nu - 2) + 2 * m0**2)
        * (nu / (nu - 2) - m0**2) ** (-1.5)
    )
    gamma2 = (
        3 * nu**2 / (nu - 2) / (nu - 4)
        - 4 * m0**2 * nu * (3 - delta**2) / (nu - 3)
        + 6 * m0**2 * nu / (nu - 2)
        - 3 * m0**4
    ) * (nu / (nu - 2) - m0**2) ** (-2) - 3.0

    return float(mean), float(sigma), float(gamma1), float(gamma2)

skew_t_pdf(x, xi, omega, eta, nu)

PDF of Azzalini's skew-t distribution.

Original: stats/pdfST.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_t_pdf(x: np.ndarray, xi: float, omega: float, eta: float, nu: float) -> np.ndarray:
    """PDF of Azzalini's skew-t distribution.

    Original: stats/pdfST.m (HSF toolbox)
    """
    xc = (np.asarray(x, dtype=float) - xi) / omega
    cdf = student_t_cdf(eta * xc * np.sqrt((nu + 1) / (xc**2 + nu)), nu + 1)
    pdf = student_t_pdf(xc, nu) / omega
    return 2 * pdf * cdf

skew_t_ppf(p, xi, omega, eta, nu, config=None)

Quantile function of Azzalini's skew-t distribution, via Newton iteration (no closed form / no scipy equivalent).

The default tolerance is looser than NewtonConfig's own default (1e-4 rather than 1e-10): skew_t_cdf is itself backed by bvt_cdf, whose underlying scipy.stats.multivariate_t.cdf uses a randomized quasi-Monte-Carlo integrator with irreducible call-to-call noise on the order of 1e-4 (verified empirically -- repeated calls with identical inputs vary by ~7e-5 to ~1e-4). A tighter tolerance (e.g. the 1e-8 this function used before) can never actually be satisfied, so the loop always burns its full max_iters budget without improving on the noise floor -- roughly a 4x slowdown for no accuracy gain, confirmed by direct timing. 1e-4 lets the loop exit as soon as it reaches that floor, matching the precision the underlying CDF can actually deliver.

Original: stats/cdfSTi.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_t_ppf(
    p: np.ndarray,
    xi: float,
    omega: float,
    eta: float,
    nu: float,
    config: NewtonConfig | None = None,
) -> np.ndarray:
    """Quantile function of Azzalini's skew-t distribution, via Newton
    iteration (no closed form / no scipy equivalent).

    The default tolerance is looser than `NewtonConfig`'s own default
    (`1e-4` rather than `1e-10`): `skew_t_cdf` is itself backed by
    `bvt_cdf`, whose underlying `scipy.stats.multivariate_t.cdf` uses a
    randomized quasi-Monte-Carlo integrator with irreducible call-to-call
    noise on the order of `1e-4` (verified empirically -- repeated calls
    with identical inputs vary by ~7e-5 to ~1e-4). A tighter tolerance
    (e.g. the `1e-8` this function used before) can never actually be
    satisfied, so the loop always burns its full `max_iters` budget
    without improving on the noise floor -- roughly a 4x slowdown for no
    accuracy gain, confirmed by direct timing. `1e-4` lets the loop exit
    as soon as it reaches that floor, matching the precision the
    underlying CDF can actually deliver.

    Original: stats/cdfSTi.m (HSF toolbox)
    """
    if config is None:
        config = NewtonConfig(tol=1e-4, max_iters=50)

    p = np.atleast_1d(np.asarray(p, dtype=float))
    x = student_t_ppf(p, nu)
    q = skew_t_cdf(x, 0.0, 1.0, eta, nu)
    e_max = q > 0.90
    e_min = q < 0.10
    if eta >= 0.0:
        x = 0.02 * e_min + (1 - e_min) * x
    else:
        x = -0.02 * e_max + (1 - e_max) * x

    for _ in range(config.max_iters):
        cdf = skew_t_cdf(x, 0.0, 1.0, eta, nu)
        pdf = skew_t_pdf(x, 0.0, 1.0, eta, nu)
        dx = (cdf - p) / pdf
        x = x - dx
        if np.max(np.abs(dx)) <= config.tol:
            break

    x = np.where(np.abs(p - skew_t_cdf(x, 0.0, 1.0, eta, nu)) >= 0.01, np.nan, x)
    return xi + omega * x

skew_t_rvs(xi, omega, eta, nu, size=1, random_state=None)

Draw samples from Azzalini's skew-t distribution: a skew-normal variate divided by sqrt(chi-square(nu) / nu).

Original: stats/rndST.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def skew_t_rvs(
    xi: float, omega: float, eta: float, nu: float, size: int = 1, random_state=None
) -> np.ndarray:
    """Draw samples from Azzalini's skew-t distribution: a skew-normal
    variate divided by sqrt(chi-square(nu) / nu).

    Original: stats/rndST.m (HSF toolbox)
    """
    rng = np.random.default_rng(random_state)
    n = skew_normal_rvs(0.0, omega, eta, size=size, random_state=rng)
    chi = rng.chisquare(nu, size=size)
    return xi + n / np.sqrt(chi / nu)

student_t_cdf(x, nu)

Original: stats/cdft.m

Source code in src/quanttoolbox/stats/distributions.py

def student_t_cdf(x: float | np.ndarray, nu: float) -> float | np.ndarray:
    """Original: stats/cdft.m"""
    return t.cdf(x, df=nu)

student_t_pdf(x, nu)

Original: stats/pdft.m (HSF toolbox)

Source code in src/quanttoolbox/stats/distributions.py

def student_t_pdf(x: np.ndarray, nu: float) -> np.ndarray:
    """Original: stats/pdft.m (HSF toolbox)"""
    return t.pdf(x, df=nu)

student_t_ppf(p, nu)

Original: stats/cdfti.m

Source code in src/quanttoolbox/stats/distributions.py

def student_t_ppf(p: float | np.ndarray, nu: float) -> float | np.ndarray:
    """Original: stats/cdfti.m"""
    return t.ppf(p, df=nu)

student_t_sf(x, nu)

Upper-tail (survival) Student's t. Original: stats/cdftc.m

Source code in src/quanttoolbox/stats/distributions.py

def student_t_sf(x: np.ndarray, nu: float) -> np.ndarray:
    """Upper-tail (survival) Student's t. Original: stats/cdftc.m"""
    return t.sf(x, df=nu)

Examples

Pinball-loss quantile estimation vs. sample and true quantiles — ects/quantile1.py

"""Translated from Examples/ects/quantile1.m -- illustrates that
minimizing the pinball (check) loss recovers the same quantile as the
standard sorted-sample quantile estimator, compared against the true
Normal quantile function (numeric core only; the original's comparison
plot is dropped).

The original draws y from MATLAB's unseeded `randn`; a fixed seed
(`np.random.default_rng(0)`) is substituted here. `quantile`/`sort`-based
quantile estimation uses plain `np.quantile` directly (not a ported
quanttoolbox function -- there's nothing toolbox-specific to translate
here, same as other examples that lean on bare NumPy/pandas idioms)."""

import numpy as np
from scipy.optimize import minimize_scalar

from quanttoolbox.stats.distributions import normal_ppf

rng = np.random.default_rng(0)
n_s = 1000
y = rng.standard_normal(n_s)

alpha = np.array([0.01, 0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95, 0.99])

# Empirical quantile via sorting
q1 = np.quantile(y, alpha)


def pinball_objective(q, data, a):
    u = data - q
    e = u > 0
    return a * np.sum(np.abs(u[e])) + (1 - a) * np.sum(np.abs(u[~e]))


q3 = np.array([minimize_scalar(pinball_objective, args=(y, a)).x for a in alpha])

# True Normal quantiles
q4_true = normal_ppf(alpha)

print(
    "alpha, empirical quantile (sorted), M-estimator quantile (pinball minimization), true Normal quantile:"
)
print(np.round(np.column_stack([alpha, q1, q3, q4_true]), 4))

stats.dose_response

Python alternatives

Keep — scipy has no dose-response-curve module. These are small, closed-form sigmoidal curves (log-logistic, log-normal, Weibull, and two "hormetic" variants) used in toxicology/ecotoxicology; no general-purpose equivalent found.

quanttoolbox.stats.dose_response

Dose-response curve models (toxicology/ecotoxicology): log-logistic, log-normal, and Weibull sigmoidal curves, plus two "hormetic" variants that add a low-dose stimulatory term to the log-logistic curve.

Ported from HSF toolbox stats/{drcHormetic1,drcHormetic2,drcLogLogistic, drcLogNormal,drcWeibull1,drcWeibull2}.m.

Translation notes:

• All six curves are small, closed-form functions with no scipy equivalent (scipy has no dose-response-curve module) -- ported algorithm-for-algorithm. drc_log_normal reuses this package's own quanttoolbox.stats.distributions.normal_cdf in place of the original's cdfn.
• alpha is the curve's inflection-point dose (ED50-style location parameter), beta its slope, y_min/y_max the lower/upper response asymptotes; gamma_/delta (hormetic variants only) control the extra low-dose stimulatory term. Parameter names and order match the originals exactly.

drc_hormetic1(x, alpha, beta, y_min, y_max, gamma_)

Hormetic (type I) dose-response curve: adds a linear-in-dose stimulatory term gamma_ * x to the log-logistic curve's numerator, producing a low-dose "hump" before the usual sigmoidal decline.

Original: stats/drcHormetic1.m

Source code in src/quanttoolbox/stats/dose_response.py

def drc_hormetic1(
    x: np.ndarray, alpha: float, beta: float, y_min: float, y_max: float, gamma_: float
) -> np.ndarray:
    """Hormetic (type I) dose-response curve: adds a linear-in-dose
    stimulatory term ``gamma_ * x`` to the log-logistic curve's numerator,
    producing a low-dose "hump" before the usual sigmoidal decline.

    Original: stats/drcHormetic1.m
    """
    x = np.asarray(x, dtype=float)
    y = 1.0 + np.exp(-beta * (np.log(x) - np.log(alpha)))
    return y_min + (y_max - y_min + gamma_ * x) / y

drc_hormetic2(x, alpha, beta, y_min, y_max, gamma_, delta)

Hormetic (type II) dose-response curve: like drc_hormetic1, but the stimulatory term decays as exp(-1 / x**delta) instead of growing linearly in x.

Original: stats/drcHormetic2.m

Source code in src/quanttoolbox/stats/dose_response.py

def drc_hormetic2(
    x: np.ndarray,
    alpha: float,
    beta: float,
    y_min: float,
    y_max: float,
    gamma_: float,
    delta: float,
) -> np.ndarray:
    """Hormetic (type II) dose-response curve: like `drc_hormetic1`, but
    the stimulatory term decays as ``exp(-1 / x**delta)`` instead of
    growing linearly in `x`.

    Original: stats/drcHormetic2.m
    """
    x = np.asarray(x, dtype=float)
    y = 1.0 + np.exp(-beta * (np.log(x) - np.log(alpha)))
    return y_min + (y_max - y_min + gamma_ * np.exp(-1.0 / x**delta)) / y

drc_log_logistic(x, alpha, beta, y_min, y_max)

Log-logistic dose-response curve.

Original: stats/drcLogLogistic.m

Source code in src/quanttoolbox/stats/dose_response.py

def drc_log_logistic(
    x: np.ndarray, alpha: float, beta: float, y_min: float, y_max: float
) -> np.ndarray:
    """Log-logistic dose-response curve.

    Original: stats/drcLogLogistic.m
    """
    x = np.asarray(x, dtype=float)
    y = 1.0 + np.exp(-beta * (np.log(x) - np.log(alpha)))
    return y_min + (y_max - y_min) / y

drc_log_normal(x, alpha, beta, y_min, y_max)

Log-normal dose-response curve.

Original: stats/drcLogNormal.m

Source code in src/quanttoolbox/stats/dose_response.py

def drc_log_normal(
    x: np.ndarray, alpha: float, beta: float, y_min: float, y_max: float
) -> np.ndarray:
    """Log-normal dose-response curve.

    Original: stats/drcLogNormal.m
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(normal_cdf(beta * (np.log(x) - np.log(alpha))))
    return y_min + (y_max - y_min) * y

drc_weibull1(x, alpha, beta, y_min, y_max)

Weibull (type I) dose-response curve.

Original: stats/drcWeibull1.m

Source code in src/quanttoolbox/stats/dose_response.py

def drc_weibull1(
    x: np.ndarray, alpha: float, beta: float, y_min: float, y_max: float
) -> np.ndarray:
    """Weibull (type I) dose-response curve.

    Original: stats/drcWeibull1.m
    """
    x = np.asarray(x, dtype=float)
    y = np.exp(-np.exp(beta * (np.log(x) - np.log(alpha))))
    return y_min + (y_max - y_min) * y

drc_weibull2(x, alpha, beta, y_min, y_max)

Weibull (type II) dose-response curve.

Original: stats/drcWeibull2.m

Source code in src/quanttoolbox/stats/dose_response.py

def drc_weibull2(
    x: np.ndarray, alpha: float, beta: float, y_min: float, y_max: float
) -> np.ndarray:
    """Weibull (type II) dose-response curve.

    Original: stats/drcWeibull2.m
    """
    x = np.asarray(x, dtype=float)
    y = 1.0 - np.exp(-np.exp(beta * (np.log(x) - np.log(alpha))))
    return y_min + (y_max - y_min) * y

stats.multivariate

Python alternatives

Switch: bvn_cdf/bvt_cdf are thin scipy.stats.multivariate_normal/multivariate_t-backed wrappers with the original's bivariate (x, y, rho[, nu]) call signature — kept for that call-site convenience, not because scipy lacks the capability. The 8 low-level Genz quasi-Monte-Carlo integrators in genz/ (qsimvn, qsimvt, qsilatmvnv, qsimvnauto, ...) are not hand-ported at all: they solve exactly the arbitrary-dimension MVN/MVT probability problem scipy.stats.multivariate_normal.cdf/multivariate_t.cdf already solve, via the same Genz algorithm family (same author). See the module docstring and Library alternatives for the full reasoning.

quanttoolbox.stats.multivariate

Bivariate normal/Student-t CDF and PDF wrappers, matching the original toolbox's call signature (explicit x, y, rho, nu arguments, broadcast element-wise) rather than requiring the caller to assemble a full covariance matrix.

Ported from HSF toolbox stats/{cdfbvn,pdfbvn,cdfbvt}.m.

Translation notes:

• genz/{bvn,bvnu}.m are Alan Genz's Drezner-Wesolowsky quadrature for the bivariate normal CDF; stats/cdfbvn.m is a thin wrapper around bvnu.m. scipy.stats.multivariate_normal.cdf already computes this (for any dimension, not just 2) via scipy.stats._mvn, which is itself built on Alan Genz's own Fortran mvndst routine -- the same author, the same underlying algorithm. bvn_cdf here is therefore a thin scipy-backed convenience wrapper with the original's (x, y, rho) call signature, not a hand-rolled reimplementation of bvnu.m's quadrature.
• stats/cdfbvt.m wraps MATLAB's Statistics Toolbox mvtcdf; ported the same way, via scipy.stats.multivariate_t.cdf (arbitrary dimension, available since SciPy 1.9).
• Both bvn_cdf and bvt_cdf clip rho a hair inside (-1, 1) before building the 2x2 correlation matrix: scipy's multivariate_normal/ multivariate_t raise LinAlgError on an exactly-singular (|rho| == 1) matrix rather than falling back to the (well-defined) degenerate limit, and rho can round to exactly +/-1 in float64 even when the caller's inputs are only extremely close to that boundary -- e.g. stats.distributions.skew_t_cdf hits this at eta near 0, where delta = (1 - eta^2) / (1 + eta^2) underflows to 1.0. Clipping is a robustness addition beyond a literal transliteration of the original, not a behavior change: the clipped and unclipped results agree to float64 precision anywhere scipy would have succeeded anyway.
• genz/{qsimvn,qsimvnv,qsimvt,qsimvtv,qsilatmvnv,qsilatmvtv,qsimvnauto, mvnrcnv}.m (8 of the 10 files in genz/) are randomized quasi-Monte-Carlo integrators solving exactly the same problem as scipy.stats.multivariate_normal.cdf/multivariate_t.cdf -- arbitrary- dimension MVN/MVT orthant probabilities -- via the same Genz (1992) algorithm family (Niederreiter/Cranley-Patterson randomized lattice rules); mvnrcnv.m is a constant-correlation special case of the same problem. None of these are hand-ported: doing so would duplicate already-available, better-tested scipy functionality (the same author's own algorithm, reimplemented) with no functional gain. See stats.distributions.mvn_cdf/mvn_pdf for the general n-dimensional case (already ported); bvn_cdf/bvt_cdf below just add the bivariate- specific call signature (x, y, rho rather than a point + covariance matrix) for parity with how the original toolbox's other modules call these.

bvn_cdf(x, y, rho)

Bivariate standard normal CDF P(X <= x, Y <= y), for (X, Y) with unit variances and correlation rho. Broadcasts element-wise over x/y/rho.

Original: stats/cdfbvn.m (via genz/bvnu.m)

Source code in src/quanttoolbox/stats/multivariate.py

def bvn_cdf(x: np.ndarray | float, y: np.ndarray | float, rho: np.ndarray | float) -> np.ndarray:
    """Bivariate standard normal CDF P(X <= x, Y <= y), for (X, Y) with unit
    variances and correlation `rho`. Broadcasts element-wise over
    `x`/`y`/`rho`.

    Original: stats/cdfbvn.m (via genz/bvnu.m)
    """
    x_arr, y_arr, rho_arr = np.broadcast_arrays(
        np.asarray(x, dtype=float), np.asarray(y, dtype=float), np.asarray(rho, dtype=float)
    )
    rho_arr = np.clip(rho_arr, -1.0 + 1e-8, 1.0 - 1e-8)
    out = np.empty(x_arr.shape, dtype=float)
    for idx in np.ndindex(x_arr.shape):
        cov = np.array([[1.0, rho_arr[idx]], [rho_arr[idx], 1.0]])
        out[idx] = multivariate_normal.cdf([x_arr[idx], y_arr[idx]], mean=[0.0, 0.0], cov=cov)
    return out[()] if out.shape == () else out

bvn_pdf(x1, x2, mu1, mu2, sigma1, sigma2, rho)

Bivariate normal PDF at (x1, x2), for (X1, X2) ~ N((mu1, mu2), Sigma) with std devs (sigma1, sigma2) and correlation rho. Closed-form, so (unlike bvn_cdf) evaluated directly rather than via scipy.

Original: stats/pdfbvn.m

Source code in src/quanttoolbox/stats/multivariate.py

def bvn_pdf(
    x1: np.ndarray | float,
    x2: np.ndarray | float,
    mu1: np.ndarray | float,
    mu2: np.ndarray | float,
    sigma1: np.ndarray | float,
    sigma2: np.ndarray | float,
    rho: np.ndarray | float,
) -> np.ndarray:
    """Bivariate normal PDF at (x1, x2), for (X1, X2) ~ N((mu1, mu2), Sigma)
    with std devs (sigma1, sigma2) and correlation rho. Closed-form, so
    (unlike `bvn_cdf`) evaluated directly rather than via `scipy`.

    Original: stats/pdfbvn.m
    """
    x1 = np.asarray(x1, dtype=float)
    x2 = np.asarray(x2, dtype=float)
    mu1 = np.asarray(mu1, dtype=float)
    mu2 = np.asarray(mu2, dtype=float)
    sigma1 = np.asarray(sigma1, dtype=float)
    sigma2 = np.asarray(sigma2, dtype=float)
    rho = np.asarray(rho, dtype=float)

    w = 1.0 - rho**2
    z1 = (x1 - mu1) / sigma1
    z2 = (x2 - mu2) / sigma2
    q = z1**2 - 2.0 * rho * z1 * z2 + z2**2
    return np.exp(-0.5 * q / w) / (2.0 * np.pi * sigma1 * sigma2 * np.sqrt(w))

bvt_cdf(x, y, rho, nu)

Bivariate Student-t CDF P(X <= x, Y <= y), for (X, Y) with unit scale, correlation rho, and nu degrees of freedom. Broadcasts element-wise over x/y/rho/nu.

Original: stats/cdfbvt.m (via MATLAB's mvtcdf)

Source code in src/quanttoolbox/stats/multivariate.py

def bvt_cdf(
    x: np.ndarray | float, y: np.ndarray | float, rho: np.ndarray | float, nu: np.ndarray | float
) -> np.ndarray:
    """Bivariate Student-t CDF P(X <= x, Y <= y), for (X, Y) with unit scale,
    correlation `rho`, and `nu` degrees of freedom. Broadcasts element-wise
    over `x`/`y`/`rho`/`nu`.

    Original: stats/cdfbvt.m (via MATLAB's mvtcdf)
    """
    x_arr, y_arr, rho_arr, nu_arr = np.broadcast_arrays(
        np.asarray(x, dtype=float),
        np.asarray(y, dtype=float),
        np.asarray(rho, dtype=float),
        np.asarray(nu, dtype=float),
    )
    rho_arr = np.clip(rho_arr, -1.0 + 1e-8, 1.0 - 1e-8)
    out = np.empty(x_arr.shape, dtype=float)
    for idx in np.ndindex(x_arr.shape):
        shape = np.array([[1.0, rho_arr[idx]], [rho_arr[idx], 1.0]])
        out[idx] = multivariate_t.cdf(
            [x_arr[idx], y_arr[idx]], loc=[0.0, 0.0], shape=shape, df=float(nu_arr[idx])
        )
    return out[()] if out.shape == () else out

stats.moments

Python alternatives

rolling_correlation/rolling_volatility: switch to pandas.DataFrame.rolling().corr()/.std() for the standard case — Cython-backed, meaningfully faster. Keep ours only for the method=2 "returns computed within each window" variant pandas doesn't offer. active_share, herfindahl_index, asynchronous_cov, weekly_cov: keep — no general-purpose equivalent exists.

quanttoolbox.stats.moments

Sample moments, dispersion, covariance/correlation, and portfolio-overlap measures.

Ported from QuantToolBox/stats/{skewness_coefficient,kurtosis_coefficient, herfindahl_index,mean_absolute_difference,cov2cor,cor2cov,corrx, pearson_correlation,active_share,active_share_upper_bound,asynchronous_cov, weekly_cov,rolling_correlation,rolling_volatility}.m

Translation notes:

• Every "column vector output for each column of x" pattern (MATLAB loops over cols(x)) is replaced by native NumPy/pandas vectorized operations (axis=0 reductions) instead of an explicit Python loop.
• packr (drop rows containing NaN before computing a statistic) is replaced by pandas' NaN-aware reductions (.mean(), .std(), ... skip NaN by default) or explicit ~np.isnan(...) masks where a NumPy-only implementation is used.
• Excess-kurtosis convention: like the original, kurtosis here returns the raw fourth standardized moment (3.0 for a normal distribution), not scipy's default excess-kurtosis convention (which subtracts 3).

active_share(x, b)

Portfolio active share: 0.5 * sum(|weight - benchmark weight|).

Original: stats/active_share.m

Source code in src/quanttoolbox/stats/moments.py

def active_share(x: np.ndarray, b: np.ndarray) -> float:
    """Portfolio active share: 0.5 * sum(|weight - benchmark weight|).

    Original: stats/active_share.m
    """
    return 0.5 * float(np.sum(np.abs(np.asarray(x) - np.asarray(b))))

active_share_upper_bound(b, x_minus=0.0, x_plus=1.0)

Compute the maximum achievable active share given per-position weight bounds [x_minus, x_plus] that must still sum to 1.

Greedy algorithm: iteratively push each position to whichever bound (x_minus or x_plus) is farthest from the benchmark weight, subject to the running total staying feasible.

Original: stats/active_share_upper_bound.m

Returns:

Name	Type	Description
as_max	the maximum active share achieved.
x_max	the weight vector achieving it (None if infeasible).
retcode	1 if converged exactly to a total weight of 1, else 0/-1.

Source code in src/quanttoolbox/stats/moments.py

def active_share_upper_bound(
    b: np.ndarray, x_minus: float | np.ndarray = 0.0, x_plus: float | np.ndarray = 1.0
) -> tuple[float, np.ndarray | None, int]:
    """Compute the maximum achievable active share given per-position weight
    bounds [x_minus, x_plus] that must still sum to 1.

    Greedy algorithm: iteratively push each position to whichever bound
    (x_minus or x_plus) is farthest from the benchmark weight, subject to
    the running total staying feasible.

    Original: stats/active_share_upper_bound.m

    Returns
    -------
    as_max : the maximum active share achieved.
    x_max : the weight vector achieving it (None if infeasible).
    retcode : 1 if converged exactly to a total weight of 1, else 0/-1.
    """
    b = np.asarray(b, dtype=float).flatten()
    n = b.shape[0]
    x_minus = (
        np.full(n, x_minus, dtype=float)
        if np.isscalar(x_minus)
        else np.asarray(x_minus, dtype=float)
    )
    x_plus = (
        np.full(n, x_plus, dtype=float) if np.isscalar(x_plus) else np.asarray(x_plus, dtype=float)
    )

    x_max = x_minus.copy()
    s = x_max.sum()
    if s > 1.0:
        return np.nan, None, -1

    y_minus = np.abs(b - x_minus)
    y_plus = np.abs(x_plus - b)
    y = np.concatenate([y_minus, y_plus])
    order = np.zeros(n, dtype=int)

    for it in range(n):
        if s == 1:
            break
        idx = int(np.nanargmax(y))
        if idx < n:  # "minus" candidate
            y[idx] = np.nan
            y[idx + n] = np.nan
            order[it] = idx
        else:  # "plus" candidate
            idx -= n
            ds = 1 - s + x_minus[idx]
            dx = min(ds, x_plus[idx])
            x_max[idx] = dx
            y[idx] = np.nan
            y[idx + n] = np.nan
            order[it] = idx
        s = x_max.sum()

    retcode = 1
    if not np.isclose(s, 1.0):
        retcode = 0
        for it in range(n - 1, -1, -1):
            idx = order[it]
            ds = 1 - s + x_max[idx]
            dx = min(ds, x_plus[idx])
            x_max[idx] = dx
            s = x_max.sum()
            if np.isclose(s, 1.0):
                break

    as_max = active_share(x_max, b)
    return as_max, x_max, retcode

asynchronous_cov(x, method=0)

Newey-West-style covariance correction for asynchronously-timed data (e.g. assets trading in different time zones), using a lag-1 autocovariance adjustment.

method=0: half-window correction (default). method=1: full-window.

Original: stats/asynchronous_cov.m

Source code in src/quanttoolbox/stats/moments.py

def asynchronous_cov(x: np.ndarray, method: int = 0) -> tuple[np.ndarray, np.ndarray]:
    """Newey-West-style covariance correction for asynchronously-timed data
    (e.g. assets trading in different time zones), using a lag-1
    autocovariance adjustment.

    method=0: half-window correction (default). method=1: full-window.

    Original: stats/asynchronous_cov.m
    """
    x = np.asarray(x, dtype=float)
    x = x[~np.isnan(x).any(axis=1)]
    r, c = x.shape

    xc = x - x.mean(axis=0)
    cov1 = (xc.T @ xc) / (r - 1)

    xc_lag = xc[:-1]
    xc_cur = xc[1:]
    cross = (xc_cur.T @ xc_lag) / (r - 1) + (xc_lag.T @ xc_cur) / (r - 1)

    weight = 1.0 if method == 1 else 0.5
    cov2 = cov1 + weight * cross

    sigma2 = np.sqrt(np.diag(cov2))
    cor2 = cov2 / sigma2[:, None] / sigma2[None, :]
    np.fill_diagonal(cor2, 1.0)

    sigma1 = np.sqrt(np.diag(cov1))
    cov2 = cor2 * sigma1[:, None] * sigma1[None, :]
    return cov1, cov2

corr_to_cov(sigma, rho)

Recombine standard deviations and a correlation matrix into a covariance matrix.

Original: stats/cor2cov.m

Source code in src/quanttoolbox/stats/moments.py

def corr_to_cov(sigma: np.ndarray, rho: np.ndarray) -> np.ndarray:
    """Recombine standard deviations and a correlation matrix into a covariance matrix.

    Original: stats/cor2cov.m
    """
    sigma = np.asarray(sigma, dtype=float).flatten()
    rho = np.asarray(rho, dtype=float)
    n = sigma.shape[0]
    if rho.shape != (n, n):
        raise ValueError("corr_to_cov: dimensions do not match")
    return rho * sigma[:, None] * sigma[None, :]

corrx(x)

Sample correlation matrix and covariance matrix of x (columns = variables, rows dropped if they contain NaN).

Original: stats/corrx.m

Source code in src/quanttoolbox/stats/moments.py

def corrx(x: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Sample correlation matrix and covariance matrix of x (columns = variables,
    rows dropped if they contain NaN).

    Original: stats/corrx.m
    """
    x = np.asarray(x, dtype=float)
    x = x[~np.isnan(x).any(axis=1)]
    n = x.shape[0]
    xc = x - x.mean(axis=0)
    sigma_matrix = (xc.T @ xc) / (n - 1)
    sigma = np.sqrt(np.diag(sigma_matrix))
    rho = sigma_matrix / sigma[:, None] / sigma[None, :]
    return rho, sigma_matrix

cov_to_corr(cov_matrix)

Decompose a covariance matrix into (std devs, correlation matrix).

Original: stats/cov2cor.m

Source code in src/quanttoolbox/stats/moments.py

def cov_to_corr(cov_matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Decompose a covariance matrix into (std devs, correlation matrix).

    Original: stats/cov2cor.m
    """
    cov_matrix = np.asarray(cov_matrix, dtype=float)
    sigma = np.sqrt(np.diag(cov_matrix))
    rho = cov_matrix / sigma[:, None] / sigma[None, :]
    return sigma, rho

herfindahl_index(x, b=None)

Herfindahl index H = sum(x^2), its inverse N = 1/H (effective number of positions), and optionally the same for a benchmark b.

Original: stats/herfindahl_index.m

Source code in src/quanttoolbox/stats/moments.py

def herfindahl_index(
    x: np.ndarray, b: np.ndarray | None = None
) -> tuple[np.ndarray, np.ndarray, np.ndarray | None]:
    """Herfindahl index H = sum(x^2), its inverse N = 1/H (effective number of
    positions), and optionally the same for a benchmark b.

    Original: stats/herfindahl_index.m
    """
    x = np.asarray(x, dtype=float)
    h = np.sum(x**2, axis=0)
    n = 1.0 / h
    if b is None:
        return h, n, None
    b = np.asarray(b, dtype=float)
    h_b = np.sum(b**2, axis=0)
    n_b = h_b / h
    return h, n, n_b

kurtosis(x)

Sample kurtosis (raw, not excess) of each column of x. Original: stats/kurtosis_coefficient.m

Source code in src/quanttoolbox/stats/moments.py

def kurtosis(x: np.ndarray) -> np.ndarray:
    """Sample kurtosis (raw, not excess) of each column of x.
    Original: stats/kurtosis_coefficient.m
    """
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x[:, None]
    out = np.full(x.shape[1], np.nan)
    for i in range(x.shape[1]):
        y = x[:, i]
        y = y[~np.isnan(y)]
        yc = y - y.mean()
        m2 = np.mean(yc**2)
        m4 = np.mean(yc**4)
        out[i] = m4 / m2**2
    return out

mean_absolute_difference(x)

Gini-style mean absolute pairwise difference of each column of x.

Original: stats/mean_absolute_difference.m

Source code in src/quanttoolbox/stats/moments.py

def mean_absolute_difference(x: np.ndarray) -> np.ndarray:
    """Gini-style mean absolute pairwise difference of each column of x.

    Original: stats/mean_absolute_difference.m
    """
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x[:, None]
    out = np.full(x.shape[1], np.nan)
    for i in range(x.shape[1]):
        y = x[:, i]
        y = y[~np.isnan(y)]
        n = y.shape[0]
        dx = y[:, None] - y[None, :]
        out[i] = np.sum(np.abs(dx)) / n**2
    return out

pearson_correlation(x, y)

Column-wise (or cross) Pearson correlation between x and y.

If x and y have the same number of columns, returns the correlation of each matching pair of columns. Otherwise returns the full nX x nY cross-correlation matrix.

Original: stats/pearson_correlation.m

Source code in src/quanttoolbox/stats/moments.py

def pearson_correlation(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    """Column-wise (or cross) Pearson correlation between x and y.

    If x and y have the same number of columns, returns the correlation of
    each matching pair of columns. Otherwise returns the full nX x nY
    cross-correlation matrix.

    Original: stats/pearson_correlation.m
    """
    x = np.atleast_2d(np.asarray(x, dtype=float))
    y = np.atleast_2d(np.asarray(y, dtype=float))
    if x.shape[0] == 1 and x.shape[1] > 1:
        x = x.T
    if y.shape[0] == 1 and y.shape[1] > 1:
        y = y.T

    n_x, n_y = x.shape[1], y.shape[1]
    d_x = x - x.mean(axis=0)
    sigma_x = np.sqrt(np.mean(d_x**2, axis=0))

    if n_x == n_y:
        d_y = y - y.mean(axis=0)
        sigma_y = np.sqrt(np.mean(d_y**2, axis=0))
        return np.mean(d_x * d_y, axis=0) / (sigma_x * sigma_y)

    rho = np.zeros((n_x, n_y))
    for i in range(n_y):
        yi = y[:, i]
        d_y = yi - yi.mean()
        sigma_y = np.sqrt(np.mean(d_y**2))
        rho[:, i] = np.mean(d_x * d_y[:, None], axis=0) / (sigma_x * sigma_y)
    return rho

rolling_correlation(x, y, n_lags, method=0)

Rolling (trailing n_lags-window) correlation, and each series' rolling volatility, column-by-column.

method=1 treats both x and y as prices and converts to returns up front. method=2 converts to returns within each window (useful for illiquid series where a fixed global return series would drop too many rows).

Original: stats/rolling_correlation.m

Source code in src/quanttoolbox/stats/moments.py

def rolling_correlation(
    x: np.ndarray, y: np.ndarray, n_lags: int, method: int = 0
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Rolling (trailing n_lags-window) correlation, and each series' rolling
    volatility, column-by-column.

    method=1 treats both x and y as prices and converts to returns up front.
    method=2 converts to returns within each window (useful for illiquid
    series where a fixed global return series would drop too many rows).

    Original: stats/rolling_correlation.m
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if x.ndim == 1:
        x = x[:, None]
    if y.ndim == 1:
        y = y[:, None]
    if x.shape != y.shape:
        raise ValueError("rolling_correlation: x and y do not match")

    r, c = x.shape
    rho = np.full((r, c), np.nan)
    sigma_x = np.full((r, c), np.nan)
    sigma_y = np.full((r, c), np.nan)

    if method == 1:
        x = x[1:] / x[:-1] - 1.0
        y = y[1:] / y[:-1] - 1.0
        r -= 1

    for i in range(n_lags, r):
        for j in range(c):
            wx = x[i - n_lags : i + 1, j]
            wy = y[i - n_lags : i + 1, j]
            valid = ~np.isnan(wx) & ~np.isnan(wy)
            wx, wy = wx[valid], wy[valid]

            if method == 2:
                wx = wx[1:] / wx[:-1] - 1.0
                wy = wy[1:] / wy[:-1] - 1.0

            min_obs = 2 if method == 2 else 0.5 * n_lags
            if wx.shape[0] >= min_obs:
                sigma_x[i, j] = wx.std(ddof=1)
                sigma_y[i, j] = wy.std(ddof=1)
                rho[i, j] = np.corrcoef(wx, wy)[0, 1]

    return rho, sigma_x, sigma_y

rolling_volatility(x, n_lags, method=0)

Rolling (trailing n_lags-window) volatility of each column of x.

method=1 treats x as prices and converts to simple returns first.

Original: stats/rolling_volatility.m

Source code in src/quanttoolbox/stats/moments.py

def rolling_volatility(x: np.ndarray, n_lags: int, method: int = 0) -> np.ndarray:
    """Rolling (trailing n_lags-window) volatility of each column of x.

    method=1 treats x as prices and converts to simple returns first.

    Original: stats/rolling_volatility.m
    """
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x[:, None]

    if method == 1:
        x = x[1:] / x[:-1] - 1.0

    r, c = x.shape
    sigma = np.full((r, c), np.nan)
    for i in range(n_lags, r):
        for j in range(c):
            window = x[i - n_lags : i + 1, j]
            window = window[~np.isnan(window)]
            if window.shape[0] >= 0.5 * n_lags:
                sigma[i, j] = window.std(ddof=1)
    return sigma

skewness(x)

Sample skewness of each column of x. Original: stats/skewness_coefficient.m

Source code in src/quanttoolbox/stats/moments.py

def skewness(x: np.ndarray) -> np.ndarray:
    """Sample skewness of each column of x. Original: stats/skewness_coefficient.m"""
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x[:, None]
    out = np.full(x.shape[1], np.nan)
    for i in range(x.shape[1]):
        y = x[:, i]
        y = y[~np.isnan(y)]
        yc = y - y.mean()
        m2 = np.mean(yc**2)
        m3 = np.mean(yc**3)
        out[i] = m3 / m2**1.5
    return out

weekly_cov(x_weekly, x_daily=None)

Covariance matrix using weekly correlations but daily-data volatilities (useful when daily vol estimates are more reliable than weekly ones but weekly correlations are less affected by asynchronous trading).

Original: stats/weekly_cov.m

Source code in src/quanttoolbox/stats/moments.py

def weekly_cov(
    x_weekly: np.ndarray, x_daily: np.ndarray | None = None
) -> tuple[np.ndarray, np.ndarray | None, np.ndarray | None]:
    """Covariance matrix using weekly correlations but daily-data volatilities
    (useful when daily vol estimates are more reliable than weekly ones but
    weekly correlations are less affected by asynchronous trading).

    Original: stats/weekly_cov.m
    """
    x = np.asarray(x_weekly, dtype=float)
    x = x[~np.isnan(x).any(axis=1)]
    r, c = x.shape
    xc = x - x.mean(axis=0)
    vcv = (xc.T @ xc) / (r - 1)

    if x_daily is None:
        return vcv, None, None

    sigma_w = np.sqrt(np.diag(vcv))
    rho = vcv / sigma_w[:, None] / sigma_w[None, :]

    xd = np.asarray(x_daily, dtype=float)
    xd = xd[~np.isnan(xd).any(axis=1)]
    sigma_d = xd.std(axis=0, ddof=1)

    vcv_out = rho * sigma_d[:, None] * sigma_d[None, :]
    return vcv_out, rho, sigma_d

Examples

Black-Litterman view sensitivity across five scenarios — rpb/test_bl2.py

"""Translated from Examples/rpb/test_bl2.m -- Black-Litterman sensitivity
analysis across 6 scenarios (base case + 5 view/uncertainty/tau variants),
each solved as a fixed-risk-aversion MVO portfolio."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.black_litterman import black_litterman_moments, implied_risk_premia
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

x0 = np.array([0.40, 0.30, 0.20, 0.10])
r = 0.03
irp = implied_risk_premia(x0, cov_matrix, sharpe_ratio=0.25)
mu_tilde = r + irp.pi
gamma0 = irp.gamma

scenarios = [
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.07, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.20**2, 0.20**2]),
        tau=1,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=0.10,
    ),
    dict(
        P=np.array([[1, 0, 0, 0], [0, 1, -1, 0]], dtype=float),
        Q=np.array([0.04, -0.01]),
        Omega=np.diag([0.10**2, 0.05**2]),
        tau=0.01,
    ),
]

results = [
    dict(weights=x0, mu=x0 @ mu_tilde, sigma=np.sqrt(x0 @ cov_matrix @ x0), alpha=0.0, te=0.0)
]
for s in scenarios:
    bl = black_litterman_moments(mu_tilde, s["tau"] * cov_matrix, s["P"], s["Q"], s["Omega"])
    mvo = mvo_portfolio(bl.mu_bar, cov_matrix, gamma=gamma0, lb=0.0, ub=1.0)
    alpha = (mvo.weights - x0) @ bl.mu_bar
    te = np.sqrt((mvo.weights - x0) @ cov_matrix @ (mvo.weights - x0))
    results.append(
        dict(weights=mvo.weights, mu=mvo.expected_return, sigma=mvo.volatility, alpha=alpha, te=te)
    )

for i, r_ in enumerate(results):
    print(
        f"scenario {i}: weights={np.round(r_['weights'],4)} mu={round(r_['mu'],5)} te={round(r_['te'],5)}"
    )

Equal risk contribution and box-constrained risk budgeting — rpb/test_box1.py

"""Translated from Examples/rpb/test_box1.m -- ERC and box-constrained
("C-ERC") risk budgeting portfolios at progressively wider bounds around
a starting position."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import erc_portfolio, solve_box_constrained
from quanttoolbox.stats.moments import corr_to_cov

x0 = np.array([0.29, 0.25, 0.23, 0.18, 0.05])
sigma = np.array([0.20, 0.20, 0.25, 0.15, 0.25])
rho = xpnd(
    np.array(
        [1.00, 0.40, 1.00, 0.70, 0.75, 1.00, 0.60, 0.55, 0.90, 1.00, 0.70, 0.60, 0.70, 0.65, 1.00]
    ),
    method=1,
)
cov_matrix = corr_to_cov(sigma, rho)

r1 = erc_portfolio(cov_matrix)
print("ERC weights:", np.round(r1.weights, 4))

for delta in [0.02, 0.07, 0.20]:
    x_minus, x_plus = x0 - delta, x0 + delta
    r = solve_box_constrained(cov_matrix, x_minus=x_minus, x_plus=x_plus, x0=x0)
    print(f"box (delta={delta}) weights:", np.round(r.weights, 4), "converged:", r.converged)

Mean-variance frontier: risk-aversion, return, and volatility targets — rpb/test_mvo2.py

"""Translated from Examples/rpb/test_mvo2.m -- Roncalli [2013], "Introduction
to Risk Parity and Budgeting", Example 1 (pages 7-8): the same 4-asset
mean-variance problem evaluated three ways -- the gamma-problem (pick a
risk-aversion, solve directly), the mu-problem (pick a target expected
return, bisect on gamma to hit it), and the sigma-problem (pick a target
volatility, bisect on gamma to hit it). All three route through
`compute_mvo_portfolio.m`'s three branches; here that's
`mvo_frontier`/`mvo_target_portfolio`.

The original passes `lb=0, ub=0` (MATLAB's "use the default -100/100 wide
bounds" sentinel); passed through explicitly here as `lb=-100, ub=100`."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import mvo_frontier, mvo_target_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

print("1. gamma-problem (page 7)")
gamma_values = np.array([0.00, 0.20, 0.50, 1.00, 2.00, 5.00])
results = mvo_frontier(mu, cov_matrix, gamma_values, lb=-100.0, ub=100.0)
for g, r in zip(gamma_values, results, strict=False):
    print(
        f"  gamma={g:5.2f}  mu={100 * r.expected_return:6.2f}  sigma={100 * r.volatility:6.2f}  "
        f"w={np.round(100 * r.weights, 2)}"
    )

print("\n2. mu-problem (page 8)")
mu_targets = np.array([5.00, 6.00, 7.00, 8.00, 9.00]) / 100
mu_results = mvo_target_portfolio(mu, cov_matrix, mu_targets, problem="mu", lb=-100.0, ub=100.0)
for target, r in zip(mu_targets, mu_results, strict=False):
    print(
        f"  target_mu={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
        f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
    )

print("\n3. sigma-problem (page 8)")
sigma_targets = np.array([15.00, 20.00, 25.00, 30.00, 35.00]) / 100
sigma_results = mvo_target_portfolio(
    mu, cov_matrix, sigma_targets, problem="sigma", lb=-100.0, ub=100.0
)
for target, r in zip(sigma_targets, sigma_results, strict=False):
    print(
        f"  target_sigma={100 * target:5.2f}  gamma={r.gamma:6.3f}  mu={100 * r.expected_return:6.2f}  "
        f"sigma={100 * r.volatility:6.2f}  w={np.round(100 * r.weights, 2)}"
    )

Mean-variance optimization plus ridge/lasso-penalized portfolios — rpb/test_lasso1.py

"""Translated from Examples/rpb/test_lasso1.m -- Roncalli [2013],
"Introduction to Risk Parity and Budgeting", Example 1 (page 53):
compares an unconstrained gamma-problem MVO portfolio against the same
problem with a budget (sum-to-1) constraint, and ridge/lasso-penalized
variants (toward 0 and toward equal-weight) of the unconstrained
problem.

`compute_mvo_portfolio`/`quadprog_ridge`/`quadprog_lasso` map to
`mvo_portfolio`/`solve_qp(..., ridge_penalty=...)`/`solve_qp(...,
lasso_penalty=...)` as established in test_lasso3.py/test_lasso5.py; the
original's `0, 0` sentinel arguments for "no equality constraint" /
"no bounds" map to simply omitting those keyword arguments."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)

# Case gamma-problem
gamma_x = 0.5

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x).weights
x2 = mvo_portfolio(
    mu, cov_matrix, gamma=gamma_x, a_eq=np.ones((1, n)), b_eq=np.array([1.0])
).weights

lambda_ridge = 0.03
s_ridge = lambda_ridge * np.eye(n)
x3 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, np.zeros(n)))
x4 = solve_qp(cov_matrix, gamma_x * mu, ridge_penalty=(s_ridge, x0))

lambda_lasso = 0.03 / 2
s_lasso = lambda_lasso * np.ones(n)
x5 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, np.zeros(n)))
x6 = solve_qp(cov_matrix, gamma_x * mu, lasso_penalty=(s_lasso, x0))

results = 100 * np.column_stack([x1, x2, x3, x4, x5, x6])
print("            x1      x2      x3      x4      x5      x6")
print(np.round(results, 2))

Minimum-variance portfolio under general linear constraints — rpb/test_minvar2.py

"""Translated from Examples/rpb/test_minvar2.m -- minimum-variance
portfolio with general linear equality/inequality constraints and box
bounds."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.mean_variance import minvar_portfolio
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)

a_eq = np.array([[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 0.0, 0.0]])
b_eq = np.array([1.0, 0.0])
c_ineq = np.array([[0.0, 0.0, 0.0, -1.0]])
d_ineq = np.array([-0.90])

r = minvar_portfolio(
    cov_matrix, a_eq=a_eq, b_eq=b_eq, c_ineq=c_ineq, d_ineq=d_ineq, lb=-1.50, ub=2.00
)
print("weights:", np.round(r.weights, 3))
print("volatility:", round(r.volatility, 5))

Mixed ridge+lasso penalties toward two different target vectors — rpb/test_lasso5.py

"""Translated from Examples/rpb/test_lasso5.m -- ridge/lasso/mixed
portfolios with two different ridge/lasso target vectors (equal-weight
and a custom 20/20/30/30 target)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
a_eq, b_eq = np.ones((1, n)), np.array([1.0])
gamma_x = 0.5

y1 = np.full(n, 1 / n)
y2 = np.array([0.20, 0.20, 0.30, 0.30])
S_ridge = np.diag(np.diag(cov_matrix))
lambda_lasso = 0.005

x1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0).weights
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, y1)
)
x4 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, y1)
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y1),
)
# mixed with DIFFERENT targets for ridge (toward y1) vs lasso (toward y2)
x8 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, y1),
    lasso_penalty=(lambda_lasso, y2),
)

for name, x in [
    ("MVO", x1),
    ("Ridge->y1", x2),
    ("Lasso->y1", x4),
    ("Mixed(ridge->y1,lasso->y1)", x6),
    ("Mixed(ridge->y1,lasso->y2)", x8),
]:
    print(f"{name}: {np.round(x, 4)}")

Ridge, lasso, and mixed-norm penalized portfolios — rpb/test_lasso3.py

"""Translated from Examples/rpb/test_lasso3.m (byte-identical to
test_lasso4.m) -- MVO/ridge/lasso/mixed-penalty portfolios via the
consolidated solve_qp (replacing the original's separate quadprog_ridge/
quadprog_lasso/quadprog_mixed calls -- see optim/quadprog.py docstring)."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.optim.quadprog import solve_qp
from quanttoolbox.portfolio.mean_variance import mvo_portfolio
from quanttoolbox.stats.moments import corr_to_cov

mu = np.array([0.05, 0.06, 0.08, 0.06])
sigma = np.array([0.15, 0.20, 0.25, 0.30])
rho = xpnd(np.array([1.00, 0.10, 1.00, 0.40, 0.70, 1.00, 0.50, 0.40, 0.80, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
n = 4
x0 = np.full(n, 1 / n)
gamma_x = 0.5
a_eq, b_eq = np.ones((1, n)), np.array([1.0])

# gamma-problem (plain MVO)
r1 = mvo_portfolio(mu, cov_matrix, gamma=gamma_x, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0)
x1 = r1.weights

# ridge-problem: penalize toward zero, scaled by each asset's own variance
S_ridge = np.diag(np.diag(cov_matrix))
x2 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, ridge_penalty=(S_ridge, x0)
)

# lasso-problem: L1 penalty toward equal weight
lambda_lasso = 0.005
x3 = solve_qp(
    cov_matrix, gamma_x * mu, a_eq=a_eq, b_eq=b_eq, lb=0.0, ub=1.0, lasso_penalty=(lambda_lasso, x0)
)

# mixed-problem: both ridge and lasso, toward various targets
x4 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x5 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, np.zeros(n)),
)
x6 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, np.zeros(n)),
    lasso_penalty=(lambda_lasso, x0),
)
x7 = solve_qp(
    cov_matrix,
    gamma_x * mu,
    a_eq=a_eq,
    b_eq=b_eq,
    lb=0.0,
    ub=1.0,
    ridge_penalty=(S_ridge, x0),
    lasso_penalty=(lambda_lasso, x0),
)

for name, x in [
    ("MVO", x1),
    ("Ridge", x2),
    ("Lasso", x3),
    ("Mixed(0,0)", x4),
    ("Mixed(EW,0)", x5),
    ("Mixed(0,EW)", x6),
    ("Mixed(EW,EW)", x7),
]:
    print(f"{name}: {np.round(x, 4)}")

Risk budgeting toward unequal target budgets — rpb/test_erc3.py

"""Translated from Examples/rpb/test_erc3.m -- Example 17 (page 123),
Roncalli (2013). RB portfolio with unequal target budgets."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import risk_contribution, solve_unconstrained
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.15, 0.20, 0.30, 0.10])
rho = xpnd(np.array([1.00, 0.50, 1.00, 0.00, 0.20, 1.00, -0.10, 0.40, 0.70, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
x = np.full(4, 0.25)

rc = risk_contribution(x, cov_matrix)
print("equal-weight risk contribution:", rc.risk, np.round(100 * rc.pct_risk_contribution, 2))

b = np.array([0.20, 0.20, 0.30, 0.30])
r = solve_unconstrained(cov_matrix, b=b, method="ccd")
print("RB weights (target budgets 20/20/30/30):", np.round(r.weights, 4))
print("converged:", r.converged, "n_iters:", r.n_iters)

Risk contribution decomposition and equal-budget risk budgeting — rpb/test_erc2.py

"""Translated from Examples/rpb/test_erc2.m -- Example 7 (page 80),
Roncalli (2013). Risk contribution decomposition + risk budgeting at
various target budgets."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.portfolio.risk_budgeting import risk_contribution, solve_unconstrained
from quanttoolbox.stats.moments import corr_to_cov

sigma = np.array([0.30, 0.20, 0.15])
rho = xpnd(np.array([1.00, 0.80, 1.00, 0.50, 0.30, 1.00]), method=1)
cov_matrix = corr_to_cov(sigma, rho)
x = np.array([0.50, 0.20, 0.30])

rc = risk_contribution(x, cov_matrix)
print("fixed-x risk contribution:", rc.risk, np.round(100 * rc.pct_risk_contribution, 2))

r_equal = solve_unconstrained(cov_matrix, b=np.full(3, 1 / 3), method="ccd")
print("equal-budget RB weights:", np.round(r_equal.weights, 4))

r_custom = solve_unconstrained(cov_matrix, b=x, method="ccd")
print("target-x-as-budget RB weights:", np.round(r_custom.weights, 4))

stats.regression.ols

Python alternatives

Hybrid: statsmodels.OLS/WLS is more complete for plain unrestricted regression. Keep this module for the restriction=(RR, r) linear-restriction parameterization, which statsmodels doesn't support as directly.

quanttoolbox.stats.regression.ols

Ordinary least squares, centering/standardization, conditional-normal regression, and principal component analysis.

Ported from QuantToolBox/stats/{regOLS,regCenter,regStandardize,regCND, regPCA}.m

Translation notes:

• regOLS uses numpy.linalg.lstsq (QR-based) internally for numerical stability, rather than MATLAB's inv(x'*x) * x'*y normal equations, but still returns the explicit inv(x'x) covariance factor since downstream code needs it for standard errors.
• regPCA decomposes a correlation matrix (following the original), using numpy.linalg.eigh (the matrix is symmetric by construction) instead of MATLAB's eig; eigenvalues/vectors are then sorted descending to match the original's convention.
• The MATLAB global Print_Results diagnostic-printing branch in regPCA is dropped -- callers should print/inspect the returned PCAResult fields themselves.

center(x)

Center each column of x by its (NaN-dropped) mean.

Original: stats/regCenter.m

Source code in src/quanttoolbox/stats/regression/ols.py

def center(x: np.ndarray) -> np.ndarray:
    """Center each column of x by its (NaN-dropped) mean.

    Original: stats/regCenter.m
    """
    x = np.asarray(x, dtype=float)
    valid = x[~np.isnan(x).any(axis=1)] if x.ndim > 1 else x[~np.isnan(x)]
    return x - valid.mean(axis=0)

conditional_normal_regression(mu_y=None, mu_x=None, sigma_yy=None, sigma_yx=None, sigma_xx=None, *, mu=None, sigma=None)

Linear regression coefficients implied by a joint normal distribution's moments (rather than estimated from a data sample).

Two calling conventions, matching the two MATLAB entry points:

• conditional_normal_regression(mu_y, mu_x, sigma_yy, sigma_yx, sigma_xx): regress a single y on x given the joint (y, x) moments.
• conditional_normal_regression(mu=mu, sigma=sigma): leave-one-out regression of each variable in mu/sigma on all the others.

Original: stats/regCND.m

Source code in src/quanttoolbox/stats/regression/ols.py

def conditional_normal_regression(
    mu_y: np.ndarray | None = None,
    mu_x: np.ndarray | None = None,
    sigma_yy: float | None = None,
    sigma_yx: np.ndarray | None = None,
    sigma_xx: np.ndarray | None = None,
    *,
    mu: np.ndarray | None = None,
    sigma: np.ndarray | None = None,
) -> ConditionalNormalResult:
    """Linear regression coefficients implied by a joint normal distribution's
    moments (rather than estimated from a data sample).

    Two calling conventions, matching the two MATLAB entry points:

    - ``conditional_normal_regression(mu_y, mu_x, sigma_yy, sigma_yx, sigma_xx)``:
      regress a single y on x given the joint (y, x) moments.
    - ``conditional_normal_regression(mu=mu, sigma=sigma)``: leave-one-out
      regression of each variable in ``mu``/``sigma`` on all the others.

    Original: stats/regCND.m
    """
    if mu is not None and sigma is not None:
        return _conditional_normal_leave_one_out(mu, sigma)

    if mu_y is None or mu_x is None or sigma_yy is None or sigma_yx is None or sigma_xx is None:
        raise ValueError(
            "conditional_normal_regression: provide either (mu, sigma) or "
            "all of (mu_y, mu_x, sigma_yy, sigma_yx, sigma_xx)"
        )

    sigma_xy = sigma_yx.T
    sigma_xx_inv = np.linalg.inv(sigma_xx)
    beta0 = mu_y - sigma_yx @ sigma_xx_inv @ mu_x
    beta = sigma_xx_inv @ sigma_xy
    sigma2 = sigma_yy - sigma_yx @ sigma_xx_inv @ sigma_xy
    r_squared = 1 - sigma2 / sigma_yy
    sigma_out = np.sqrt(sigma2)
    return ConditionalNormalResult(beta0=beta0, beta=beta, sigma=sigma_out, r_squared=r_squared)

ols(y, x)

Ordinary least squares regression, dropping rows with missing y or x.

Original: stats/regOLS.m

Source code in src/quanttoolbox/stats/regression/ols.py

def ols(y: np.ndarray, x: np.ndarray) -> OLSResult:
    """Ordinary least squares regression, dropping rows with missing y or x.

    Original: stats/regOLS.m
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    nobs = y.shape[0]
    nvar = x.shape[1]

    valid = ~np.isnan(y) & ~np.isnan(x).any(axis=1)
    y_valid, x_valid = y[valid], x[valid]

    xx = x_valid.T @ x_valid
    inv_xx = np.linalg.inv(xx)
    beta = inv_xx @ (x_valid.T @ y_valid)

    u = y_valid - x_valid @ beta
    sigma = u.std(ddof=1)
    vcv = sigma**2 * inv_xx
    stderr = np.sqrt(np.diag(vcv))

    residuals = np.full(nobs, np.nan)
    residuals[valid] = u

    return OLSResult(
        beta=beta, stderr=stderr, vcv=vcv, residuals=residuals, nobs=nobs, nvar=nvar, sigma=sigma
    )

pca(x, num_factors=None, normalize=False)

Principal component analysis on a data matrix (standardized first) or directly on a pre-computed correlation matrix.

Original: stats/regPCA.m

Source code in src/quanttoolbox/stats/regression/ols.py

def pca(x: np.ndarray, num_factors: int | None = None, normalize: bool = False) -> PCAResult:
    """Principal component analysis on a data matrix (standardized first) or
    directly on a pre-computed correlation matrix.

    Original: stats/regPCA.m
    """
    x = np.asarray(x, dtype=float)
    if x.shape[0] == x.shape[1]:
        corr = x
    else:
        x_std = (x - x.mean(axis=0)) / x.std(axis=0, ddof=1)
        corr = (x_std.T @ x_std) / x_std.shape[0]

    n = corr.shape[0]
    if num_factors is None or num_factors == 0:
        num_factors = n

    eigenvalues, eigenvectors = np.linalg.eigh(corr)
    order = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[order][:num_factors]
    eigenvectors = eigenvectors[:, order][:, :num_factors]

    quality = eigenvalues / np.sum(np.diag(corr))
    cum_quality = np.cumsum(quality)

    saturation = eigenvectors * np.sqrt(eigenvalues)
    variable_quality = saturation**2
    variable_contribution = variable_quality / eigenvalues

    if normalize:
        variable_quality = variable_quality / variable_quality.sum(axis=1, keepdims=True)
        variable_contribution = variable_contribution / variable_contribution.sum(
            axis=0, keepdims=True
        )

    return PCAResult(
        loadings=eigenvectors,
        eigenvalues=eigenvalues,
        quality=quality,
        cum_quality=cum_quality,
        saturation=saturation,
        variable_quality=variable_quality,
        variable_contribution=variable_contribution,
    )

standardize(x)

Center and scale each column of x to unit variance (using the NaN-dropped mean/std).

Original: stats/regStandardize.m

Source code in src/quanttoolbox/stats/regression/ols.py

def standardize(x: np.ndarray) -> np.ndarray:
    """Center and scale each column of x to unit variance (using the
    NaN-dropped mean/std).

    Original: stats/regStandardize.m
    """
    x = np.asarray(x, dtype=float)
    valid = x[~np.isnan(x).any(axis=1)] if x.ndim > 1 else x[~np.isnan(x)]
    mean = valid.mean(axis=0)
    std = valid.std(axis=0, ddof=1)
    return (x - mean) / std

Examples

Elastic-net regression path — stats/elasticnet1.py

"""Translated from Examples/stats/elasticnet1.m and elasticnet2.m --
elastic net regression path at two different alpha (L1/L2 mix) values,
on the same 15-observation dataset used in ridge1.m/lasso1.m."""

import numpy as np

from quanttoolbox.stats.regression.lasso import elastic_net_ccd
from quanttoolbox.stats.regression.ols import standardize

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])

for alpha, lam in [(0.50, 1.0), (0.25, 1.0)]:
    beta, path = elastic_net_ccd(y, x, lambda_=lam, alpha=alpha, n_iters=200)
    print(f"alpha={alpha}, lambda={lam}: beta={np.round(beta,4)}")

Lasso path with R-squared, degrees of freedom, and complexity — stats/lasso1.py

"""Translated from Examples/stats/lasso1.m -- penalized-form lasso path at
5 lambda values (numeric core only; the original's 501-point tau-sweep and
its multi-line coefficient-path plot are dropped -- the same lambda
values and dataset are already used by lasso2.py, which this extends with
R^2/degrees-of-freedom/complexity reporting to match lasso1.m's own
output table)."""

import numpy as np

from quanttoolbox.stats.regression.lasso import lasso_ccd
from quanttoolbox.stats.regression.ols import standardize

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])

beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)
tss = np.mean(y**2)

print("lambda  |beta|         tau=sum|beta| rss     R2      df  complexity")
for lam in [0.0, 0.9, 2.5, 5.5, 7.5]:
    beta, _ = lasso_ccd(y, x, lambda_=lam, n_iters=200)
    u = y - x @ beta
    rss = np.mean(u**2)
    r2 = 1 - rss / tss
    tau = np.sum(np.abs(beta))
    df = int(np.sum(np.abs(beta) >= 1e-6))
    complexity = 1.0 / df if df > 0 else np.inf
    print(
        f"{lam:6.2f}  {np.round(beta, 4)}  tau={tau:.4f}  rss={rss:.4f}  "
        f"R2={r2:.4f}  df={df}  complexity={complexity:.4f}"
    )

Norm-budget (tau-targeted) ridge regression vs. fixed-lambda ridge — stats/ridge2.py

"""Translated from Examples/stats/ridge2.m -- tau-targeted ridge regression
(`ridge_tau_targeted`) cross-checked against fixed-lambda ridge (`ridge`),
in both absolute and OLS-relative tau modes, on the same 15-observation
dataset used throughout stats/."""

import numpy as np

from quanttoolbox.stats.regression.ols import standardize
from quanttoolbox.stats.regression.ridge import ridge, ridge_tau_targeted

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])

beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)
tau_ols = np.sum(beta_ols**2)
print("tau(ols) =", round(tau_ols, 4))

# `ridge()` returns (beta, df, complexity) -- it doesn't report the
# resulting L2 budget, so tau_ridge (the quantity `ridge_tau_targeted`
# was aiming for) is recomputed here as sum(beta**2), to cross-check that
# fixed-lambda ridge at `lambda_out` reproduces the same beta/tau
# `ridge_tau_targeted` found by its grid search.

# Absolute tau targets
tau = np.array([0.7, 0.9, 1.1])
beta_ridge2, lambda_out, df_ridge2, complexity2 = ridge_tau_targeted(y, x, tau)
beta_ridge, df_ridge, complexity = ridge(y, x, lambda_out)
tau_ridge = np.sum(beta_ridge**2, axis=1)

print("\nAbsolute-tau analysis: tau target, tau achieved (fixed-lambda ridge), lambda")
print(np.round(np.column_stack([tau, tau_ridge, lambda_out]), 4))

# Relative tau targets (fraction of tau_ols)
tau = np.array([0.25, 0.50, 0.75, 1.00])
beta_ridge2, lambda_out, df_ridge2, complexity2 = ridge_tau_targeted(y, x, tau, relative=True)
beta_ridge, df_ridge, complexity = ridge(y, x, lambda_out)
tau_ridge_rel = np.sum(beta_ridge**2, axis=1) / tau_ols

print("\nRelative-tau analysis: tau target, tau achieved / tau_ols, lambda")
print(np.round(np.column_stack([tau, tau_ridge_rel, lambda_out]), 4))
print("\nbeta (fixed-lambda ridge) at each relative-tau lambda:")
print(np.round(beta_ridge, 4))
print("\nbeta (tau-targeted ridge):")
print(np.round(beta_ridge2, 4))

# Finer lambda search grid
lambda_search = np.arange(0, 15, 0.001)
beta_ridge2, lambda_out, df_ridge2, complexity2 = ridge_tau_targeted(
    y, x, tau, lambda_search=lambda_search, relative=True
)
beta_ridge, df_ridge, complexity = ridge(y, x, lambda_out)
tau_ridge_rel = np.sum(beta_ridge**2, axis=1) / tau_ols

print("\nRelative-tau analysis, finer lambda grid: tau target, tau achieved / tau_ols, lambda")
print(np.round(np.column_stack([tau, tau_ridge_rel, lambda_out]), 4))

Penalized-form lasso path via coordinate descent — stats/lasso2.py

"""Translated from Examples/stats/lasso2.m -- penalized-form lasso
(regLasso2) at several lambda values, on the same 15-obs dataset."""

import numpy as np

from quanttoolbox.stats.regression.lasso import lasso_ccd
from quanttoolbox.stats.regression.ols import standardize

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])
beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)

for lam in [0.0, 0.9, 2.5, 5.5, 7.5]:
    beta, _ = lasso_ccd(y, x, lambda_=lam, n_iters=200)
    rss = np.mean((y - x @ beta) ** 2)
    print(f"lambda={lam}: beta={np.round(beta,4)} rss={round(rss,4)}")

Principal component analysis of a 3-asset correlation matrix — stats/pca1.py

"""Translated from Examples/stats/pca1.m -- PCA on a 3-asset correlation
matrix."""

import numpy as np

from quanttoolbox.linalg.special_matrices import xpnd
from quanttoolbox.stats.regression.ols import pca

C = xpnd(np.array([1.00, 0.80, 1.00, 0.80, 0.80, 1.00]), method=1)
result = pca(C)
print("eigenvalues:", np.round(result.eigenvalues, 4))
print("quality (variance share):", np.round(result.quality, 4))
print("loadings:\n", np.round(result.loadings, 4))

stats.regression.ridge

Python alternatives

Hybrid: sklearn.linear_model.Ridge/RidgeCV is better optimized for large/sparse problems. Keep ridge_tau_targeted — its L2-norm-budget (not penalty) parameterization has no sklearn equivalent.

quanttoolbox.stats.regression.ridge

Ridge regression: fixed-lambda and tau-targeted (L2-budget) variants.

Ported from QuantToolBox/stats/{regRidge,regRidge2}.m

Translation notes:

• Ridge regression's closed-form solution is a simple 3-line NumPy expression, so no external library (not even scikit-learn) is needed here -- numpy.linalg is sufficient and keeps this dependency-free.
• ridge_tau_targeted (the original regRidge2.m) reproduces the original's grid-search approach (scan a dense lambda grid, pick the closest match for the desired L2-norm budget tau) rather than solving for lambda analytically, to preserve exact behavioral parity.

ridge(y, x, lambda_)

Ridge regression at one or more penalty values lambda.

Original: stats/regRidge.m

Returns:

Name	Type	Description
beta	(n_lambda, n_vars) array of coefficients (or (n_vars,) if scalar lambda_).
df	effective degrees of freedom at each lambda.
complexity	1 / df.

Source code in src/quanttoolbox/stats/regression/ridge.py

def ridge(
    y: np.ndarray, x: np.ndarray, lambda_: np.ndarray | float
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Ridge regression at one or more penalty values lambda.

    Original: stats/regRidge.m

    Returns
    -------
    beta : (n_lambda, n_vars) array of coefficients (or (n_vars,) if scalar lambda_).
    df : effective degrees of freedom at each lambda.
    complexity : 1 / df.
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    lambdas = np.atleast_1d(np.asarray(lambda_, dtype=float))

    xx = x.T @ x
    xy = x.T @ y
    m = x.shape[1]
    identity = np.eye(m)

    n_lambda = lambdas.shape[0]
    beta = np.zeros((m, n_lambda))
    df = np.zeros(n_lambda)
    for i, lam in enumerate(lambdas):
        xx_inv = np.linalg.inv(xx + lam * identity)
        beta[:, i] = xx_inv @ xy
        df[i] = np.sum(np.diag(xx_inv @ xx))

    complexity = 1.0 / df
    if n_lambda == 1:
        return beta[:, 0], df[0], complexity[0]
    return beta.T, df, complexity

ridge_tau_targeted(y, x, tau, lambda_search=None, relative=False)

Ridge regression targeting a desired sum-of-squared-coefficients budget tau, found by grid search over a candidate lambda range.

Original: stats/regRidge2.m

Returns:

Name	Type	Description
beta	(n_tau, n_vars) array of coefficients (or (n_vars,) if scalar tau).
lambda_out	the lambda achieving the closest match to each tau.
df	effective degrees of freedom at each match.
complexity	1 / df.

Source code in src/quanttoolbox/stats/regression/ridge.py

def ridge_tau_targeted(
    y: np.ndarray,
    x: np.ndarray,
    tau: np.ndarray | float,
    lambda_search: np.ndarray | None = None,
    relative: bool = False,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Ridge regression targeting a desired sum-of-squared-coefficients
    budget tau, found by grid search over a candidate lambda range.

    Original: stats/regRidge2.m

    Returns
    -------
    beta : (n_tau, n_vars) array of coefficients (or (n_vars,) if scalar tau).
    lambda_out : the lambda achieving the closest match to each tau.
    df : effective degrees of freedom at each match.
    complexity : 1 / df.
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)

    xx = x.T @ x
    xy = x.T @ y
    m = x.shape[1]
    identity = np.eye(m)

    all_lambda = np.arange(0, 100, 0.01) if lambda_search is None else np.atleast_1d(lambda_search)

    if relative:
        beta_ols = np.linalg.inv(xx) @ xy
        tau_ols = np.sum(beta_ols**2)
    else:
        tau_ols = 1.0

    n_lambda = all_lambda.shape[0]
    all_beta = np.zeros((m, n_lambda))
    all_df = np.zeros(n_lambda)
    for i, lam in enumerate(all_lambda):
        xx_inv = np.linalg.inv(xx + lam * identity)
        all_beta[:, i] = xx_inv @ xy
        all_df[i] = np.sum(np.diag(xx_inv @ xx))

    all_tau = np.sum(all_beta**2, axis=0) / tau_ols

    tau_arr = np.atleast_1d(np.asarray(tau, dtype=float))
    n_tau = tau_arr.shape[0]
    beta = np.zeros((m, n_tau))
    df = np.zeros(n_tau)
    lambda_out = np.zeros(n_tau)
    for i, t in enumerate(tau_arr):
        idx = int(np.argmin(np.abs(all_tau - t)))
        beta[:, i] = all_beta[:, idx]
        df[i] = all_df[idx]
        lambda_out[i] = all_lambda[idx]

    complexity = 1.0 / df
    if n_tau == 1:
        return beta[:, 0], lambda_out[0], df[0], complexity[0]
    return beta.T, lambda_out, df, complexity

Examples

Norm-budget (tau-targeted) ridge regression vs. fixed-lambda ridge — stats/ridge2.py

"""Translated from Examples/stats/ridge2.m -- tau-targeted ridge regression
(`ridge_tau_targeted`) cross-checked against fixed-lambda ridge (`ridge`),
in both absolute and OLS-relative tau modes, on the same 15-observation
dataset used throughout stats/."""

import numpy as np

from quanttoolbox.stats.regression.ols import standardize
from quanttoolbox.stats.regression.ridge import ridge, ridge_tau_targeted

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])

beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)
tau_ols = np.sum(beta_ols**2)
print("tau(ols) =", round(tau_ols, 4))

# `ridge()` returns (beta, df, complexity) -- it doesn't report the
# resulting L2 budget, so tau_ridge (the quantity `ridge_tau_targeted`
# was aiming for) is recomputed here as sum(beta**2), to cross-check that
# fixed-lambda ridge at `lambda_out` reproduces the same beta/tau
# `ridge_tau_targeted` found by its grid search.

# Absolute tau targets
tau = np.array([0.7, 0.9, 1.1])
beta_ridge2, lambda_out, df_ridge2, complexity2 = ridge_tau_targeted(y, x, tau)
beta_ridge, df_ridge, complexity = ridge(y, x, lambda_out)
tau_ridge = np.sum(beta_ridge**2, axis=1)

print("\nAbsolute-tau analysis: tau target, tau achieved (fixed-lambda ridge), lambda")
print(np.round(np.column_stack([tau, tau_ridge, lambda_out]), 4))

# Relative tau targets (fraction of tau_ols)
tau = np.array([0.25, 0.50, 0.75, 1.00])
beta_ridge2, lambda_out, df_ridge2, complexity2 = ridge_tau_targeted(y, x, tau, relative=True)
beta_ridge, df_ridge, complexity = ridge(y, x, lambda_out)
tau_ridge_rel = np.sum(beta_ridge**2, axis=1) / tau_ols

print("\nRelative-tau analysis: tau target, tau achieved / tau_ols, lambda")
print(np.round(np.column_stack([tau, tau_ridge_rel, lambda_out]), 4))
print("\nbeta (fixed-lambda ridge) at each relative-tau lambda:")
print(np.round(beta_ridge, 4))
print("\nbeta (tau-targeted ridge):")
print(np.round(beta_ridge2, 4))

# Finer lambda search grid
lambda_search = np.arange(0, 15, 0.001)
beta_ridge2, lambda_out, df_ridge2, complexity2 = ridge_tau_targeted(
    y, x, tau, lambda_search=lambda_search, relative=True
)
beta_ridge, df_ridge, complexity = ridge(y, x, lambda_out)
tau_ridge_rel = np.sum(beta_ridge**2, axis=1) / tau_ols

print("\nRelative-tau analysis, finer lambda grid: tau target, tau achieved / tau_ols, lambda")
print(np.round(np.column_stack([tau, tau_ridge_rel, lambda_out]), 4))

stats.regression.lasso

Python alternatives

Switch the penalized-form solvers (lasso_ccd, lasso_admm) to sklearn.linear_model.Lasso/ElasticNet — Cython-compiled, extensively battle-tested. Keep lasso_tau_constrained — sklearn has no L1-budget interface.

quanttoolbox.stats.regression.lasso

Lasso and elastic-net regression via coordinate descent / ADMM.

Ported from QuantToolBox/stats/{regLassoCCD,regLassoADMM,regLassoADMM2, regLasso,regElasticNet,selectLasso,regLasso2}.m

Translation notes:

• The original MATLAB toolbox has two Lasso parameterizations that are easy to conflate:
  1. Penalized form (regLassoCCD.m, regLassoADMM.m): minimize ||y - X*beta||^2 + lambda * ||beta||_1. This is what coordinate descent and ADMM solve directly, and what scikit-learn's own Lasso (also coordinate-descent-based) solves too.
  2. Constrained/budget form (regLasso.m, selectLasso.m, via quadprog): minimize ||y - X*beta||^2 subject to ||beta||_1 <= tau. These are Lagrangian duals of each other but need different solvers. This module ports the penalized form directly (lasso_ccd, lasso_admm) using plain NumPy coordinate descent -- the same algorithm scikit-learn uses internally, so no extra dependency is needed. lasso_tau_constrained recovers the budget-form behavior of regLasso.m/selectLasso.m by bisecting on lambda until the resulting L1 norm matches the target tau, rather than porting the original's quadprog-based QP formulation (MATLAB Optimization Toolbox has no free-standing Python equivalent worth adding as a dependency just for this).
• regElasticNet.m's alpha/lambda parameterization (lambda*(1-alpha)*L2 + lambda*alpha*L1) is preserved in elastic_net_ccd for call-site compatibility, rather than switching to scikit-learn's l1_ratio naming.
• MATLAB's global ADMM_* settings blocks are replaced by the quanttoolbox.config.ADMMConfig dataclass.
• regLasso2.m builds a regularization path by looping its solver over a vector of lambdas, returning [beta, tau, df, complexity]. It uses quadprog on a doubled positive/negative-part formulation of beta (like regLasso.m/selectLasso.m); this is ported as lasso_regularization_path, which loops lasso_ccd instead -- the same quadprog-avoiding substitution docs/examples/stats/lasso2.py already makes for Examples/stats/lasso2.m's single-lambda regLasso2 calls.

elastic_net_ccd(y, x, lambda_, alpha, beta_init=None, n_iters=100)

Elastic net via cyclical coordinate descent: minimize ||y - Xbeta||^2 + lambda(1-alpha)||beta||_2^2 + lambdaalpha*||beta||_1.

Original: stats/regElasticNet.m (ported as CCD rather than the original's quadprog QP formulation -- see module docstring.)

Source code in src/quanttoolbox/stats/regression/lasso.py

def elastic_net_ccd(
    y: np.ndarray,
    x: np.ndarray,
    lambda_: float,
    alpha: float,
    beta_init: np.ndarray | None = None,
    n_iters: int = 100,
) -> tuple[np.ndarray, np.ndarray]:
    """Elastic net via cyclical coordinate descent:
    minimize ||y - X*beta||^2 + lambda*(1-alpha)*||beta||_2^2 + lambda*alpha*||beta||_1.

    Original: stats/regElasticNet.m (ported as CCD rather than the
    original's quadprog QP formulation -- see module docstring.)
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    n, p = x.shape

    beta = (
        np.linalg.inv(x.T @ x) @ (x.T @ y)
        if beta_init is None or np.asarray(beta_init).shape[0] != p
        else np.asarray(beta_init, dtype=float).copy()
    )

    l1_penalty = lambda_ * alpha
    l2_penalty = lambda_ * (1 - alpha)

    beta_path = np.zeros((p, n_iters))
    for it in range(n_iters):
        beta_path[:, it] = beta
        for j in range(p):
            x_j = x[:, j]
            x_minus_j = x.copy()
            x_minus_j[:, j] = 0.0
            v = x_j @ (y - x_minus_j @ beta)
            denom = x_j @ x_j + l2_penalty
            beta[j] = np.sign(v) * max(abs(v) - l1_penalty, 0.0) / denom

    return beta, beta_path.T

lasso_admm(y, x, lambda_, beta_init=None, config=None)

Lasso via ADMM (penalized form) with a fixed step-size (varphi).

Note: the original MATLAB also supported an adaptive-varphi mode (ADMM_varphi_mtd == 2, residual-balancing per Boyd et al. 2011); that branch is not ported here -- only the fixed-varphi path, which is what the default config uses anyway. Add adaptive step-size updates here if convergence speed on your problem needs it.

Original: stats/regLassoADMM.m

Returns:

Name	Type	Description
beta	final coefficient estimate.
beta_path	(n_iters_run, n_vars) iterate history.
converged	whether the convergence tolerance was reached.
n_iters_run	number of iterations actually performed.

Source code in src/quanttoolbox/stats/regression/lasso.py

def lasso_admm(
    y: np.ndarray,
    x: np.ndarray,
    lambda_: float,
    beta_init: np.ndarray | None = None,
    config: ADMMConfig | None = None,
) -> tuple[np.ndarray, np.ndarray, bool, int]:
    """Lasso via ADMM (penalized form) with a fixed step-size (varphi).

    Note: the original MATLAB also supported an adaptive-varphi mode
    (``ADMM_varphi_mtd == 2``, residual-balancing per Boyd et al. 2011);
    that branch is not ported here -- only the fixed-varphi path, which
    is what the default config uses anyway. Add adaptive step-size
    updates here if convergence speed on your problem needs it.

    Original: stats/regLassoADMM.m

    Returns
    -------
    beta : final coefficient estimate.
    beta_path : (n_iters_run, n_vars) iterate history.
    converged : whether the convergence tolerance was reached.
    n_iters_run : number of iterations actually performed.
    """
    if config is None:
        config = ADMMConfig()

    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    n, p = x.shape

    xx = x.T @ x
    xy = x.T @ y
    identity = np.eye(p)

    beta = (
        np.linalg.inv(xx) @ xy
        if beta_init is None or np.asarray(beta_init).shape[0] != p
        else np.asarray(beta_init, dtype=float).copy()
    )

    varphi = config.varphi
    beta_prev = beta.copy()
    beta_bar = beta.copy()
    beta_bar_prev = beta_bar.copy()
    u = np.zeros(p)

    beta_path = np.full((config.max_iters, p), np.nan)
    converged = False
    n_iters_run = config.max_iters

    for it in range(config.max_iters):
        beta_path[it] = beta

        v = beta_bar - u
        beta = np.linalg.inv(xx + varphi * identity) @ (xy + varphi * v)

        v = beta + u
        beta_bar = soft_threshold(v, lambda_ / varphi)

        r = beta - beta_bar
        u = u + r

        cvg1 = np.sum((beta - beta_prev) ** 2)
        cvg2 = np.sum(r**2)
        cvg3 = np.sum((beta_bar - beta_bar_prev) ** 2)
        cvg = max(cvg1, cvg2, cvg3)
        if cvg <= config.tol:
            converged = True
            n_iters_run = it + 1
            break

        beta_prev = beta.copy()
        beta_bar_prev = beta_bar.copy()

    return beta, beta_path[:n_iters_run], converged, n_iters_run

lasso_ccd(y, x, lambda_, beta_init=None, n_iters=100)

Lasso via cyclical coordinate descent (penalized form).

Original: stats/regLassoCCD.m

Source code in src/quanttoolbox/stats/regression/lasso.py

def lasso_ccd(
    y: np.ndarray,
    x: np.ndarray,
    lambda_: float,
    beta_init: np.ndarray | None = None,
    n_iters: int = 100,
) -> tuple[np.ndarray, np.ndarray]:
    """Lasso via cyclical coordinate descent (penalized form).

    Original: stats/regLassoCCD.m
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    n, p = x.shape

    beta = (
        np.linalg.inv(x.T @ x) @ (x.T @ y)
        if beta_init is None or np.asarray(beta_init).shape[0] != p
        else np.asarray(beta_init, dtype=float).copy()
    )

    beta_path = np.zeros((p, n_iters))
    for it in range(n_iters):
        beta_path[:, it] = beta
        for j in range(p):
            x_j = x[:, j]
            x_minus_j = x.copy()
            x_minus_j[:, j] = 0.0
            v = x_j @ (y - x_minus_j @ beta)
            denom = x_j @ x_j
            if lambda_ > 0:
                beta[j] = np.sign(v) * max(abs(v) - lambda_, 0.0) / denom
            else:
                beta[j] = v / denom

    return beta, beta_path.T

lasso_regularization_path(y, x, lambda_grid, n_iters=1000)

Build the full lasso regularization path by looping lasso_ccd (penalized-form coordinate descent) over a grid of lambda values.

At each lambda, records the fitted coefficients, their L1 norm, the number of nonzero coefficients (degrees of freedom), and the resulting model complexity (1 / df).

Original: stats/regLasso2.m -- regLasso2.m solves the same [beta, tau, df, complexity] = f(y, x, lambda) grid problem via quadprog on a doubled positive/negative-part formulation of beta; this is translated using lasso_ccd's coordinate descent instead, the same substitution docs/examples/stats/lasso2.py already makes for Examples/stats/lasso2.m's single-lambda regLasso2 calls (see module docstring). Promoted from HSF-Notebooks chapter 15b (reg_lasso2), which uses this exact grid-loop as a documented example of driving quanttoolbox's own lasso_ccd solver across a lambda grid, rather than an external one (e.g. scikit-learn's lasso_path).

Parameters:

Name	Type	Description	Default
y	response vector, shape (n,).		required
x	design matrix, shape (n, p).		required
lambda_grid	penalty values to solve at, shape (n_lambda,).		required
n_iters	coordinate-descent iterations passed through to	lasso_ccd at each lambda.	1000

Returns:

Type	Description
LassoPathResult with, for each lambda in ``lambda_grid``:	beta: coefficient matrix, shape (n_lambda, p). l1_norm: sum(abs(beta)) per row. df: count of coefficients with abs(beta) >= 1e-6 per row. complexity: 1 / df per row (+inf where every coefficient has been shrunk to zero -- expected at heavy enough penalization, not an error).

Source code in src/quanttoolbox/stats/regression/lasso.py

def lasso_regularization_path(
    y: np.ndarray,
    x: np.ndarray,
    lambda_grid: np.ndarray,
    n_iters: int = 1000,
) -> LassoPathResult:
    """Build the full lasso regularization path by looping ``lasso_ccd``
    (penalized-form coordinate descent) over a grid of lambda values.

    At each lambda, records the fitted coefficients, their L1 norm, the
    number of nonzero coefficients (degrees of freedom), and the resulting
    model complexity (``1 / df``).

    Original: stats/regLasso2.m -- ``regLasso2.m`` solves the same
    ``[beta, tau, df, complexity] = f(y, x, lambda)`` grid problem via
    ``quadprog`` on a doubled positive/negative-part formulation of beta;
    this is translated using ``lasso_ccd``'s coordinate descent instead,
    the same substitution ``docs/examples/stats/lasso2.py`` already makes
    for ``Examples/stats/lasso2.m``'s single-lambda ``regLasso2`` calls
    (see module docstring). Promoted from HSF-Notebooks chapter 15b
    (``reg_lasso2``), which uses this exact grid-loop as a documented
    example of driving quanttoolbox's own ``lasso_ccd`` solver across a
    lambda grid, rather than an external one (e.g. scikit-learn's
    ``lasso_path``).

    Parameters
    ----------
    y : response vector, shape (n,).
    x : design matrix, shape (n, p).
    lambda_grid : penalty values to solve at, shape (n_lambda,).
    n_iters : coordinate-descent iterations passed through to
        ``lasso_ccd`` at each lambda.

    Returns
    -------
    LassoPathResult with, for each lambda in ``lambda_grid``:

        - ``beta``: coefficient matrix, shape (n_lambda, p).
        - ``l1_norm``: ``sum(abs(beta))`` per row.
        - ``df``: count of coefficients with ``abs(beta) >= 1e-6`` per row.
        - ``complexity``: ``1 / df`` per row (``+inf`` where every
          coefficient has been shrunk to zero -- expected at heavy enough
          penalization, not an error).
    """
    x = np.asarray(x, dtype=float)
    lambda_grid = np.atleast_1d(np.asarray(lambda_grid, dtype=float))
    n_lambda = len(lambda_grid)
    p = x.shape[1]

    beta = np.zeros((n_lambda, p))
    for i, lam in enumerate(lambda_grid):
        b, _ = lasso_ccd(y, x, lam, n_iters=n_iters)
        beta[i, :] = b

    l1_norm = np.sum(np.abs(beta), axis=1)
    df = np.sum(np.abs(beta) >= 1e-6, axis=1)
    # At heavy enough shrinkage all coefficients are zeroed (df=0); complexity=1/df is
    # then +inf at that boundary case, which is expected, not an error.
    with np.errstate(divide="ignore"):
        complexity = 1.0 / df

    return LassoPathResult(
        lambda_grid=lambda_grid, beta=beta, l1_norm=l1_norm, df=df, complexity=complexity
    )

lasso_tau_constrained(y, x, tau, lambda_search=None, n_iters=200)

Lasso in budget-constrained form: minimize ||y - X*beta||^2 subject to ||beta||_1 <= tau, found via bisection on the penalized-form lambda.

Reproduces the tau-based interface of stats/regLasso.m and stats/selectLasso.m without needing a QP solver.

Returns:

Name	Type	Description
beta	coefficients whose L1 norm is closest to tau.
lambda_matched	the lambda value that achieved it.
df	number of non-zero coefficients.

Source code in src/quanttoolbox/stats/regression/lasso.py

def lasso_tau_constrained(
    y: np.ndarray,
    x: np.ndarray,
    tau: float,
    lambda_search: np.ndarray | None = None,
    n_iters: int = 200,
) -> tuple[np.ndarray, float, int]:
    """Lasso in budget-constrained form: minimize ||y - X*beta||^2 subject to
    ||beta||_1 <= tau, found via bisection on the penalized-form lambda.

    Reproduces the tau-based interface of stats/regLasso.m and
    stats/selectLasso.m without needing a QP solver.

    Returns
    -------
    beta : coefficients whose L1 norm is closest to tau.
    lambda_matched : the lambda value that achieved it.
    df : number of non-zero coefficients.
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)

    beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)
    if np.sum(np.abs(beta_ols)) <= tau:
        # unconstrained OLS already satisfies the budget
        df = int(np.sum(np.abs(beta_ols) >= 1e-10))
        return beta_ols, 0.0, df

    lo, hi = 0.0, np.max(np.abs(x.T @ y)) * 2  # hi large enough to drive beta to 0
    beta = beta_ols
    for _ in range(n_iters):
        mid = 0.5 * (lo + hi)
        beta, _ = lasso_ccd(y, x, mid, beta_init=beta.copy(), n_iters=50)
        l1_norm = np.sum(np.abs(beta))
        if l1_norm > tau:
            lo = mid
        else:
            hi = mid
        if abs(l1_norm - tau) < 1e-6:
            break

    df = int(np.sum(np.abs(beta) >= 1e-10))
    return beta, mid, df

soft_threshold(v, threshold)

Elementwise soft-thresholding operator: sign(v) * max(|v| - threshold, 0).

Original: optim/soft_thresholding.m (used by regLassoADMM.m)

Source code in src/quanttoolbox/stats/regression/lasso.py

def soft_threshold(v: np.ndarray, threshold: float) -> np.ndarray:
    """Elementwise soft-thresholding operator: sign(v) * max(|v| - threshold, 0).

    Original: optim/soft_thresholding.m (used by regLassoADMM.m)
    """
    return np.sign(v) * np.maximum(np.abs(v) - threshold, 0.0)

Examples

Elastic-net regression path — stats/elasticnet1.py

"""Translated from Examples/stats/elasticnet1.m and elasticnet2.m --
elastic net regression path at two different alpha (L1/L2 mix) values,
on the same 15-observation dataset used in ridge1.m/lasso1.m."""

import numpy as np

from quanttoolbox.stats.regression.lasso import elastic_net_ccd
from quanttoolbox.stats.regression.ols import standardize

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])

for alpha, lam in [(0.50, 1.0), (0.25, 1.0)]:
    beta, path = elastic_net_ccd(y, x, lambda_=lam, alpha=alpha, n_iters=200)
    print(f"alpha={alpha}, lambda={lam}: beta={np.round(beta,4)}")

Lasso path with R-squared, degrees of freedom, and complexity — stats/lasso1.py

"""Translated from Examples/stats/lasso1.m -- penalized-form lasso path at
5 lambda values (numeric core only; the original's 501-point tau-sweep and
its multi-line coefficient-path plot are dropped -- the same lambda
values and dataset are already used by lasso2.py, which this extends with
R^2/degrees-of-freedom/complexity reporting to match lasso1.m's own
output table)."""

import numpy as np

from quanttoolbox.stats.regression.lasso import lasso_ccd
from quanttoolbox.stats.regression.ols import standardize

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])

beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)
tss = np.mean(y**2)

print("lambda  |beta|         tau=sum|beta| rss     R2      df  complexity")
for lam in [0.0, 0.9, 2.5, 5.5, 7.5]:
    beta, _ = lasso_ccd(y, x, lambda_=lam, n_iters=200)
    u = y - x @ beta
    rss = np.mean(u**2)
    r2 = 1 - rss / tss
    tau = np.sum(np.abs(beta))
    df = int(np.sum(np.abs(beta) >= 1e-6))
    complexity = 1.0 / df if df > 0 else np.inf
    print(
        f"{lam:6.2f}  {np.round(beta, 4)}  tau={tau:.4f}  rss={rss:.4f}  "
        f"R2={r2:.4f}  df={df}  complexity={complexity:.4f}"
    )

Penalized-form lasso path via coordinate descent — stats/lasso2.py

"""Translated from Examples/stats/lasso2.m -- penalized-form lasso
(regLasso2) at several lambda values, on the same 15-obs dataset."""

import numpy as np

from quanttoolbox.stats.regression.lasso import lasso_ccd
from quanttoolbox.stats.regression.ols import standardize

data = np.array(
    [
        [3.1, 2.8, 4.3, 0.3, 2.2, 3.5],
        [24.9, 5.9, 3.6, 3.2, 0.7, 6.4],
        [27.3, 6.0, 9.6, 7.6, 9.5, 0.9],
        [25.4, 8.4, 5.4, 1.8, 1.0, 7.1],
        [46.1, 5.2, 7.6, 8.3, 0.6, 4.5],
        [45.7, 6.0, 7.0, 9.6, 0.6, 0.6],
        [47.4, 6.1, 1.0, 8.5, 9.6, 8.6],
        [-1.8, 1.2, 9.6, 2.7, 4.8, 5.8],
        [20.8, 3.2, 5.0, 4.2, 2.7, 3.6],
        [6.8, 0.5, 9.2, 6.9, 9.3, 0.7],
        [12.9, 7.9, 9.1, 1.0, 5.9, 5.4],
        [37.0, 1.8, 1.3, 9.2, 6.1, 8.3],
        [14.7, 7.4, 5.6, 0.9, 5.6, 3.9],
        [-3.2, 2.3, 6.6, 0.0, 3.6, 6.4],
        [44.3, 7.7, 2.2, 6.5, 1.3, 0.7],
    ]
)
y = standardize(data[:, 0])
x = standardize(data[:, 1:6])
beta_ols = np.linalg.inv(x.T @ x) @ (x.T @ y)

for lam in [0.0, 0.9, 2.5, 5.5, 7.5]:
    beta, _ = lasso_ccd(y, x, lambda_=lam, n_iters=200)
    rss = np.mean((y - x @ beta) ** 2)
    print(f"lambda={lam}: beta={np.round(beta,4)} rss={round(rss,4)}")

stats.regression.kernel

Python alternatives

Hybrid: statsmodels.nonparametric.KernelReg supports automatic bandwidth selection via cross-validation and both local-constant/local-linear estimators. Switch for general use; keep ours where the exact original bandwidth formula must match existing MATLAB-based results.

quanttoolbox.stats.regression.kernel

Nonparametric kernel density estimation and local-polynomial kernel regression.

Ported from QuantToolBox/stats/{regKernelDensity,regKernelMean, regKernelQuantile,regKernelPayoff}.m

Translation notes:

• regKernelDensity is a manual Gaussian-kernel density/CDF estimator using scipy.stats.norm as the kernel -- this is equivalent to scipy.stats.gaussian_kde but the original's Silverman-type bandwidth rule (1.364 * std * n^-0.2) differs slightly from scipy's default bandwidth factor, so it's reproduced explicitly here rather than delegating to gaussian_kde (which would silently change bandwidths).
• regKernelMean/regKernelPayoff are local polynomial regression (Nadaraya-Watson generalized to include local polynomial terms), solved via a small per-evaluation-point weighted least squares -- there's no single scipy/sklearn function that matches this directly, so it's ported as a direct weighted-least-squares loop.
• regKernelQuantile depends on quantile_regression (scipy.optimize.linprog-based); see quanttoolbox.stats.regression.quantile.

kernel_density(data, x, h=None)

Gaussian-kernel density (and CDF) estimate of each column of data, evaluated at points x.

Original: stats/regKernelDensity.m

Returns:

Name	Type	Description
density	(n_x, n_cols) estimated density values.
cdf	(n_x, n_cols) estimated CDF values.
bandwidths	(n_cols,) bandwidth used per column (auto-computed if h is None).

Source code in src/quanttoolbox/stats/regression/kernel.py

def kernel_density(
    data: np.ndarray, x: np.ndarray, h: float | np.ndarray | None = None
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Gaussian-kernel density (and CDF) estimate of each column of `data`,
    evaluated at points `x`.

    Original: stats/regKernelDensity.m

    Returns
    -------
    density : (n_x, n_cols) estimated density values.
    cdf : (n_x, n_cols) estimated CDF values.
    bandwidths : (n_cols,) bandwidth used per column (auto-computed if h is None).
    """
    data = np.atleast_2d(np.asarray(data, dtype=float))
    if data.shape[0] == 1:
        data = data.T
    n_cols = data.shape[1]

    x = np.atleast_2d(np.asarray(x, dtype=float))
    if x.shape[0] == 1 and x.shape[1] != n_cols:
        x = x.T
    if x.shape[1] == 1 and n_cols > 1:
        x = np.tile(x, (1, n_cols))

    bandwidths = np.zeros(n_cols) if h is None else np.broadcast_to(h, (n_cols,)).astype(float)

    n_x = x.shape[0]
    density = np.zeros((n_x, n_cols))
    cdf = np.zeros((n_x, n_cols))

    for col in range(n_cols):
        y = data[:, col]
        y = y[~np.isnan(y)]
        bw = bandwidths[col]
        if bw == 0:
            bw = _silverman_bandwidth(y)
            bandwidths[col] = bw

        for i in range(n_x):
            u = (x[i, col] - y) / bw
            density[i, col] = np.mean(norm.pdf(u)) / bw
            cdf[i, col] = np.mean(norm.cdf(u))

    return density, cdf, bandwidths

kernel_mean_regression(y, x, z, order=1, h=None)

Local polynomial kernel regression of y on x, evaluated at points z.

order=1 is standard local-linear (Nadaraya-Watson-generalized) regression; higher orders fit a local polynomial of that degree.

Original: stats/regKernelMean.m

Source code in src/quanttoolbox/stats/regression/kernel.py

def kernel_mean_regression(
    y: np.ndarray, x: np.ndarray, z: np.ndarray, order: int = 1, h: float | None = None
) -> np.ndarray:
    """Local polynomial kernel regression of y on x, evaluated at points z.

    order=1 is standard local-linear (Nadaraya-Watson-generalized)
    regression; higher orders fit a local polynomial of that degree.

    Original: stats/regKernelMean.m
    """
    x = np.asarray(x, dtype=float).flatten()
    y = np.asarray(y, dtype=float).flatten()
    valid = ~np.isnan(x) & ~np.isnan(y)
    x, y = x[valid], y[valid]
    n_x = x.shape[0]

    if h is None:
        h = 1.364 * x.std(ddof=1) * n_x ** (-0.20)

    z = np.atleast_1d(np.asarray(z, dtype=float))
    n_z = z.shape[0]
    m = np.zeros(n_z)

    design = np.ones((n_x, order + 1))
    for i in range(n_z):
        dz = x - z[i]
        w = norm.pdf(dz / h)
        for k in range(1, order + 1):
            design[:, k] = dz**k
        weighted_design = design * w[:, None]
        beta = np.linalg.solve(weighted_design.T @ design, weighted_design.T @ y)
        m[i] = beta[0]

    return m

kernel_payoff_regression(x, y, z=None, n_points=100, order=1, h=None)

Convenience wrapper around kernel_mean_regression: builds an evaluation grid z automatically (from data range, or from a (min, max) pair) if one isn't supplied directly.

Original: stats/regKernelPayoff.m

Source code in src/quanttoolbox/stats/regression/kernel.py

def kernel_payoff_regression(
    x: np.ndarray,
    y: np.ndarray,
    z: np.ndarray | tuple[float, float] | None = None,
    n_points: int = 100,
    order: int = 1,
    h: float | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Convenience wrapper around kernel_mean_regression: builds an
    evaluation grid `z` automatically (from data range, or from a
    (min, max) pair) if one isn't supplied directly.

    Original: stats/regKernelPayoff.m
    """
    x = np.asarray(x, dtype=float).flatten()
    y = np.asarray(y, dtype=float).flatten()
    valid = ~np.isnan(x) & ~np.isnan(y)
    x, y = x[valid], y[valid]

    if z is None:
        rg_min, rg_max = x.min(), x.max()
        z = np.linspace(rg_min, rg_max, n_points)
    elif np.asarray(z).shape[0] == 2:
        rg_min, rg_max = z
        z = np.linspace(rg_min, rg_max, n_points)
    else:
        z = np.asarray(z, dtype=float)

    q = kernel_mean_regression(y, x, z, order=order, h=h)
    return z, q

kernel_quantile_regression(y, x, tau, z, order=1, h=1.0)

Local polynomial kernel quantile regression of y on x at quantile tau, evaluated at points z.

Original: stats/regKernelQuantile.m

Source code in src/quanttoolbox/stats/regression/kernel.py

def kernel_quantile_regression(
    y: np.ndarray,
    x: np.ndarray,
    tau: float,
    z: np.ndarray,
    order: int = 1,
    h: float = 1.0,
) -> np.ndarray:
    """Local polynomial kernel quantile regression of y on x at quantile tau,
    evaluated at points z.

    Original: stats/regKernelQuantile.m
    """
    x = np.asarray(x, dtype=float).flatten()
    y = np.asarray(y, dtype=float).flatten()
    n_x = x.shape[0]

    h1 = 1.364 * x.std(ddof=1) * n_x ** (-0.20)
    h2 = (tau * (1 - tau) * norm.pdf(norm.ppf(tau)) ** (-2)) ** 0.20
    bandwidth = h * h1 * h2

    z = np.atleast_1d(np.asarray(z, dtype=float))
    n_z = z.shape[0]
    q = np.zeros(n_z)

    const = np.ones(n_x)
    for i in range(n_z):
        dz = z[i] - x
        w = norm.pdf(dz / bandwidth)
        design = np.column_stack([const] + [dz**k for k in range(1, order + 1)])
        beta, _, _ = quantile_regression(y, design, tau, weights=w)
        q[i] = beta[0]

    return q

Examples

Gaussian kernel density of two correlated series — stats/kernel1.py

"""Translated from Examples/stats/kernel1.m -- Gaussian-kernel density
estimate of two correlated series (numeric core only; the original's plot
of both density curves is dropped).

The original draws x from MATLAB's unseeded `randn`; a fixed seed
(`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.stats.regression.kernel import kernel_density

rng = np.random.default_rng(0)
x1 = rng.standard_normal(100)
x2 = 0.8 * x1 + 1.0
x = np.column_stack([x1, x2])

z = np.linspace(-5, 5, 101)

density, cdf, bandwidths = kernel_density(x, z)

print("Bandwidths (auto-computed, Silverman-type rule):", np.round(bandwidths, 4))
print("\nDensity at z in {-5,-2.5,0,2.5,5} (columns: series 1, series 2):")
idx = [0, 25, 50, 75, 100]
print(np.round(np.column_stack([z[idx], density[idx]]), 4))

Kernel mean and quantile regression vs. true population curves — stats/qreg2.py

"""Translated from Examples/stats/qreg2.m -- local-polynomial kernel mean
and quantile regression (orders 1 and 2) of a nonlinear y=f(x)+noise
relationship, compared against the known population mean/quantile curves
(numeric core only; the original's scatter+curve plot is dropped).

The original explicitly seeds MATLAB's RNG (`rng(123)`); NumPy's generator
is seeded the same way for a comparable (not bit-identical) run."""

import numpy as np

from quanttoolbox.stats.regression.kernel import kernel_mean_regression, kernel_quantile_regression

rng = np.random.default_rng(123)
n = 500
x = rng.random(n)
y = rng.random(n) * (np.cos(2 * np.pi * x - np.pi) + 1)

tau = 0.95
p = 50
z = np.arange(0, 1 + 1e-9, 1 / (p - 1))[:p]

m0 = 0.5 * (np.cos(2 * np.pi * z - np.pi) + 1)
m1 = kernel_mean_regression(y, x, z, order=1)
m2 = kernel_mean_regression(y, x, z, order=2)

q0 = tau * (np.cos(2 * np.pi * z - np.pi) + 1)
q1 = kernel_quantile_regression(y, x, tau, z, order=1)
q2 = kernel_quantile_regression(y, x, tau, z, order=2)

print("z, population mean m0, local-linear m1, local-quadratic m2 (every 10th point):")
print(np.round(np.column_stack([z, m0, m1, m2])[::10], 4))
print("\nz, population q(0.95) q0, local-linear q1, local-quadratic q2 (every 10th point):")
print(np.round(np.column_stack([z, q0, q1, q2])[::10], 4))

stats.regression.quantile

Python alternatives

Switch to statsmodels.regression.quantile_regression.QuantReg or sklearn.linear_model.QuantileRegressor — both more battle-tested than this module's scipy.optimize.linprog-based implementation. Keep ours only if the LP slack-variable (u, v) outputs are needed downstream.

quanttoolbox.stats.regression.quantile

Quantile regression (via linear programming) and quantile-regression copulas.

Ported from QuantToolBox/stats/{quantile_regression,qrCopulaNormal, qrCopulaStudent}.m

Translation notes:

• MATLAB's linprog (Optimization Toolbox, interior-point) maps directly onto scipy.optimize.linprog (also supports an interior-point method) -- both are standard-library-adjacent linear programming solvers, so no extra dependency is needed.
• The original solves each quantile tau independently in a loop; this is preserved here (rather than vectorizing across tau) since each is an independent LP.

qr_copula_normal(u1, rho, alpha)

Conditional quantile u2 = Q(alpha | u1) implied by a bivariate Gaussian copula with correlation rho.

Original: stats/qrCopulaNormal.m

Source code in src/quanttoolbox/stats/regression/quantile.py

def qr_copula_normal(
    u1: float | np.ndarray, rho: float, alpha: float | np.ndarray
) -> float | np.ndarray:
    """Conditional quantile u2 = Q(alpha | u1) implied by a bivariate Gaussian
    copula with correlation rho.

    Original: stats/qrCopulaNormal.m
    """
    return normal_cdf(rho * normal_ppf(u1) + np.sqrt(1 - rho**2) * normal_ppf(alpha))

qr_copula_student(u1, rho, nu, alpha)

Conditional quantile u2 = Q(alpha | u1) implied by a bivariate Student-t copula with correlation rho and nu degrees of freedom.

Original: stats/qrCopulaStudent.m

Source code in src/quanttoolbox/stats/regression/quantile.py

def qr_copula_student(
    u1: float | np.ndarray, rho: float, nu: float, alpha: float | np.ndarray
) -> float | np.ndarray:
    """Conditional quantile u2 = Q(alpha | u1) implied by a bivariate Student-t
    copula with correlation rho and nu degrees of freedom.

    Original: stats/qrCopulaStudent.m
    """
    t1 = student_t_ppf(u1, nu)
    t2 = rho * t1 + np.sqrt((1 - rho**2) * (nu + t1**2) / (1 + nu)) * student_t_ppf(alpha, nu + 1)
    return student_t_cdf(t2, nu)

quantile_regression(y, x, tau, weights=None)

Linear quantile regression at one or more quantile levels tau, solved via linear programming (Koenker & Bassett's LP formulation).

Original: stats/quantile_regression.m

Returns:

Name	Type	Description
beta	(n_vars, n_tau) coefficients (or (n_vars,) if scalar tau).
u	(n_obs, n_tau) positive residual parts.
v	(n_obs, n_tau) negative residual parts.

Source code in src/quanttoolbox/stats/regression/quantile.py

def quantile_regression(
    y: np.ndarray, x: np.ndarray, tau: float | np.ndarray, weights: np.ndarray | None = None
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Linear quantile regression at one or more quantile levels tau, solved
    via linear programming (Koenker & Bassett's LP formulation).

    Original: stats/quantile_regression.m

    Returns
    -------
    beta : (n_vars, n_tau) coefficients (or (n_vars,) if scalar tau).
    u : (n_obs, n_tau) positive residual parts.
    v : (n_obs, n_tau) negative residual parts.
    """
    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    n, m = x.shape
    tau_arr = np.atleast_1d(np.asarray(tau, dtype=float))
    p = tau_arr.shape[0]

    w = np.ones(n) if weights is None else np.asarray(weights, dtype=float)
    y_w = y * w
    x_w = x * w[:, None]

    a_eq = np.hstack([x_w, np.eye(n), -np.eye(n)])
    b_eq = y_w
    bounds = [(None, None)] * m + [(0, None)] * n + [(0, None)] * n

    beta = np.full((m, p), np.nan)
    u = np.full((n, p), np.nan)
    v = np.full((n, p), np.nan)

    for i, t in enumerate(tau_arr):
        c = np.concatenate([np.zeros(m), t * np.ones(n), (1 - t) * np.ones(n)])
        result = linprog(c, A_eq=a_eq, b_eq=b_eq, bounds=bounds, method="highs")
        if result.success:
            z = result.x
            beta[:, i] = z[:m]
            u[:, i] = z[m : m + n]
            v[:, i] = z[m + n :]

    if p == 1:
        return beta[:, 0], u[:, 0], v[:, 0]
    return beta, u, v

Examples

Linear quantile regression at nine quantile levels — stats/qreg1.py

"""Translated from Examples/stats/qreg1.m -- linear quantile regression at
9 quantile levels (tau=0.1..0.9) on a simulated 3-predictor dataset.

The original draws x/b/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.stats.regression.quantile import quantile_regression

rng = np.random.default_rng(0)
n, m = 100, 3
tau = np.arange(0.1, 0.91, 0.1)
x = 10 * rng.random((n, m))
b = 10 * rng.standard_normal(m)
sigma = 5.20
u = sigma * rng.standard_normal(n)
y = x @ b + u

beta, u_pos, v_neg = quantile_regression(y, x, tau)
res = y[:, None] - x @ beta

print("tau:", np.round(tau, 2))
print("beta (rows=predictors, cols=tau):")
print(np.round(beta, 4))
print("\nresidual std by tau:", np.round(np.std(res, axis=0), 4))

Monte Carlo OLS vs. LAD regression under heteroskedastic noise — ects/quantile2.py

"""Translated from Examples/ects/quantile2.m -- Monte Carlo comparison of
OLS vs. LAD (median quantile regression) under heteroskedastic noise:
mean and standard error of each estimator's slope coefficient across many
simulated datasets (numeric core only; the original's kernel-density plot
of the two estimators' sampling distributions is dropped).

The original draws x/sigma/u from MATLAB's unseeded `rand`/`randn`; a
fixed seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.stats.regression.quantile import quantile_regression

rng = np.random.default_rng(0)
n_s = 2000
n_t = 200

b0, b1 = 0.1, 0.2
alpha = 0.50

beta_ols = np.zeros(n_s)
beta_lad = np.zeros(n_s)

for it in range(n_s):
    x = rng.standard_normal(n_t)
    sigma = 0.20 + 0.60 * rng.random(n_t)
    sigma = sigma**2.5
    u = sigma * rng.standard_normal(n_t)
    y = b0 + b1 * x + u

    design = np.column_stack([np.ones(n_t), x])
    beta = np.linalg.inv(design.T @ design) @ (design.T @ y)
    beta_ols[it] = beta[1]

    beta_q, _, _ = quantile_regression(y, design, alpha)
    beta_lad[it] = beta_q[1]

print("parameters      estimate      stderr")
print(f"OLS             {np.mean(beta_ols):.4f}        {np.std(beta_ols, ddof=1):.4f}")
print(f"LAD             {np.mean(beta_lad):.4f}        {np.std(beta_lad, ddof=1):.4f}")

OLS, LAD, quantile, and SVM regression compared — svm/svm6.py

"""Translated from Examples/svm/svm6.m -- comparison of OLS, LAD, and
SVM-regression (LS-SVM and epsilon-SVM, at two different C values)
estimates on the same 10-obs dataset used in ols1.m/robust1.m."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression
from quanttoolbox.svm.svm import svm_regression_primal

data = np.array(
    [
        [1.5, 1.0, 2.4, 3.6, 0.3],
        [20.4, 1.0, 1.1, 3.8, 5.9],
        [17.1, 1.0, 5.1, 6.3, 6.1],
        [30.9, 1.0, 2.7, 2.4, 9.5],
        [22.2, 1.0, 3.3, 3.0, 7.4],
        [9.1, 1.0, 1.0, 5.4, 4.9],
        [39.2, 1.0, 9.6, 2.8, 8.1],
        [3.1, 1.0, 2.9, 4.4, 1.0],
        [7.2, 1.0, 4.2, 5.6, 1.7],
        [27.6, 1.0, 8.1, 1.7, 5.4],
    ]
)
y = data[:, 0]
x_full = data[:, 1:5]

beta_ols = ols_estimation(y, x_full).beta
beta_lad = lad_regression(y, x_full).beta
beta_lad2, _, _ = quantile_regression(y, x_full, tau=0.5)

x = x_full[:, 1:4]  # drop the intercept and 4th column for SVM (3 predictors)

r_ls1 = svm_regression_primal(y, x, c=1)
r_eps1 = svm_regression_primal(y, x, c=1, epsilon=1)
r_ls2 = svm_regression_primal(y, x, c=100000)
r_eps2 = svm_regression_primal(y, x, c=100000, epsilon=0)

print("OLS:", np.round(beta_ols, 3))
print("LAD:", np.round(beta_lad, 3))
print("Quantile(0.5):", np.round(beta_lad2, 3))
print("SVM LS (C=1):", round(r_ls1.beta0, 3), np.round(r_ls1.beta, 3))
print("SVM eps (C=1, eps=1):", round(r_eps1.beta0, 3), np.round(r_eps1.beta, 3))
print("SVM LS (C=1e5):", round(r_ls2.beta0, 3), np.round(r_ls2.beta, 3))
print("SVM eps (C=1e5, eps=0):", round(r_eps2.beta0, 3), np.round(r_eps2.beta, 3))

OLS, median, LAD, and Huber regression compared — ects/robust2.py

"""Translated from Examples/ects/robust2.m -- OLS vs. median (quantile,
alpha=0.5) regression, LAD regression, and Huber regression, on a
simulated 3-predictor dataset.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import huber_regression, lad_regression

rng = np.random.default_rng(0)
n, k = 100, 3
beta_true = np.arange(1, k + 1, dtype=float)
sigma = 0.20

x = rng.random((n, k))
y = x @ beta_true + sigma * rng.standard_normal(n)

r_ols = ols_estimation(y, x)
u_ols = r_ols.residuals
print("OLS beta:", np.round(r_ols.beta, 4))

alpha = 0.50
beta_med, _, _ = quantile_regression(y, x, alpha)
print("\nMedian regression beta:", np.round(beta_med, 4))

r_lad = lad_regression(y, x)
print("\nLAD regression beta:", np.round(r_lad.beta, 4), "converged:", r_lad.converged)

c = np.quantile(np.abs(u_ols), 0.90)
r_huber = huber_regression(y, x, c=c)
print(
    f"\nHuber regression (c={c:.4f}) beta:",
    np.round(r_huber.beta, 4),
    "converged:",
    r_huber.converged,
)

OLS/median/LAD plus quantile regression via IRLS and exact LP — ects/robust3.py

"""Translated from Examples/ects/robust3.m -- OLS vs. median regression,
LAD regression, and quantile regression at alpha=0.90, both via the IRLS
M-estimator (`quantile_m_regression`) and the exact LP formulation
(`quantile_regression`), on a simulated 3-predictor dataset.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression, quantile_m_regression

rng = np.random.default_rng(0)
n, k = 100, 3
beta_true = np.arange(1, k + 1, dtype=float)
sigma = 0.20

x = rng.random((n, k))
y = x @ beta_true + sigma * rng.standard_normal(n)

r_ols = ols_estimation(y, x)
u_ols = r_ols.residuals
print("OLS beta:", np.round(r_ols.beta, 4))

alpha = 0.50
beta_med, _, _ = quantile_regression(y, x, alpha)
print("\nMedian regression beta:", np.round(beta_med, 4))

r_lad = lad_regression(y, x)
print("\nLAD regression beta:", np.round(r_lad.beta, 4), "converged:", r_lad.converged)

alpha = 0.90
r_qm = quantile_m_regression(y, x, alpha)
print(f"\nQuantile (alpha={alpha}) regression via IRLS M-estimator, beta:", np.round(r_qm.beta, 4))

beta_q, _, _ = quantile_regression(y, x, alpha)
print(f"Quantile (alpha={alpha}) regression via exact LP, beta:        ", np.round(beta_q, 4))

OLS/SVM-LS and quantile/SVM-epsilon regression on synthetic data — svm/svm8.py

"""Translated from Examples/svm/svm8.m -- OLS, LAD, quantile, and SVM
regression (primal, at two very different C values, LS- and
epsilon-insensitive) compared on a larger (n=1000) simulated dataset.

The original draws x/beta/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here for reproducibility
-- this is the fixed-seed substitution the tracker notes svm8.m as still
needing."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression
from quanttoolbox.svm.svm import svm_regression_primal

rng = np.random.default_rng(0)
n, k = 1000, 4
x = rng.random((n, k))
beta_true = 5 * rng.random(k)
u = 0.20 * rng.standard_normal(n)
beta0_true = -3.0
y = beta0_true + x @ beta_true + u

x_design = np.column_stack([np.ones(n), x])

beta_ols = ols_estimation(y, x_design).beta
beta_lad = lad_regression(y, x_design).beta
beta_lad2, _, _ = quantile_regression(y, x_design, tau=0.5)

# SVM regression (primal), C=1
r_ls1 = svm_regression_primal(y, x, c=1)  # LS-SVM
beta_svm_ls = np.concatenate([[r_ls1.beta0], r_ls1.beta])
r_eps1 = svm_regression_primal(y, x, c=1, epsilon=1)  # epsilon-SVM
beta_svm_epsilon = np.concatenate([[r_eps1.beta0], r_eps1.beta])

# SVM regression (primal), C=1000 (effectively unregularized)
r_ls2 = svm_regression_primal(y, x, c=1000)  # LS-SVM
beta_svm_ls2 = np.concatenate([[r_ls2.beta0], r_ls2.beta])
r_eps2 = svm_regression_primal(y, x, c=1000, epsilon=0)  # epsilon-SVM, epsilon=0
beta_svm_epsilon2 = np.concatenate([[r_eps2.beta0], r_eps2.beta])

results = np.column_stack(
    [beta_ols, beta_lad, beta_lad2, beta_svm_ls, beta_svm_epsilon, beta_svm_ls2, beta_svm_epsilon2]
)
print("Comparison of OLS, LAD, quantile(0.5), and SVM estimates")
print(
    "columns: OLS, LAD, Quantile(0.5), SVM-LS(C=1), SVM-eps(C=1,eps=1), SVM-LS(C=1000), SVM-eps(C=1000,eps=0)"
)
print(np.round(results, 3))

print("\nOLS vs. SVM-LS at large C (expect near-identical -- unregularized limit):")
print(np.round(np.column_stack([beta_ols, beta_svm_ls2]), 3))

print(
    "\nQuantile(0.5) vs. SVM-eps at large C, epsilon=0 (expect near-identical -- both approach LAD):"
)
print(np.round(np.column_stack([beta_lad2, beta_svm_epsilon2]), 3))

stats.regression.robust

Python alternatives

Hybrid: statsmodels.robust.robust_linear_model.RLM covers Huber, Tukey biweight, Andrew's wave, Hampel, and trimmed-mean M-estimators via IRLS — a superset of this module's Huber implementation. It does not cover LAD or quantile M-estimation directly (though QuantReg(q=0.5) is exactly LAD). Use RLM for general M-estimation; keep lad_regression/quantile_m_regression/inverse_quantile_m_regression for those specific losses.

quanttoolbox.stats.regression.robust

Robust M-estimator regression via Iteratively Reweighted Least Squares (IRLS).

Ported from QuantToolBox/ects/{robust_regression,robust_huber_regression, robust_lad_regression,robust_quantile_regression, robust_inverse_quantile_regression}.m

Translation notes:

• All four specific estimators (Huber, LAD, quantile, inverse-quantile) are thin wrappers supplying a (rho, rho_prime) loss-function pair to one shared IRLS core (robust_regression), exactly mirroring the original's structure.
• quantile_m_regression/inverse_quantile_m_regression here solve the quantile-loss problem via IRLS (a smooth M-estimation approximation), which is a different algorithm from quanttoolbox.stats.regression.quantile.quantile_regression (which solves the exact Koenker-Bassett LP formulation). IRLS is faster and matches the original MATLAB toolbox's approach here, but the LP version is the numerically exact one if precision matters more than speed.
• MATLAB's global ROBUST_eps convergence tolerance and global Print_Results are replaced by the RobustRegressionConfig dataclass (the original never actually set a value for ROBUST_eps in the files reviewed, so a conventional 1e-6 default is used here -- pass RobustRegressionConfig(eps=...) to match a specific original run).
• The original's rho_prime numerical-gradient fallback (numerical_gradient(rho, beta) when only rho is supplied) is reproduced with a simple central-difference approximation.
• MATLAB's cdftc/cdffc (upper-tail Student-t / F complementary CDFs) map directly to scipy.stats.t.sf / scipy.stats.f.sf.

RobustRegressionConfig(eps=1e-06, max_iters=500) dataclass

Replaces the ROBUST_eps global. See module docstring for the default's provenance.

huber_regression(y, x, c=1.345, config=None)

Huber M-estimator regression: quadratic loss for |residual| < c, linear loss beyond. c=1.345 is the conventional choice giving ~95% efficiency under normality.

Original: ects/robust_huber_regression.m

Source code in src/quanttoolbox/stats/regression/robust.py

def huber_regression(
    y: np.ndarray,
    x: np.ndarray,
    c: float = 1.345,
    config: RobustRegressionConfig | None = None,
) -> RobustRegressionResult:
    """Huber M-estimator regression: quadratic loss for |residual| < c,
    linear loss beyond. c=1.345 is the conventional choice giving ~95%
    efficiency under normality.

    Original: ects/robust_huber_regression.m
    """
    if config is None:
        config = RobustRegressionConfig()

    def rho(u: np.ndarray) -> np.ndarray:
        return (u**2) * (np.abs(u) < c) + c * np.abs(u) * (np.abs(u) >= c)

    def rho_prime(u: np.ndarray) -> np.ndarray:
        return 2 * u * (np.abs(u) < c) + c * np.sign(u) * (np.abs(u) >= c)

    return robust_regression(y, x, rho, rho_prime, config)

inverse_quantile_m_regression(y, x, alpha, config=None)

"Inverse" quantile regression (loss weighted by which side of zero the residual falls on, complementary to quantile_m_regression).

Original: ects/robust_inverse_quantile_regression.m

Source code in src/quanttoolbox/stats/regression/robust.py

def inverse_quantile_m_regression(
    y: np.ndarray,
    x: np.ndarray,
    alpha: float,
    config: RobustRegressionConfig | None = None,
) -> RobustRegressionResult:
    """ "Inverse" quantile regression (loss weighted by which side of zero
    the residual falls on, complementary to quantile_m_regression).

    Original: ects/robust_inverse_quantile_regression.m
    """
    if config is None:
        config = RobustRegressionConfig()

    def rho(u: np.ndarray) -> np.ndarray:
        return u * ((u > 0).astype(float) - alpha)

    def rho_prime(u: np.ndarray) -> np.ndarray:
        return (u > 0).astype(float) - alpha

    return robust_regression(y, x, rho, rho_prime, config)

lad_regression(y, x, config=None)

Least Absolute Deviations (median) regression.

Note: LAD's weight function 1/(u+eps) is not smooth near u=0, so the IRLS loop typically doesn't settle below a tight eps tolerance -- it will often run to max_iters even after beta has already stabilized at the optimum. Check that the coefficients look stable across the last few iterations rather than relying solely on result.converged for this estimator.

Original: ects/robust_lad_regression.m

Source code in src/quanttoolbox/stats/regression/robust.py

def lad_regression(
    y: np.ndarray, x: np.ndarray, config: RobustRegressionConfig | None = None
) -> RobustRegressionResult:
    """Least Absolute Deviations (median) regression.

    Note: LAD's weight function 1/(u+eps) is not smooth near u=0, so the
    IRLS loop typically doesn't settle below a tight ``eps`` tolerance --
    it will often run to ``max_iters`` even after beta has already
    stabilized at the optimum. Check that the coefficients look stable
    across the last few iterations rather than relying solely on
    ``result.converged`` for this estimator.

    Original: ects/robust_lad_regression.m
    """
    if config is None:
        config = RobustRegressionConfig()

    def rho(u: np.ndarray) -> np.ndarray:
        return np.abs(u)

    def rho_prime(u: np.ndarray) -> np.ndarray:
        return np.sign(u)

    return robust_regression(y, x, rho, rho_prime, config)

quantile_m_regression(y, x, alpha, config=None)

Quantile regression at level alpha via IRLS M-estimation (see module docstring for how this differs from the exact LP-based quantile.quantile_regression).

Original: ects/robust_quantile_regression.m

Source code in src/quanttoolbox/stats/regression/robust.py

def quantile_m_regression(
    y: np.ndarray,
    x: np.ndarray,
    alpha: float,
    config: RobustRegressionConfig | None = None,
) -> RobustRegressionResult:
    """Quantile regression at level alpha via IRLS M-estimation (see module
    docstring for how this differs from the exact LP-based
    ``quantile.quantile_regression``).

    Original: ects/robust_quantile_regression.m
    """
    if config is None:
        config = RobustRegressionConfig()

    def rho(u: np.ndarray) -> np.ndarray:
        return u * (alpha - (u < 0))

    def rho_prime(u: np.ndarray) -> np.ndarray:
        return alpha - (u < 0).astype(float)

    return robust_regression(y, x, rho, rho_prime, config)

robust_regression(y, x, rho, rho_prime=None, config=None)

Generic M-estimation regression via IRLS, for an arbitrary loss function rho (with derivative rho_prime, numerically approximated if not supplied).

Original: ects/robust_regression.m

Source code in src/quanttoolbox/stats/regression/robust.py

def robust_regression(
    y: np.ndarray,
    x: np.ndarray,
    rho: Callable[[np.ndarray], np.ndarray],
    rho_prime: Callable[[np.ndarray], np.ndarray] | None = None,
    config: RobustRegressionConfig | None = None,
) -> RobustRegressionResult:
    """Generic M-estimation regression via IRLS, for an arbitrary loss
    function rho (with derivative rho_prime, numerically approximated if
    not supplied).

    Original: ects/robust_regression.m
    """
    if config is None:
        config = RobustRegressionConfig()

    y = np.asarray(y, dtype=float).flatten()
    x = np.asarray(x, dtype=float)
    n_obs = y.shape[0]
    n_params = x.shape[1]

    if rho_prime is None:
        h = 1e-6

        def rho_prime(u: np.ndarray, _rho=rho, _h=h) -> np.ndarray:
            return (_rho(u + _h) - _rho(u - _h)) / (2 * _h)

    valid = ~np.isnan(y) & ~np.isnan(x).any(axis=1)
    valid_idx = np.where(valid)[0]
    n_missing = int(np.sum(~valid))
    y_v, x_v = y[valid], x[valid]
    n_valid = y_v.shape[0]
    n_valid_params = n_params

    beta = np.zeros(n_params)
    new_beta = np.linalg.solve(x_v.T @ x_v, x_v.T @ y_v)
    u = y_v - x_v @ new_beta
    w = rho_prime(u) / (u + 1e-10)
    xw = x_v * w[:, None]

    converged = False
    diff = np.max(np.abs(new_beta - beta))
    n_iter = 1
    while diff > config.eps:
        new_beta = np.linalg.solve(xw.T @ x_v, xw.T @ y_v)
        u = y_v - x_v @ new_beta
        w = rho_prime(u) / (u + 1e-5)
        xw = x_v * w[:, None]
        diff = np.max(np.abs(new_beta - beta))
        beta = new_beta
        n_iter += 1
        if n_iter > config.max_iters:
            break
    else:
        converged = True

    xx_w = xw.T @ x_v
    inv_xx_w = np.linalg.inv(xx_w)

    rss = float(np.sum(u**2))
    df_y = n_valid - 1
    df_x = n_valid_params - 1
    df_u = df_y - df_x

    sigma2 = rss / df_u
    sigma = np.sqrt(sigma2)
    vcv = sigma2 * inv_xx_w
    stderr = np.diag(vcv)
    stderr = np.where(stderr < 0, np.nan, stderr)
    stderr = np.sqrt(stderr)

    t_stat = beta / np.where(stderr == 0, np.nan, stderr)
    p_value = 2 * t_dist.sf(np.abs(t_stat), df_u)

    tss = float(np.sum(y_v**2))
    ess = tss - rss
    r_squared = 1 - rss / tss
    r_squared_adj = 1 - (rss / df_u) / (tss / df_y)

    yc = y_v - y_v.mean()
    tss_c = float(np.sum(yc**2))
    ess_c = tss_c - rss
    r_squared_c = 1 - rss / tss_c
    r_squared_c_adj = 1 - (rss / df_u) / (tss_c / df_y)

    f_stat = (r_squared_c / df_x) / ((1 - r_squared_c) / df_u)
    f_pvalue = f_dist.sf(f_stat, df_x, df_u)

    residuals = np.full(n_obs, np.nan)
    residuals[valid_idx] = u

    return RobustRegressionResult(
        beta=beta,
        stderr=stderr,
        vcv=vcv,
        residuals=residuals,
        converged=converged,
        n_iters=n_iter,
        weights=w,
        t_stat=t_stat,
        p_value=p_value,
        r_squared=r_squared,
        r_squared_adj=r_squared_adj,
        r_squared_centered=r_squared_c,
        r_squared_centered_adj=r_squared_c_adj,
        sigma2=sigma2,
        sigma=sigma,
        rss=rss,
        tss=tss,
        tss_centered=tss_c,
        ess=ess,
        ess_centered=ess_c,
        f_stat=f_stat,
        f_pvalue=f_pvalue,
        df_residual=df_u,
        n_obs=n_obs,
        n_obs_valid=n_valid,
        n_obs_missing=n_missing,
    )

Examples

OLS, LAD, quantile, and SVM regression compared — svm/svm6.py

"""Translated from Examples/svm/svm6.m -- comparison of OLS, LAD, and
SVM-regression (LS-SVM and epsilon-SVM, at two different C values)
estimates on the same 10-obs dataset used in ols1.m/robust1.m."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression
from quanttoolbox.svm.svm import svm_regression_primal

data = np.array(
    [
        [1.5, 1.0, 2.4, 3.6, 0.3],
        [20.4, 1.0, 1.1, 3.8, 5.9],
        [17.1, 1.0, 5.1, 6.3, 6.1],
        [30.9, 1.0, 2.7, 2.4, 9.5],
        [22.2, 1.0, 3.3, 3.0, 7.4],
        [9.1, 1.0, 1.0, 5.4, 4.9],
        [39.2, 1.0, 9.6, 2.8, 8.1],
        [3.1, 1.0, 2.9, 4.4, 1.0],
        [7.2, 1.0, 4.2, 5.6, 1.7],
        [27.6, 1.0, 8.1, 1.7, 5.4],
    ]
)
y = data[:, 0]
x_full = data[:, 1:5]

beta_ols = ols_estimation(y, x_full).beta
beta_lad = lad_regression(y, x_full).beta
beta_lad2, _, _ = quantile_regression(y, x_full, tau=0.5)

x = x_full[:, 1:4]  # drop the intercept and 4th column for SVM (3 predictors)

r_ls1 = svm_regression_primal(y, x, c=1)
r_eps1 = svm_regression_primal(y, x, c=1, epsilon=1)
r_ls2 = svm_regression_primal(y, x, c=100000)
r_eps2 = svm_regression_primal(y, x, c=100000, epsilon=0)

print("OLS:", np.round(beta_ols, 3))
print("LAD:", np.round(beta_lad, 3))
print("Quantile(0.5):", np.round(beta_lad2, 3))
print("SVM LS (C=1):", round(r_ls1.beta0, 3), np.round(r_ls1.beta, 3))
print("SVM eps (C=1, eps=1):", round(r_eps1.beta0, 3), np.round(r_eps1.beta, 3))
print("SVM LS (C=1e5):", round(r_ls2.beta0, 3), np.round(r_ls2.beta, 3))
print("SVM eps (C=1e5, eps=0):", round(r_eps2.beta0, 3), np.round(r_eps2.beta, 3))

OLS, median, LAD, and Huber regression compared — ects/robust2.py

"""Translated from Examples/ects/robust2.m -- OLS vs. median (quantile,
alpha=0.5) regression, LAD regression, and Huber regression, on a
simulated 3-predictor dataset.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import huber_regression, lad_regression

rng = np.random.default_rng(0)
n, k = 100, 3
beta_true = np.arange(1, k + 1, dtype=float)
sigma = 0.20

x = rng.random((n, k))
y = x @ beta_true + sigma * rng.standard_normal(n)

r_ols = ols_estimation(y, x)
u_ols = r_ols.residuals
print("OLS beta:", np.round(r_ols.beta, 4))

alpha = 0.50
beta_med, _, _ = quantile_regression(y, x, alpha)
print("\nMedian regression beta:", np.round(beta_med, 4))

r_lad = lad_regression(y, x)
print("\nLAD regression beta:", np.round(r_lad.beta, 4), "converged:", r_lad.converged)

c = np.quantile(np.abs(u_ols), 0.90)
r_huber = huber_regression(y, x, c=c)
print(
    f"\nHuber regression (c={c:.4f}) beta:",
    np.round(r_huber.beta, 4),
    "converged:",
    r_huber.converged,
)

OLS/median/LAD plus quantile regression via IRLS and exact LP — ects/robust3.py

"""Translated from Examples/ects/robust3.m -- OLS vs. median regression,
LAD regression, and quantile regression at alpha=0.90, both via the IRLS
M-estimator (`quantile_m_regression`) and the exact LP formulation
(`quantile_regression`), on a simulated 3-predictor dataset.

The original draws x/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression, quantile_m_regression

rng = np.random.default_rng(0)
n, k = 100, 3
beta_true = np.arange(1, k + 1, dtype=float)
sigma = 0.20

x = rng.random((n, k))
y = x @ beta_true + sigma * rng.standard_normal(n)

r_ols = ols_estimation(y, x)
u_ols = r_ols.residuals
print("OLS beta:", np.round(r_ols.beta, 4))

alpha = 0.50
beta_med, _, _ = quantile_regression(y, x, alpha)
print("\nMedian regression beta:", np.round(beta_med, 4))

r_lad = lad_regression(y, x)
print("\nLAD regression beta:", np.round(r_lad.beta, 4), "converged:", r_lad.converged)

alpha = 0.90
r_qm = quantile_m_regression(y, x, alpha)
print(f"\nQuantile (alpha={alpha}) regression via IRLS M-estimator, beta:", np.round(r_qm.beta, 4))

beta_q, _, _ = quantile_regression(y, x, alpha)
print(f"Quantile (alpha={alpha}) regression via exact LP, beta:        ", np.round(beta_q, 4))

OLS/SVM-LS and quantile/SVM-epsilon regression on synthetic data — svm/svm8.py

"""Translated from Examples/svm/svm8.m -- OLS, LAD, quantile, and SVM
regression (primal, at two very different C values, LS- and
epsilon-insensitive) compared on a larger (n=1000) simulated dataset.

The original draws x/beta/u from MATLAB's unseeded `rand`/`randn`; a fixed
seed (`np.random.default_rng(0)`) is substituted here for reproducibility
-- this is the fixed-seed substitution the tracker notes svm8.m as still
needing."""

import numpy as np

from quanttoolbox.econometrics.estimation import ols_estimation
from quanttoolbox.stats.regression.quantile import quantile_regression
from quanttoolbox.stats.regression.robust import lad_regression
from quanttoolbox.svm.svm import svm_regression_primal

rng = np.random.default_rng(0)
n, k = 1000, 4
x = rng.random((n, k))
beta_true = 5 * rng.random(k)
u = 0.20 * rng.standard_normal(n)
beta0_true = -3.0
y = beta0_true + x @ beta_true + u

x_design = np.column_stack([np.ones(n), x])

beta_ols = ols_estimation(y, x_design).beta
beta_lad = lad_regression(y, x_design).beta
beta_lad2, _, _ = quantile_regression(y, x_design, tau=0.5)

# SVM regression (primal), C=1
r_ls1 = svm_regression_primal(y, x, c=1)  # LS-SVM
beta_svm_ls = np.concatenate([[r_ls1.beta0], r_ls1.beta])
r_eps1 = svm_regression_primal(y, x, c=1, epsilon=1)  # epsilon-SVM
beta_svm_epsilon = np.concatenate([[r_eps1.beta0], r_eps1.beta])

# SVM regression (primal), C=1000 (effectively unregularized)
r_ls2 = svm_regression_primal(y, x, c=1000)  # LS-SVM
beta_svm_ls2 = np.concatenate([[r_ls2.beta0], r_ls2.beta])
r_eps2 = svm_regression_primal(y, x, c=1000, epsilon=0)  # epsilon-SVM, epsilon=0
beta_svm_epsilon2 = np.concatenate([[r_eps2.beta0], r_eps2.beta])

results = np.column_stack(
    [beta_ols, beta_lad, beta_lad2, beta_svm_ls, beta_svm_epsilon, beta_svm_ls2, beta_svm_epsilon2]
)
print("Comparison of OLS, LAD, quantile(0.5), and SVM estimates")
print(
    "columns: OLS, LAD, Quantile(0.5), SVM-LS(C=1), SVM-eps(C=1,eps=1), SVM-LS(C=1000), SVM-eps(C=1000,eps=0)"
)
print(np.round(results, 3))

print("\nOLS vs. SVM-LS at large C (expect near-identical -- unregularized limit):")
print(np.round(np.column_stack([beta_ols, beta_svm_ls2]), 3))

print(
    "\nQuantile(0.5) vs. SVM-eps at large C, epsilon=0 (expect near-identical -- both approach LAD):"
)
print(np.round(np.column_stack([beta_lad2, beta_svm_epsilon2]), 3))