|
|
@@ -2,7 +2,6 @@ import numpy as np
|
|
|
import scipy
|
|
|
import scipy.stats
|
|
|
|
|
|
-
|
|
|
"""
|
|
|
Model function
|
|
|
|
|
|
@@ -178,6 +177,8 @@ def logit_poly_pars_quantiles_normal(probs, mean_b, cov):
|
|
|
|
|
|
b ~ N(mean_b, cov)
|
|
|
|
|
|
+ This distribution is asymptotic MLE distribution of parameters.
|
|
|
+
|
|
|
Input:
|
|
|
x: array of n float
|
|
|
probs: array of m floats, probabilities
|
|
|
@@ -212,13 +213,14 @@ def logit_poly_model_quantiles_normal(x, probs, mean_b, cov):
|
|
|
F(x|b) = sum_{i=0}^degree x^i b_i
|
|
|
|
|
|
at given probabilities p and values x assuming
|
|
|
- normal distribution of parameters:
|
|
|
+ normal distribution of parameters :
|
|
|
|
|
|
b ~ N(mean_b, cov)
|
|
|
|
|
|
+ This distribution is asymptotic MLE distribution of parameters.
|
|
|
We approximate exact model with linear expansion
|
|
|
|
|
|
- p(x|b) = p(x|b_mean) + dp/db (x| b_mean) (b - b_mean)
|
|
|
+ p(x|b) = p(x|b_mean) + dp/db(x| b_mean) (b - b_mean)
|
|
|
|
|
|
and the last term is normally distributed.
|
|
|
|
|
|
@@ -242,4 +244,127 @@ def logit_poly_model_quantiles_delta(x, probs, mean_b, cov):
|
|
|
# computing quantiles of logit
|
|
|
Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)
|
|
|
|
|
|
- return np.clip(Q, a_min = 0, a_max = 1)
|
|
|
+ return np.clip(Q, a_min = 0, a_max = 1)
|
|
|
+
|
|
|
+
|
|
|
+"""
|
|
|
+ Generate m parameters via non-parametric bootstrapping with a minimal constraint
|
|
|
+ that both groups should be present in the sampled data.
|
|
|
+
|
|
|
+ Input:
|
|
|
+ lm: linear_model.LogisticRegression
|
|
|
+ x : array of n floats
|
|
|
+ y : array of n int in {0,1}
|
|
|
+ m: integer, number of samples
|
|
|
+ degree: int, degree of decision function
|
|
|
+ seed : int, seed for the random generator
|
|
|
+
|
|
|
+ Return:
|
|
|
+ array of mx(degree + 1)
|
|
|
+
|
|
|
+ Return:
|
|
|
+ array of mx(degree + 1)
|
|
|
+"""
|
|
|
+def get_nonparam_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
|
|
|
+
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
+
|
|
|
+ # fitting original data
|
|
|
+ pars = logit_poly_fit(lm, x, y, degree = degree)
|
|
|
+
|
|
|
+ n = len(x)
|
|
|
+
|
|
|
+ # generate parameters
|
|
|
+ lst = [pars]
|
|
|
+ while True:
|
|
|
+
|
|
|
+ # create set indices for sampling with replacement + constraint
|
|
|
+ idx = rng.choice(n, n)
|
|
|
+ if np.sum(y[idx]) in [0, n]: continue
|
|
|
+
|
|
|
+ lst.append(logit_poly_fit(lm, x[idx], y[idx], degree=degree))
|
|
|
+ if len(lst) == m: break
|
|
|
+
|
|
|
+ return np.array(lst)
|
|
|
+
|
|
|
+
|
|
|
+"""
|
|
|
+ Generate m parameters via non-parametric stratified bootstrapping.
|
|
|
+
|
|
|
+ Input:
|
|
|
+ lm: linear_model.LogisticRegression
|
|
|
+ x : array of n floats
|
|
|
+ y : array of n int in {0,1}
|
|
|
+ m: integer, number of samples
|
|
|
+ degree: int, degree of decision function
|
|
|
+ seed : int, seed for the random generator
|
|
|
+
|
|
|
+ Return:
|
|
|
+ array of mx(degree + 1)
|
|
|
+"""
|
|
|
+def get_nonparam_stratified_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
|
|
|
+
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
+
|
|
|
+ # pars of original data
|
|
|
+ pars = logit_poly_fit(lm, x, y, degree = degree)
|
|
|
+
|
|
|
+ # 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)])
|
|
|
+
|
|
|
+ # generate parameters
|
|
|
+ lst = [pars]
|
|
|
+ for _ in range(m):
|
|
|
+ # stratified sampling with replacement
|
|
|
+ xb = np.concatenate([rng.choice(xs[i], ns[i]) for i in range(2)])
|
|
|
+ # do fitting
|
|
|
+ lst.append(logit_poly_fit(lm, xb, yb, degree = degree))
|
|
|
+
|
|
|
+ return np.array(lst)
|
|
|
+
|
|
|
+"""
|
|
|
+ Generate m parameters via parametric bootstrapping.
|
|
|
+
|
|
|
+ Input:
|
|
|
+ lm: linear_model.LogisticRegression
|
|
|
+ x : array of n floats
|
|
|
+ y : array of n int in {0,1}
|
|
|
+ m: integer, number of samples
|
|
|
+ degree: int, degree of decision function
|
|
|
+ seed : int, seed for the random generator
|
|
|
+
|
|
|
+ Return:
|
|
|
+ array of mx(degree + 1)
|
|
|
+
|
|
|
+ Ref:
|
|
|
+
|
|
|
+ * https://www.scirp.org/journal/paperinformation?paperid=70962
|
|
|
+"""
|
|
|
+def get_parametric_boots_pars(lm, x, y, m, degree = 1, seed = 1977):
|
|
|
+
|
|
|
+ # first discuss original dataset
|
|
|
+ pars = logit_poly_fit(lm, x, y, degree = degree)
|
|
|
+ p = logit_poly_model(x, pars)
|
|
|
+
|
|
|
+ rng = np.random.default_rng(seed)
|
|
|
+
|
|
|
+ n = len(x)
|
|
|
+
|
|
|
+ # generate parameters
|
|
|
+ lst = [pars]
|
|
|
+ while True:
|
|
|
+
|
|
|
+ # Generate new binary outcomes from Bernoulli(p_i)
|
|
|
+ y_sim = np.random.binomial(n = 1, p = p)
|
|
|
+
|
|
|
+ if np.sum(y_sim) in [0, n]: continue
|
|
|
+
|
|
|
+ lst.append(logit_poly_fit(lm, x, y_sim, degree = degree))
|
|
|
+
|
|
|
+ if len(lst) == m: break
|
|
|
+
|
|
|
+ return np.array(lst)
|
