qra/python/logistic/logit_utils_gen.py

"""
    The file provides a class for logistic regression utilities with polynomial 
    logit (log odds) function:

        p(x|pars) = 1/(1 + exp(-F(x|beta(pars))))
    
    with log odds F of polynomial form:

        F(x|beta) = sum_{j=0}^degree beta_j x^j
                
    with decision function (aka logit) F and parameters 
        
        beta(pars) = [beta(pars)_j]_{j=0}^degree
    
    where pars are regression parameters and len(pars) = degree + 1.
    
    Coefficients beta(pars) can be constrained to be monotonic 
    function of x by using monotonic cubic transformation:
        
        beta(pars) = mc.forward_map(pars)
    
    where mc is module mono_cubic2.

    NOTES: The model defines the conditional probability 

        Prob(Y = y|x, pars) = 1/(1 + exp(-s(y) F(x| beta(pars))))
    
    where
    
        s(y) = 2*y - 1
    
    with x in R and y in {0,1}. For degree = 1 this is standard logistic regression 
    
        beta(pars) = (pars[0], pars[1]). 
    
    and generally without monotonicity condition

        beta(pars) = [pars_i]_{i=0}^degree

    For degree = 3 and mono = True the coefficients beta(pars) are constrained to be 
    monotonic by using monotonic cubic transformation.

    We have data

        {(x_i, y_i) : x_i in R, y_i in {0,1}, i = 0, ..., n-1 }
    
    and the model is fit to data by minimizing negative log-likelihood function:

        nlff = -sum_i log(Prob(Y = y_i| x_i, pars))  : neg. log likelihood

    with respect to parameters pars. We can have regularization term in cost function
    and this case we minimize cost function:

        cost(pars) = nllf(pars) + lambda ||pars||^2

    Author: Martin Horvat, January 2026
"""


import numpy as np
import scipy
import scipy.optimize

import mono_cubic2 as mc

def resize_with_const(v, n, val=0):
    """
    Resize a 1D vector to a specified length `n`.

    If the input vector `v` is longer than `n`, it is truncated.
    If it is shorter, it is padded with the constant value `val`.
    If it is already of length `n`, it is returned unchanged.

    Parameters:
        v (array-like): Input 1D vector (list or NumPy array).
        n (int): Target length of the output vector.
        val (scalar, optional): Value used to pad if `v` is shorter than `n`. Default is 0.

    Returns:
        np.ndarray: Resized 1D NumPy array of length `n`.
    """
    v = np.asarray(v)
    
    if len(v) == n: return v
    
    if len(v) > n: return v[:n]
    
    return np.concatenate([v, np.full(n - len(v), val)])

def safe_exp(x, max_exp = 700):
    """
    A numerically robust version of np.exp that avoids overflow by clipping the input.
    
    Parameters:
        x : array_like
            Input value or array.
        max_exp : float
            Maximum allowed exponent value. np.exp(709) ≈ 8.2e307 (close to float64 max).
    
    Returns:
        array_like
            The exponential of the input with overflow protection.
    """

    return np.exp(np.clip(x, -max_exp, max_exp))

def safe_expit(x, max_exp = 700): return 1/(1 + safe_exp(-x, max_exp))

