|
@@ -0,0 +1,443 @@
|
|
|
|
|
+import numpy as np
|
|
|
|
|
+import scipy
|
|
|
|
|
+
|
|
|
|
|
+"""
|
|
|
|
|
+ Bayesian model describing conditional probability
|
|
|
|
|
+
|
|
|
|
|
+ Prob(X|Y = 1)
|
|
|
|
|
+ = Prob(Y=1)p(X|Y=1)/(Prob(Y=0) p(X|Y=0) + Prob(Y=1) p(X|Y=1))
|
|
|
|
|
+ = 1 /(1 + O p(X|Y=0)/p(X|Y=1))
|
|
|
|
|
+ = 1/(1 + exp( -F))
|
|
|
|
|
+
|
|
|
|
|
+ where F is the decision function
|
|
|
|
|
+
|
|
|
|
|
+ F = log(p(X|Y=1)) - log(p(X|Y=0)) - log(O)
|
|
|
|
|
+
|
|
|
|
|
+ and O are the odds
|
|
|
|
|
+
|
|
|
|
|
+ O = Prob(Y=0)/Prob(Y=1)
|
|
|
|
|
+"""
|
|
|
|
|
+
|
|
|
|
|
+class BayesianModelRegression:
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Constructor
|
|
|
|
|
+
|
|
|
|
|
+ Input:
|
|
|
|
|
+ odds: float, ratio Prob(Y=0)/Prob(Y=1)
|
|
|
|
|
+ log_pdfs: function (x, pars, choice = 2)
|
|
|
|
|
+ match choice:
|
|
|
|
|
+ case 0: return log_pdf0
|
|
|
|
|
+ case 1: return log_pdf1
|
|
|
|
|
+ case _: return (log_pdf0, log_pdf1)
|
|
|
|
|
+ bounds: tuple of bounds, (bounds0, bound1)
|
|
|
|
|
+
|
|
|
|
|
+ distr_sample: function (rng, x, pars, n, choice):
|
|
|
|
|
+ generate n sampled os points using pdfs(choice, pars)
|
|
|
|
|
+ """
|
|
|
|
|
+ def __init__(self, odds, log_pdfs, bounds, distr_sample = None):
|
|
|
|
|
+
|
|
|
|
|
+ self.odds = odds
|
|
|
|
|
+ self.log_odds = np.log(odds)
|
|
|
|
|
+ self.log_pdfs = log_pdfs
|
|
|
|
|
+ self.bounds = bounds
|
|
|
|
|
+ self.distr_sample = distr_sample
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Calculate decision function:
|
|
|
|
|
+
|
|
|
|
|
+ decision = log(p(X|Y=1)) - log(p(X|Y=0)) - log(odds)
|
|
|
|
|
+
|
|
|
|
|
+ Input:
|
|
|
|
|
+ x: float or array of floats
|
|
|
|
|
+ pars: parameters for log_pdfs
|
|
|
|
|
+ Return:
|
|
|
|
|
+ float or array of float
|
|
|
|
|
+ """
|
|
|
|
|
+ def decision(self, x, pars):
|
|
|
|
|
+
|
|
|
|
|
+ # log of pdf for each group
|
|
|
|
|
+ lf0, lf1 = self.log_pdfs(x, pars)
|
|
|
|
|
+
|
|
|
|
|
+ return lf1 - lf0 - self.log_odds
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Calculate model of the conditional probability Prob(X|Y = 1)
|
|
|
|
|
+ Input:
|
|
|
|
|
+ x: float or array of floats
|
|
|
|
|
+ pars: parameters for log_pdfs
|
|
|
|
|
+ """
|
|
|
|
|
+ def model(self, x, pars):
|
|
|
|
|
+
|
|
|
|
|
+ # decision function
|
|
|
|
|
+ F = self.decision(x, pars)
|
|
|
|
|
+
|
|
|
|
|
+ # calculating model
|
|
|
|
|
+ return 1/(1 + np.exp(-F))
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Negative Log Likelihood function:
|
|
|
|
|
+
|
|
|
|
|
+ neg. log likelihood = -sum_i log(Prob(X = x_i, Y = y_i))
|
|
|
|
|
+
|
|
|
|
|
+ where
|
|
|
|
|
+
|
|
|
|
|
+ Prob(X, Y = 1) = 1/(1 + exp(-F))
|
|
|
|
|
+ Prob(X, Y = 0) = 1 - Prob(X, Y = 1) = 1/(1 + exp(+F))
|
|
|
|
|
+
|
|
|
|
|
+ with
|
|
|
|
|
+
|
|
|
|
|
+ log Prob(X, Y = y) = -log(1 + exp(-S(y) F))
|
|
|
|
|
+ S(y) = [ +1 : y = 1
|
|
|
|
|
+ [ -1 : y = 0
|
|
|
|
|
+
|
|
|
|
|
+ Input:
|
|
|
|
|
+ x: array of floats
|
|
|
|
|
+ y: array of ints in {0,1}
|
|
|
|
|
+ pars: array of floats, model parameters
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ float: negative log likelihood
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ def nllf(self, x, y, pars):
|
|
|
|
|
+
|
|
|
|
|
+ # signs for the groups:
|
|
|
|
|
+ # group 1 has + sign and group 0 has - sign
|
|
|
|
|
+ S = 2.0*y - 1
|
|
|
|
|
+
|
|
|
|
|
+ # decision function
|
|
|
|
|
+ F = self.decision(x, pars)
|
|
|
|
|
+
|
|
|
|
|
+ return np.sum(np.log(1 + np.exp(-S*F)))
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ MLE fitting of a distribution, given by log_pdf, to data x associated
|
|
|
|
|
+ to the group 0 or 1 by maximizing
|
|
|
|
|
+
|
|
|
|
|
+ loglikehood_{single group} = sum_i log_pdf(x | pars)
|
|
|
|
|
+
|
|
|
|
|
+ Input:
|
|
|
|
|
+ x: array of floats
|
|
|
|
|
+ choice: int in {0,1}, selecting the group
|
|
|
|
|
+ method: string in ["local", "diff_evol", "anneal"]
|
|
|
|
|
+ seed: int, seed of the random generator
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ {"pars": pars_MLE, "cost": NLLF at pars_MLE}
|
|
|
|
|
+ """
|
|
|
|
|
+
|
|
|
|
|
+ def fit_distr(self, x, choice, method = "local", seed = 1977):
|
|
|
|
|
+ fname = "fit_distr"
|
|
|
|
|
+
|
|
|
|
|
+ cost = lambda pars: -np.sum(self.log_pdfs(x, pars, choice))
|
|
|
|
|
+ bnds = self.bounds[choice]
|
|
|
|
|
+
|
|
|
|
|
+ match method:
|
|
|
|
|
+ case "local":
|
|
|
|
|
+ # random parameters from boundaries
|
|
|
|
|
+ pars0 = np.random.default_rng(seed).uniform(*zip(*bnds))
|
|
|
|
|
+ # use local optimizer
|
|
|
|
|
+ res = scipy.optimize.minimize(cost, pars0, bounds = bnds, method = "L-BFGS-B")
|
|
|
|
|
+ case "diff_evol":
|
|
|
|
|
+ res = scipy.optimize.differential_evolution(cost, bounds = bnds)
|
|
|
|
|
+ case "annel":
|
|
|
|
|
+ res = scipy.optimize.dual_annealing(cost, bounds = bnds)
|
|
|
|
|
+ case _:
|
|
|
|
|
+ assert False, f"{fname}::this method does not exist"
|
|
|
|
|
+
|
|
|
|
|
+ return {"pars": res.x, "success": res.success, "cost": res.fun}
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ MLE fitting of the Bayesian model:
|
|
|
|
|
+
|
|
|
|
|
+ pars_MLE = argmin_pars NLLF(pars| x, y)
|
|
|
|
|
+
|
|
|
|
|
+ Input:
|
|
|
|
|
+ x: array of floats
|
|
|
|
|
+ y: array of ints in {0,1}
|
|
|
|
|
+ method: string in ["local", "diff_evol", "anneal"], optimizer
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ {"pars": pars_MLE, "cost": NLLF at pars_MLE}
|
|
|
|
|
+ """
|
|
|
|
|
+ def fit(self, x, y, pars0 = None, method = "local"):
|
|
|
|
|
+ fname = "fit"
|
|
|
|
|
+
|
|
|
|
|
+ # defined nllf as function of parameters, data is already included
|
|
|
|
|
+ cost = lambda pars: self.nllf(x, y, pars)
|
|
|
|
|
+
|
|
|
|
|
+ # joint bounds of two groups
|
|
|
|
|
+ bnds = np.concatenate(self.bounds)
|
|
|
|
|
+
|
|
|
|
|
+ match method:
|
|
|
|
|
+ case "local":
|
|
|
|
|
+ # estimate initial guess of parameters (for local method)
|
|
|
|
|
+ if pars0 is None:
|
|
|
|
|
+ get_pars = lambda choice: self.fit_distr(x[y == choice], choice)["pars"]
|
|
|
|
|
+ pars0 = np.r_[get_pars(0), get_pars(1)]
|
|
|
|
|
+
|
|
|
|
|
+ # optimize using local optimizer
|
|
|
|
|
+ res = scipy.optimize.minimize(cost, pars0, bounds = bnds, method = "L-BFGS-B")
|
|
|
|
|
+ case "diff_evol":
|
|
|
|
|
+ res = scipy.optimize.differential_evolution(cost, bounds = bnds)
|
|
|
|
|
+ case "anneal":
|
|
|
|
|
+ res = scipy.optimize.dual_annealing(cost, bounds = bnds)
|
|
|
|
|
+ case _:
|
|
|
|
|
+ assert False, f"{fname}::this method does not exist"
|
|
|
|
|
+
|
|
|
|
|
+ return {"pars": res.x, "success": res.success, "cost": res.fun}
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Producing goodness of fit measures:
|
|
|
|
|
+
|
|
|
|
|
+ LLF = log_likelihood function
|
|
|
|
|
+ AIC = Akaike information criterion
|
|
|
|
|
+ BIC = Bayesian information criterion
|
|
|
|
|
+
|
|
|
|
|
+ Input:
|
|
|
|
|
+ x: array of floats
|
|
|
|
|
+ y: array of ints in {0,1}
|
|
|
|
|
+ pars: array of floats, model parameters
|
|
|
|
|
+ thresh: float, default 0.5, threshold value for classification
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ {"n": n, "k":k, "dof":n-k,
|
|
|
|
|
+ "LLF": log_likelihood,
|
|
|
|
|
+ "AIC": AIC,
|
|
|
|
|
+ "BIC": BIC,
|
|
|
|
|
+ "A": classification accuracy (threshold values = 0.5 prob)}
|
|
|
|
|
+ """
|
|
|
|
|
+ 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))
|
|
|
|
|
+ chi2 = np.sum(r**2)
|
|
|
|
|
+ p_val = scipy.stats.chi2.sf(chi2, dof)
|
|
|
|
|
+
|
|
|
|
|
+ # using model as classifier
|
|
|
|
|
+ matches = y == np.heaviside(p - thresh, 1)
|
|
|
|
|
+
|
|
|
|
|
+ return {"LLF": llf, "AIC": AIC, "BIC": BIC,
|
|
|
|
|
+ "A" : np.count_nonzero (matches)/n,
|
|
|
|
|
+ "chi2": chi2, "p-value(chi2)": p_val, # not very useful
|
|
|
|
|
+ "n": n, "k": k, "dof": dof}
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Generate m parameters via non-parametric bootstrapping with minimal constraint
|
|
|
|
|
+
|
|
|
|
|
+ bootstrapped sampled = (xb, yb) sampled with replacement from (x,y)
|
|
|
|
|
+
|
|
|
|
|
+ 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
|
|
|
|
|
+ method: string in ["local", "diff_evol", "anneal"], optimizer
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ array of m x len(pars) floats
|
|
|
|
|
+ """
|
|
|
|
|
+ def get_nonparam_boots_pars(self, x, y, m, seed = 1, method = "local"):
|
|
|
|
|
+ fname = "get_nonparam_boots_pars"
|
|
|
|
|
+
|
|
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
|
|
+
|
|
|
|
|
+ # discussing original data
|
|
|
|
|
+ res = self.fit(x, y, method = method)
|
|
|
|
|
+ assert res["success"], f"{fname}::fitting original data failed."
|
|
|
|
|
+
|
|
|
|
|
+ # generate bootstrapped parameters
|
|
|
|
|
+ pars = res["pars"]
|
|
|
|
|
+ lst = [pars]
|
|
|
|
|
+ n = len(x)
|
|
|
|
|
+
|
|
|
|
|
+ while True:
|
|
|
|
|
+
|
|
|
|
|
+ # sampling with replacement with restrictions
|
|
|
|
|
+ idx = rng.choice(n, n)
|
|
|
|
|
+ if np.sum(y[idx]) in [0, n]: continue
|
|
|
|
|
+
|
|
|
|
|
+ res = self.fit(x[idx], y[idx], pars0 = pars, method = method)
|
|
|
|
|
+ if not res["success"]: continue
|
|
|
|
|
+
|
|
|
|
|
+ lst.append(res["pars"])
|
|
|
|
|
+ if len(lst) == m: break
|
|
|
|
|
+
|
|
|
|
|
+ return np.array(lst)
|
|
|
|
|
+
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Generate m parameters via non-parametric stratified bootstrapping:
|
|
|
|
|
+
|
|
|
|
|
+ bootstrapped sampled = (xb, yb) sampled with replacement from (x,y) for each groups separately
|
|
|
|
|
+
|
|
|
|
|
+ meaning
|
|
|
|
|
+
|
|
|
|
|
+ xb = (sampled with replacement from x0, sampled with replacement from x1)
|
|
|
|
|
+ yb = (0 ... 0, 1 ... 1)
|
|
|
|
|
+ n0 n1
|
|
|
|
|
+
|
|
|
|
|
+ Note samples from each group in yb is constant and same as in y.
|
|
|
|
|
+
|
|
|
|
|
+ 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
|
|
|
|
|
+ method: string in ["local", "diff_evol", "anneal"], optimizer
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ array of m x len(pars) floats
|
|
|
|
|
+ """
|
|
|
|
|
+ def get_nonparam_strat_boots_pars(self, x, y, m, seed = 1, method = "local"):
|
|
|
|
|
+ fname = "get_nonparam_strat_boots_pars"
|
|
|
|
|
+
|
|
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
|
|
+
|
|
|
|
|
+ # discussing original data
|
|
|
|
|
+ res = self.fit(x, y, method = method)
|
|
|
|
|
+ assert res["success"], f"{fname}::fitting original data failed."
|
|
|
|
|
+
|
|
|
|
|
+ # separate data of both 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)])
|
|
|
|
|
+
|
|
|
|
|
+ # generate bootstrapped parameters
|
|
|
|
|
+ pars = res["pars"]
|
|
|
|
|
+ lst = [pars]
|
|
|
|
|
+ while True:
|
|
|
|
|
+ # stratified sampling with replacement
|
|
|
|
|
+ xb = np.concatenate([rng.choice(xs[i], ns[i]) for i in range(2)])
|
|
|
|
|
+
|
|
|
|
|
+ # do fitting
|
|
|
|
|
+ res = self.fit(xb, yb, pars0 = pars, method = method)
|
|
|
|
|
+ if not res["success"]: continue
|
|
|
|
|
+
|
|
|
|
|
+ lst.append(res["pars"])
|
|
|
|
|
+ if len(lst) == m: break
|
|
|
|
|
+
|
|
|
|
|
+ return np.array(lst)
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Generate m parameters via parametric bootstrapping:
|
|
|
|
|
+
|
|
|
|
|
+ bootstrapped sample = (x, yb) yb ~ B(ymodel)
|
|
|
|
|
+
|
|
|
|
|
+ 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
|
|
|
|
|
+ method:
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ array of m x len(pars) floats
|
|
|
|
|
+ """
|
|
|
|
|
+ def get_param_boots_pars(self, x, y, m, seed = 1, method = "local"):
|
|
|
|
|
+ fname = "get_param_boots_pars"
|
|
|
|
|
+
|
|
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
|
|
+
|
|
|
|
|
+ # discussing original data
|
|
|
|
|
+ res = self.fit(x, y, method = method)
|
|
|
|
|
+ assert res["success"], f"{fname}::fitting original data failed."
|
|
|
|
|
+
|
|
|
|
|
+ # calculate predicted conditional probabilities
|
|
|
|
|
+ pars = res["pars"]
|
|
|
|
|
+ p = self.model(x, pars)
|
|
|
|
|
+
|
|
|
|
|
+ # generate bootstrapped parameters
|
|
|
|
|
+ lst = [pars]
|
|
|
|
|
+ n = len(x)
|
|
|
|
|
+
|
|
|
|
|
+ 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
|
|
|
|
|
+
|
|
|
|
|
+ # fit and get new parameter
|
|
|
|
|
+ res = self.fit(x, y_sim, pars0 = pars, method=method)
|
|
|
|
|
+ if not res["success"]: continue
|
|
|
|
|
+
|
|
|
|
|
+ lst.append(res["pars"])
|
|
|
|
|
+ if len(lst) == m: break
|
|
|
|
|
+
|
|
|
|
|
+ return np.array(lst)
|
|
|
|
|
+
|
|
|
|
|
+ """
|
|
|
|
|
+ Generate m parameters via parametric stratified bootstrapping by
|
|
|
|
|
+ sampling x from parametrized distributions associated to individual groups:
|
|
|
|
|
+
|
|
|
|
|
+ bootstrapped sample = (xb, yb)
|
|
|
|
|
+ xb = (sampled from distr for x0, sampled from distr for x1)
|
|
|
|
|
+ yb = (0 ... 0, 1 ... 1)
|
|
|
|
|
+ n0 n1
|
|
|
|
|
+
|
|
|
|
|
+ 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
|
|
|
|
|
+ method:
|
|
|
|
|
+
|
|
|
|
|
+ Return:
|
|
|
|
|
+ array of m x len(pars) floats
|
|
|
|
|
+ """
|
|
|
|
|
+ def get_param_strat_boots_pars(self, x, y, m, seed = 1, method = "local"):
|
|
|
|
|
+ fname = "get_param_strat_boots_pars"
|
|
|
|
|
+
|
|
|
|
|
+ assert self.distr_sample is not None, f"{fname}::distr_sample is not defined"
|
|
|
|
|
+
|
|
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
|
|
+
|
|
|
|
|
+ # separate data of both groups
|
|
|
|
|
+ xs = [x[y == i] for i in range(2)]
|
|
|
|
|
+ ns = [len(e) for e in xs]
|
|
|
|
|
+
|
|
|
|
|
+ # separate data of both groups
|
|
|
|
|
+ pars_g = np.concatenate([self.fit_distr(e, i)["pars"] for i, e in enumerate(xs)])
|
|
|
|
|
+
|
|
|
|
|
+ # discussing original data
|
|
|
|
|
+ res = self.fit(x, y, method = method)
|
|
|
|
|
+ assert res["success"], f"{fname}::fitting original data failed."
|
|
|
|
|
+
|
|
|
|
|
+ # common vector states
|
|
|
|
|
+ yb = np.concatenate([np.full(ns[i], i) for i in range(2)])
|
|
|
|
|
+
|
|
|
|
|
+ # generate bootstrapped parameters
|
|
|
|
|
+ pars = res["pars"]
|
|
|
|
|
+ lst = [pars]
|
|
|
|
|
+
|
|
|
|
|
+ while True:
|
|
|
|
|
+
|
|
|
|
|
+ # Generate new sample of points for each group
|
|
|
|
|
+ xb = self.distr_sample(rng, pars_g, ns)
|
|
|
|
|
+
|
|
|
|
|
+ # fit and get new parameter
|
|
|
|
|
+ res = self.fit(xb, yb, pars0 = pars, method = method)
|
|
|
|
|
+ if not res["success"]: continue
|
|
|
|
|
+
|
|
|
|
|
+ lst.append(res["pars"])
|
|
|
|
|
+ if len(lst) == m: break
|
|
|
|
|
+
|
|
|
|
|
+ return np.array(lst)
|
