diff --git a/glmdisc/_fitSEM.py b/glmdisc/_fitSEM.py index 957b193..da02a23 100644 --- a/glmdisc/_fitSEM.py +++ b/glmdisc/_fitSEM.py @@ -29,39 +29,70 @@ def _calculate_criterion(self, emap, model_emap, current_encoder_emap): ------- criterion value """ + rows = self.validation_rows if self.validation else self.train_rows if self.criterion in ['aic', 'bic']: - loglik = -sk.metrics.log_loss(self.labels[self.train_rows], + loglik = -sk.metrics.log_loss(self.labels[rows], model_emap.predict_proba( X=current_encoder_emap.transform( - emap[self.train_rows, :].astype(str))), + emap[rows, :].astype(str))), normalize=False) if self.validation: performance = loglik logger.info("Current likelihood on validation set: " + str(performance)) + elif self.criterion == 'aic': + performance = -(2 * model_emap.coef_.shape[1] - 2 * loglik) + logger.info("Current AIC on train set: " + str(- performance)) + else: + performance = -(log(len(self.train_rows)) * model_emap.coef_.shape[1] - 2 * loglik) + logger.info("Current BIC on train set: " + str(- performance)) + else: + performance = sk.metrics.roc_auc_score( + self.labels[rows], model_emap.predict_proba( + X=current_encoder_emap.transform( + emap[rows, :].astype(str)))[:, 1:]) + logger.info("Current Gini on " + ("validation" if self.validation else "training") + " set: " + + str(performance)) - if self.criterion == 'aic' and not self.validation: - performance = -(2 * model_emap.coef_.shape[1] - 2 * loglik) - logger.info("Current AIC on train set: " + str(- performance)) + return performance - if self.criterion == 'bic' and not self.validation: - performance = -(log(self.n) * model_emap.coef_.shape[1] - 2 * loglik) - logger.info("Current BIC on train set: " + str(- performance)) - if self.criterion == 'gini' and self.validation: - performance = sk.metrics.roc_auc_score( - self.labels[self.validation_rows], model_emap.predict_proba( - X=current_encoder_emap.transform( - emap[self.validation_rows, :].astype(str)))[:, 1:]) - logger.info("Current Gini on validation set: " + str(performance)) +class _SingleClassLink(sk.linear_model.LogisticRegression): + """ + Link p(e|x) when all complete training rows carry the same label, which a logistic regression cannot fit + """ - if self.criterion == 'gini' and not self.validation: - performance = sk.metrics.roc_auc_score( - self.labels[self.train_rows], model_emap.predict_proba( - X=current_encoder_emap.transform( - emap[self.train_rows, :].astype(str)))[:, 1:]) - logger.info("Current Gini on training set: " + str(performance)) + def fit(self, X, y, sample_weight=None): + self.classes_ = np.unique(y) + return self - return performance + def predict_proba(self, X): + return np.ones((X.shape[0], 1)) + + +def _new_link(): + return sk.linear_model.LogisticRegression(C=1e40, + solver='newton-cg', + max_iter=25, + tol=0.001, + warm_start=False) + + +def _fit_link(x, y): + """ + Fit p(e|x) of one continuous feature, falling back to a constant when y has a single label + """ + link = _SingleClassLink() if np.unique(y).size < 2 else _new_link() + return link.fit(X=x, y=y) + + +def _likelihood_of_labels(proba, labels): + """ + Probability given to the observed label of each row: p(y = 1) where y is 1, p(y = 0) elsewhere + + :param numpy.ndarray proba: predict_proba output of the logistic regression, columns (y = 0, y = 1) + :param numpy.ndarray labels: observed labels + """ + return np.where(np.ravel(labels) == 1, proba[:, 1], proba[:, 0]) def _init_fit_sem(self, edisc): @@ -91,11 +122,7 @@ def _init_fit_sem(self, edisc): m = [None] * (self.d_cont + self.d_qual) for j in range(self.d_cont): - link[j] = sk.linear_model.LogisticRegression(C=1e40, - solver='newton-cg', - max_iter=25, - tol=0.001, - warm_start=False) + link[j] = _new_link() return current_best, emap, model_edisc, model_emap, current_encoder_edisc, current_encoder_emap, link, m @@ -160,9 +187,10 @@ def _fit_sem(self, edisc, predictors_trans, continu_complete_case, **kwargs): # On commence par les quantitatives if j < self.d_cont: # On apprend q_j | x_j - link[j].fit(y=edisc[self.train_rows, :][continu_complete_case[self.train_rows, j], j], - X=self.predictors_cont[self.train_rows, :][continu_complete_case[self.train_rows, - j], j].reshape(-1, 1)) + link[j] = _fit_link( + x=self.predictors_cont[self.train_rows, :][continu_complete_case[self.train_rows, + j], j].reshape(-1, 1), + y=edisc[self.train_rows, :][continu_complete_case[self.train_rows, j], j]) y_p = np.zeros((self.n, len(m[j]))) @@ -170,12 +198,12 @@ def _fit_sem(self, edisc, predictors_trans, continu_complete_case, **kwargs): for k in range(len(m[j])): modalites = np.zeros((self.n, len(m[j]))) modalites[:, k] = np.ones((self.n, )) - y_p[:, k] = model_edisc.predict_proba(np.column_stack( + y_p[:, k] = _likelihood_of_labels(model_edisc.predict_proba(np.column_stack( (base_disjonctive[:, 0:(sum(list(map(len, m[0:j]))))], modalites, base_disjonctive[:, (sum(list(map(len, - m[0:(j + 1)])))):(sum(list(map(len, m))))])) - )[:, 1] * (2 * np.ravel(self.labels) - 1) - np.ravel(self.labels) + 1 + m[0:(j + 1)])))):(sum(list(map(len, m))))]))), + self.labels) # On calcule q_j | x_j sur tout le monde t = link[j].predict_proba( @@ -220,12 +248,12 @@ def _fit_sem(self, edisc, predictors_trans, continu_complete_case, **kwargs): modalites = np.zeros((self.n, len(m[j]))) modalites[:, k] = np.ones((self.n, )) - y_p[:, k] = model_edisc.predict_proba(np.column_stack( + y_p[:, k] = _likelihood_of_labels(model_edisc.predict_proba(np.column_stack( (base_disjonctive[:, 0:(sum(list(map(len, m[0:j]))))], modalites, base_disjonctive[:, - (sum(list(map(len, m[0:(j + 1)])))):(sum(list(map(len, m))))])) - )[:, 1] * (2 * np.ravel(self.labels) - 1) - np.ravel(self.labels) + 1 + (sum(list(map(len, m[0:(j + 1)])))):(sum(list(map(len, m))))]))), + self.labels) t = np.zeros((self.n, int(len(m[j])))) diff --git a/tests/tests_SEM/test_SEM_fit.py b/tests/tests_SEM/test_SEM_fit.py index 41c0263..f0ddcc3 100644 --- a/tests/tests_SEM/test_SEM_fit.py +++ b/tests/tests_SEM/test_SEM_fit.py @@ -158,7 +158,7 @@ def test_calculate_criterion(): np.random.seed(1) model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11) assert math.isclose(glmdisc._fitSEM._calculate_criterion(model, emap, model_emap, current_encoder_emap), - modele_bic + (math.log(model.n) - 2) * model_emap.coef_.shape[1]) + modele_bic + (math.log(len(model.train_rows)) - 2) * model_emap.coef_.shape[1]) model = glmdisc.Glmdisc(criterion="gini") model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11) @@ -169,6 +169,72 @@ def test_calculate_criterion(): assert 0 <= glmdisc._fitSEM._calculate_criterion(model, emap, model_emap, current_encoder_emap) <= 1 +def test_calculate_criterion_bic_penalty_uses_train_size(): + n = 100 + x, y, theta = glmdisc.Glmdisc.generate_data(n, 2) + np.random.seed(1) + model = glmdisc.Glmdisc(criterion="bic", validation=False) + model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11) + assert len(model.train_rows) < n + emap = np.random.randint(0, 3, size=(n, 2)) + encoder = sk.preprocessing.OneHotEncoder().fit(X=emap.astype(str)) + reglog = sk.linear_model.LogisticRegression(C=1e40, tol=0.001, max_iter=25) + reglog.fit(encoder.transform(emap[model.train_rows].astype(str)), y[model.train_rows]) + loglik = -sk.metrics.log_loss(y[model.train_rows], + reglog.predict_proba(encoder.transform(emap[model.train_rows].astype(str))), + normalize=False) + expected = -(math.log(len(model.train_rows)) * reglog.coef_.shape[1] - 2 * loglik) + assert math.isclose(glmdisc._fitSEM._calculate_criterion(model, emap, reglog, encoder), expected) + + +def test_calculate_criterion_validation_likelihood_uses_validation_rows(): + n = 100 + x, y, theta = glmdisc.Glmdisc.generate_data(n, 2) + np.random.seed(1) + model = glmdisc.Glmdisc(criterion="bic", validation=True) + model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11) + emap = np.random.randint(0, 3, size=(n, 2)) + encoder = sk.preprocessing.OneHotEncoder().fit(X=emap.astype(str)) + reglog = sk.linear_model.LogisticRegression(C=1e40, tol=0.001, max_iter=25) + reglog.fit(encoder.transform(emap[model.train_rows].astype(str)), y[model.train_rows]) + rows = model.validation_rows + expected = -sk.metrics.log_loss(y[rows], reglog.predict_proba(encoder.transform(emap[rows].astype(str))), + normalize=False) + assert math.isclose(glmdisc._fitSEM._calculate_criterion(model, emap, reglog, encoder), expected) + + +def test_fit_link_single_class(): + x = np.linspace(0, 1, 20).reshape(-1, 1) + link = glmdisc._fitSEM._fit_link(x, np.full(20, 3)) + assert isinstance(link, sk.linear_model.LogisticRegression) + np.testing.assert_array_equal(link.classes_, [3]) + np.testing.assert_array_equal(link.predict_proba(x), np.ones((20, 1))) + + link = glmdisc._fitSEM._fit_link(x, np.arange(20) % 2) + assert link.predict_proba(x).shape == (20, 2) + + +def test_fit_sem_collapsed_continuous_feature(): + np.random.seed(1) + n = 100 + x, y, theta = glmdisc.Glmdisc.generate_data(n, 2) + model = glmdisc.Glmdisc(criterion="bic", validation=False, test=False, m_start=2) + model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11) + model.predict(predictors_cont=x, predictors_qual=None) + edisc = np.zeros((n, 2), dtype=int) + edisc[:, 1] = np.random.randint(0, 2, size=n) + model._calculate_shape() + complete = np.ones((n, 2), dtype=bool) + glmdisc._fitSEM._fit_sem(model, edisc, None, complete) + assert isinstance(model.best_link[0], glmdisc._fitSEM._SingleClassLink) + + +def test_likelihood_of_labels(): + proba = np.array([[0.7, 0.3], [0.2, 0.8], [0.6, 0.4], [0.1, 0.9]]) + for labels in (np.array([0, 1, 1, 0]), np.array([[0], [1], [1], [0]]), np.array([False, True, True, False])): + np.testing.assert_allclose(glmdisc._fitSEM._likelihood_of_labels(proba, labels), [0.7, 0.8, 0.4, 0.1]) + + def test_nan(): n = 100 d = 2