"""
    Fitting data

        {(x_i, y_i) : x_i in R, y_i in {0,1}, i = 0, ..., n-1 }
    
    to model function 

        f(x|pars) = 1/(1 + exp(-F(x|beta(pars))))

    with log odds of polynomial form:

        log(f(x|pars)/(1 - f(x|pars))) = F(x| beta(pars)) 

    where F is decision function (aka logit)
        
        F(x|beta) = sum_{i=0}^degree beta_i x^i
                
    and coefficients
        
        beta(pars) = [beta_i(pars)]_{i=0}^degrees
    
    Conditional probability

        Prob(Y = y_i|x, pars) = 1/(1 + exp(-s(y) F(x|beta(pars))) 

    with

        s(y) = 2*y -1
"""
class LogisticPolyRegression:

    """
        Class constructor. 

        Input:
            degree: int, degree of polynomial
            mono: boolean, default False
            lambda: None or tuple float, L2 regularization

    """
    def __init__(self, degree = 1, mono = False, lam = None):

        self.degree = degree
        self.mono = mono
        self.big = 200
        self.small = 1e-8
        self.lam = lam

        if mono:
            assert self.degree in [1, 3], f"Degree {self.degree} not supported in mono!"

        self.mono3 = self.mono and (self.degree == 3)

    """
        Mapping regression parameters pars to coefficients beta 
            
            beta = beta(pars)
        
        used in decision function:

            F(x|beta) = sum_{i=0}^degree beta_i x^i

        Input:
            pars
        
        Return:
            beta
    """
    def get_beta(self, pars):
        
        beta = mc.forward_map(pars) if self.mono3 else pars
        return np.array(beta)
    
    """
        Mapping beta to regression parameters used in decision function.

        Input:
            beta: coefficient beta
        
        Return:
            pars
    """
    def get_pars(self, beta):

        pars = mc.backward_map(beta) if self.mono3 else beta
        return np.array(pars)

    """
        Calculate jacobian between beta and  regression parameters
            
            J = d(beta)/d(pars)
              = [d(beta_i)/d(pars_a)]_{i,a}
        and
            
            H = [d^2 beta_i/(d(pars_a) d(pars_b))]_{i,a,b}

        Input:
            pars
            hess: boolean, False
        
        Return:
            J       if hess = True
            (J, H)  if hess = False
    """
    def get_jac_beta(self, pars, hess = False):
        
        n = len(pars)
            
        J = mc.forward_map_jacobian(pars) if self.mono3 else np.eye(n)
            
        if not hess: return J
        
        H = mc.forward_map_hessian(pars) if self.mono3 else np.zeros(shape = (n, n, n))

        return (J, H)


    """
        Model function 

            f(x|pars) = 1/(1 + exp(-F(x|beta)))  beta = beta(pars)

        Input:
            x: scalar value or a array of values
            pars: array of r = degree + 1 floats, model parameters array of floats

        Return:
            model function values 
    """
    def model(self, x, pars):

        beta = self.get_beta(pars)  # beta
        X = np.column_stack([x**i for i in range(len(beta))])
        F = X @ beta              # decision function, X beta

        return safe_expit(F)

    """
        Calculate negative log-likelihood function

            nllf = -sum_i log(Prob(Y = y_i| x_i, pars))   : neg. log likelihood
            
            grad = [d(nllf)/d(pars_a)]_a                  : jacobian
        
        where

            Prob(Y = y_i| x_i, pars) = 1/(1 + exp(-s_i F(x_i| beta)))  beta=beta(pars)
            s_i = 2*y_i  - 1
     
        Input:
            x: array of n floats
            y: array of n int in {0,1}  
            pars: array of r = degree + 1 floats, model parameters array of floats
            jac: boolean, default False, if jacobian is needed
        
        Return:
            nllf           : if jac is false
            (nllf, grad)   : if jac is true
    """
    def nllf(self, x, y, pars, jac = False):

        beta = self.get_beta(pars)
        X = np.column_stack([x**i for i in range(len(beta))])
    
        # signs
        s = 2.0*y - 1
        
        # decision function for conditional probability Prob(Y = y| x)
        F = s*(X @ beta)

        nllf = np.sum(np.log(1 + safe_exp(-F))) 

        if not jac: return nllf
        
        J = self.get_jac_beta(pars)

        grad = -(s*safe_expit(-F)) @ (X @ J)

        return (nllf, grad)

    """
        Penalty function

         Input:
            pars: array of r = degree + 1 floats, model parameters array of floats
            jac: boolean, default False, if jacobian is needed
        
        Return:
            val           : if jac is false
            (val, grad)   : if jac is true
    """
    def penalty(self, pars, jac = False):
        
        val = self.lam[0]*np.sum(np.abs(pars)) + self.lam[1]*np.sum(pars**2)

        if jac:
            grad = self.lam[0]*np.sign(pars) + 2*self.lam[1]*pars
            return (val, grad)
           
        return val

    """
        Cost function
    """
    def cost(self, x, y, pars, jac = False):

        val = self.nllf(x, y, pars, jac)
        
        if self.lam is not None:
            pen = self.penalty(pars, jac)
            return (val[0] + pen[0], val[1] + pen[1]) if jac else val + pen
        
        return val

    """
        Estimate parameters.

        Input:
            x: array of n floats
            y: array of n int in {0,1}
        
        Return:
            pars0
    """
    def get_est_pars(self, x, y):

        L = np.log(2*len(x) + 1)

        if self.mono3:    
            z = L*np.polyfit(x, 2.0*y - 1, 1)[::-1]
            beta = resize_with_const(z, self.degree + 1, 1e-8)    
        else:
            beta = L*np.polyfit(x, 2.0*y - 1, self.degree)[::-1]

        return self.get_pars(beta)
    
    """
        Performing logistic regression with log odds of polynomial form:

            log(f(x|pars)/(1 - f(x|pars))) = F(x|beta)  beta = beta(pars)
        
        and this gives           
            
            f(x|pars) = 1/(1 + exp(-F(x|beta)))

        where coefficient beta = [beta_i(pars)]_{i=0}^degree, with decision 
        function (aka logit)

            F(x|beta) = sum_{i=0}^degree beta_i x^i
                
        Input:
            x: array of n floats
            y: array of n int in {0, 1}
        
        Return:
            pars
    """
    def fit(self, x, y, pars0 = None, method = "local"):
        
        bnds = [(-self.big, self.big)]*(self.degree + 1)
  
        if method == "local":

            pars0 = self.get_est_pars(x, y)
            cf = lambda pars: self.cost(x, y, pars, jac = True)
            res = scipy.optimize.minimize(cf, x0 = pars0, method = 'L-BFGS-B', 
                                          jac = True, bounds = bnds, tol=1e-12)

        elif method == "diff_evol":
            
            cf = lambda pars: self.cost(x, y, pars, jac = False)
            res = scipy.optimize.differential_evolution(cf, bounds = bnds, 
                                                        tol = 1e-8, polish = False)

            cf = lambda pars: self.cost(x, y, pars, jac = True)
            res = scipy.optimize.minimize(cf, x0 = res.x, method = 'L-BFGS-B', 
                                          jac = True, bounds = bnds, tol=1e-12)

        elif method == "anneal":

            cf = lambda pars: self.cost(x, y, pars, jac = False)
            res = scipy.optimize.dual_annealing(cf, bounds = bnds)
        
        else:
            assert False, "This method is not supported."
        
        return {"pars": res.x, "cost": res.fun, "success": res.success}

    """
        Producing goodness of fit measures:
        
            LLF = log_likelihood function
            AIC = Akaike information criterion
            BIC = Bayesian information criterion

        Input:
            x: array of n floats
            y: array of n int in {0,1}
            pars: array of r = degree+1 floats, model parameters
            thresh: float, default 0.5, threshold value for classification

        Return:
            {"n": n, 
            "k": r, 
            "dof": n - r, 
            "LLF": log_likelihood, 
            "AIC": AIC, 
            "BIC": BIC, 
            "A": classification accuracy (threshold values = 0.5 prob)}
        
        Ref:
            https://en.wikipedia.org/wiki/Logistic_regression
            https://en.wikipedia.org/wiki/Akaike_information_criterion
            https://www.medicine.mcgill.ca/epidemiology/joseph/courses/epib-621/logfit.pdf
    """
    def goodness_of_fit(self, x, y, pars, thresh = 0.5):
        
        # model probabilities
        p = self.model(x, pars)

        # log likelihood
        llf = -self.nllf(x, y, pars)
        
        # information criteria
        k, n = len(pars), len(x)
        AIC = 2*k - 2*llf
        BIC = k*np.log(n) - 2*llf

        # chi2
        dof = n - k
        r = (y - p)/np.sqrt(p*(1-p) + self.small)
        chi2 = np.sum(r**2)
        p_val = scipy.stats.chi2.sf(chi2, dof)

        # using model as classifier
        matches = y == np.heaviside(p - thresh, 1)

        # accuracy A
        A = np.count_nonzero(matches)/n

        return {"LLF": llf, 
                "AIC": AIC, 
                "BIC": BIC, 
                "A" : A, 
                "chi2": chi2, 
                "p-value(chi2)": p_val,  # not very useful
                "n": n, "k": k, "dof": dof}

    """
        Calculation of asymptotic variance-covariance matrix of regression 
        parameters pars

           cov_{asymp}[pars] = H^{-1}

        where H is hessian of nllf 

           H = [d^2(nllf)/(d(pars)_a d(pars)_b ]_{a,b}
        
        for the logistic regression of the polynomial model:
        
            log(f(x)/(1 - f(x))) ~ sum_{i=0}^degree b_i(pars) x^i

        Input:
            x: array of n floats
            pars: array of r = degree+1 floats, model parameters
            
        Return:
            array of rxr floats; r = degree + 1
        
        Ref:
          * https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients
          * https://goodboychan.github.io/machine_learning/2020/09/14/02-Regularized-likelihood-methods.html
    """
    def cov(self, x, y, pars):
        
        # coefficients
        beta = self.get_beta(pars)
        
        # design matrix -- add column of 1's at the beginning of your X_train matrix
        X = np.column_stack([x**i for i in range(len(beta))])

        # Jacobian J = [dbeta_i/dpars_j]_{ij}
        J, H = self.get_jac_beta(pars, hess = True)

        # signs
        s = 2.0*y - 1
        
        # decision function for conditional probability Prob(Y = y| x)
        F = s*(X @ beta)

        # probabilities p_i = P(Y=y_i | x_i)
        p = safe_expit(F)
        q = 1 - p

        # calculate hessian
        L = X @ J
        H = (L.T*(q*p))@L - np.tensordot((s*q)@X, H, axes = ([0], [0]))

        if self.lam is not None: 
            Hp = H + 2*self.lam[1]*np.eye(len(pars))   # H' = H + lambda id 
            iHp = np.linalg.inv(Hp)                    # inv(H')

            return iHp@H@iHp
        
        # covariance matrix C_params = H^-1
        return np.linalg.inv(H)

    """
        Calculating quantiles of the model parameters at given probabilities p 
        for normal distribution of parameters:

            pars ~ N(mean_pars, cov_pars)
        
        Input:
            probs: array of m floats, probabilities
            mean_pars: array of r = degree+1 floats, mean model parameters
            cov_pars: array of rxr floats, variance-covariance matrix of parameters
        
        Return:
            array of mxr floats
    """
    def get_pars_quantiles_normal(self, probs, mean_pars, cov_pars):
        
        # mean and standard variance parameters
        locs = mean_pars
        scales = np.sqrt(np.diag(cov_pars))

        # computing quantiles of parameters
        return locs + np.outer(scipy.stats.norm.ppf(probs), scales)

    """
        Calculating quantiles of the model values 
        
            f(x|pars) = 1/(1 + exp(-F(x|beta)))     beta = beta(pars)
        
        with

            F(x|beta) = sum_{i=0}^degree x^i beta_i
            
        at given probabilities p and values x assuming 
        normal distribution of parameters:

            pars ~ N(mean_pars, cov_pars)
        
        This distribution is asymptotic MLE distribution of parameters.

        Input:
            x: array of n float
            probs: array of m floats, probabilities
            mean_pars: array of r = degree+1 floats, mean model parameters
            cov_pars: array of rxr floats, variance-covariance matrix of parameters
        
        Return:
            array of mxn floats
    """

    def get_model_quantiles_normal(self, x, probs, mean_pars, cov_pars, 
                                   exact = True, seed = 1977, m = 10**5):
        
        mean_beta = self.get_beta(mean_pars)
        X = np.column_stack([x**i for i in range(len(mean_beta))])
        
        if exact and self.mono3:

            # init random generator
            rng = np.random.default_rng(seed)
            pars = rng.multivariate_normal(mean_pars, cov_pars, size = m)

            # get betas
            beta = np.apply_along_axis(self.get_beta, 1, pars)
            
            # quantiles of decision function
            Q = np.quantile(X@beta.T, probs, axis = 1)
            
            return np.apply_along_axis(safe_expit, 1, Q)

        # J = d(beta)/d(pars)
        J = self.get_jac_beta(mean_pars)
        
        # transform data
        S = X@J

        # mean and standard variance of logit (aka log of odds)
        locs = X@mean_beta
        scales = np.sqrt(np.diag(S@cov_pars@S.T))

        # computing quantiles of logit
        Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)
        
        # convert logit to expit
        return safe_expit(Q)

    """
        Calculating quantiles using delta method of the model values 
        
            f(x|pars) = 1/(1 + exp(-F(x|beta))) beta = beta(pars)

        with

            F(x|beta) = sum_{i=0}^degree x^i beta_i

        at given probabilities p and values x assuming 
        normal distribution of parameters :

            pars ~ N(mean_pars, cov_pars)

        This distribution is asymptotic MLE distribution of parameters.
        We approximate exact model with linear expansion

            f(x|pars) = p(x|pars_mean) + dp/db(x|pars_mean) (pars - pars_mean)
        
        and the last term is normally distributed.

        Input:
            x: array of n float
            probs: array of m floats, probabilities
            mean_pars: array of r = degree+1 floats, mean model parameters
            cov_pars: array of rxr floats, variance-covariance matrix of parameters 
        
        Return:
            array of mxn floats
    """
    def get_model_quantiles_delta(self, x, probs, mean_pars, cov_pars):
        
        mean_beta = self.get_beta(mean_pars)
        X = np.column_stack([x**i for i in range(len(mean_beta))])
        F = X@mean_beta

        # J = d(beta)/d(pars)
        J = self.get_jac_beta(mean_pars)

        # S = d(F)/d(pars)
        S = X@J

        # attributes of normal distribution of model values
        locs = safe_expit(F)
        scales = np.sqrt(np.diag(S@cov_pars@S.T))/(4*np.cosh(F/2)**2)

        # computing quantiles of logit
        Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)

        return np.clip(Q, a_min = 0, a_max = 1)

    """
        Generate parameters assuming normal distribution.

        Input:
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator
    
        Return:
            array of mx(degree + 1)
    """
    def get_normal_pars(self, x, y, m, seed = 1977):
        
        res = self.fit(x, y, method = "diff_evol")

        assert res["success"], "Fit did not succeed."

        # optimal parameters
        pars = res["pars"]

        # covariance matrix of parameters
        cov = self.cov(x, y, pars)

        # init random generator
        rng = np.random.default_rng(seed)

        return rng.multivariate_normal(pars, cov, size = m)

    """
        Generate m parameters via non-parametric bootstrapping with a minimal constraint 
        that both groups should be present in the sampled data:

        boostrapped sample = (xb, yb)  by sampling with replacement pairs (x_i, y_i) 
                                    with condition that yb can not be just 0 or just 1 

        Input:
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator
    
        Return:
            array of mx(degree + 1)
    """
    def get_nonparam_boots_pars(self, x, y, m, seed = 1977):
        
        # fitting original data
        res_fit = self.fit(x, y, method = "diff_evol")
        assert res_fit["success"]

        pars0 = res_fit["pars"]

        n = len(x)

        rng = np.random.default_rng(seed)

        # generate parameters
        lst = [pars0]
        while True:

            # create set indices for sampling with replacement + constraint
            idx = rng.choice(n, n)
            if np.sum(y[idx]) in [0, n]: continue

            # do fitting 
            res_fit = self.fit(x[idx], y[idx], pars0, method="local")
            if not res_fit["success"]: continue
            
            # store pars
            lst.append(res_fit["pars"])
            if len(lst) == m: break

        return np.array(lst)


    """
        Generate m parameters via non-parametric stratified bootstrapping:

            boostrapped sample = (xb, yb)  
                
                xb = (sampled with replacement from x0, sampled with replacement from x1)
                yb = (0...0, 1...1)

        Input: 
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator
        
        Return:
            array of mx(degree + 1)
    """
    def get_nonparam_stratified_boots_pars(self, x, y, m, seed = 1977):

        # pars of original data
        res_fit = self.fit(x, y, method = "diff_evol")
        assert res_fit["success"]

        pars0 = res_fit["pars"]
        
        # statistics about groups
        xs = [x[y == i] for i in range(2)]
        ns = [len(e) for e in xs]

        # common vector states
        yb = np.concatenate([np.full(ns[i], i) for i in range(2)])

        rng = np.random.default_rng(seed)

        # generate parameters
        lst = [pars0]
        for _ in range(m):
            # stratified sampling with replacement
            xb = np.concatenate([rng.choice(xs[i], ns[i]) for i in range(2)])

            # do fitting 
            res_fit = self.fit(xb, yb, pars0, method="local")
            if not res_fit["success"]: continue
            
            # store pars
            lst.append(res_fit["pars"])
            if len(lst) == m: break

        return np.array(lst)

    """
        Generate m parameters via parametric bootstrapping:

            boostrapped sample = (x, yb) yb ~ B(model(x, fitted pars))  
        
        where B is Bernoulli distribution

        Input: 
            x : array of n floats
            y : array of n int in {0,1}
            m: integer, number of samples
            seed : int, seed for the random generator
        
        Return:
            array of mx(degree + 1)

        Ref:

        * https://www.scirp.org/journal/paperinformation?paperid=70962
        * https://en.wikipedia.org/wiki/Bernoulli_distribution
    """
    def get_parametric_boots_pars(self, x, y, m, seed = 1977):

        # first discuss original dataset
        res_fit = self.fit(x, y, method = "diff_evol")
        assert res_fit["success"]

        pars0 = res_fit["pars"]

        p = self.model(x, pars0)
    
        n = len(x)
        
        rng = np.random.default_rng(seed)

        # generate parameters
        lst = [pars0]
        while True:

            # Generate new binary outcomes from Bernoulli(p_i)
            y_sim = rng.binomial(n = 1, p = p)
            if np.sum(y_sim) in [0, n]: continue
            
            # do fitting 
            res_fit = self.fit(x, y_sim, pars0, method="local")
            if not res_fit["success"]: continue
            
            # store pars
            lst.append(res_fit["pars"])
            if len(lst) == m: break

        return np.array(lst)