Logit Adjustment for Modeling Probabilities Using Scipy

logit_adjustment.py

from scipy.special import loggamma
from scipy.special import expit, logit
import numpy as np

from scipy.optimize import minimize

def logLikelihood(params, y, X):
    b = np.array(params[0:-1])      # the beta parameters of the regression model
    phi = params[-1]                # the phi parameter
    mu = expit(np.dot(X,b))
   
    eps = 1e-6                      # used for safety of the gamma and log functions avoiding inf
    res = - np.sum(loggamma(phi+eps) # the log likelihood
                   - loggamma(mu*phi+eps) 
                   - loggamma((1-mu)*phi+eps) 
                   + (mu*phi-1)*np.log(y+eps) 
                   + ((1-mu)*phi-1)*np.log(1-y+eps))

    return res


class Logit_Adj():

    def fit(self, X, y):
        # initial parameters for optimization

        n = X.shape[1]

        bounds = [(None, None)] * n + [(0,None)]

        phi = 1
        b0 = 1
        x0 = np.array([b0] * n + [phi])

        self.res = minimize(logLikelihood, x0=x0, args=(y,X), bounds=bounds)
        
    def predict(self, X):
        b = np.array(res.x[0:X.shape[1]])   # optimal regression parameters
        return expit(np.dot(X, b))