From b0696c808bf2f2968ae35053c65e32268ff41da4 Mon Sep 17 00:00:00 2001 From: adimajo Date: Sun, 1 Nov 2026 13:55:36 +0100 Subject: [PATCH 1/5] Chain the criterion branches in SEM fit The five independent ifs are one if/elif chain, so exactly one branch sets the performance and the Gini branch computes it once, on the validation or training rows. --- glmdisc/_fitSEM.py | 32 ++++++++++++-------------------- 1 file changed, 12 insertions(+), 20 deletions(-) diff --git a/glmdisc/_fitSEM.py b/glmdisc/_fitSEM.py index 957b193..5bba1b1 100644 --- a/glmdisc/_fitSEM.py +++ b/glmdisc/_fitSEM.py @@ -38,28 +38,20 @@ def _calculate_criterion(self, emap, model_emap, current_encoder_emap): if self.validation: performance = loglik logger.info("Current likelihood on validation 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)) - - 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)) - - if self.criterion == 'gini' and not self.validation: + 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(self.n) * model_emap.coef_.shape[1] - 2 * loglik) + logger.info("Current BIC on train set: " + str(- performance)) + else: + rows = self.validation_rows if self.validation else self.train_rows performance = sk.metrics.roc_auc_score( - self.labels[self.train_rows], model_emap.predict_proba( + self.labels[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)) + emap[rows, :].astype(str)))[:, 1:]) + logger.info("Current Gini on " + ("validation" if self.validation else "training") + " set: " + + str(performance)) return performance From 8779af3b2fc6cbb28868b11aefce55cc987aef80 Mon Sep 17 00:00:00 2001 From: adimajo Date: Sun, 1 Nov 2026 13:55:36 +0100 Subject: [PATCH 2/5] Use the training size in the BIC penalty The log-likelihood is computed on the training rows, so the penalty log(n) k must use their number, not the size of the whole sample. --- glmdisc/_fitSEM.py | 2 +- tests/tests_SEM/test_SEM_fit.py | 20 +++++++++++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/glmdisc/_fitSEM.py b/glmdisc/_fitSEM.py index 5bba1b1..766710e 100644 --- a/glmdisc/_fitSEM.py +++ b/glmdisc/_fitSEM.py @@ -42,7 +42,7 @@ def _calculate_criterion(self, emap, model_emap, current_encoder_emap): performance = -(2 * model_emap.coef_.shape[1] - 2 * loglik) logger.info("Current AIC on train set: " + str(- performance)) else: - performance = -(log(self.n) * model_emap.coef_.shape[1] - 2 * loglik) + performance = -(log(len(self.train_rows)) * model_emap.coef_.shape[1] - 2 * loglik) logger.info("Current BIC on train set: " + str(- performance)) else: rows = self.validation_rows if self.validation else self.train_rows diff --git a/tests/tests_SEM/test_SEM_fit.py b/tests/tests_SEM/test_SEM_fit.py index 41c0263..7bf7366 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,24 @@ 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_nan(): n = 100 d = 2 From dccd572ec8787776062f1f25efea391ad24695e8 Mon Sep 17 00:00:00 2001 From: adimajo Date: Sun, 1 Nov 2026 13:56:12 +0100 Subject: [PATCH 3/5] Guard the continuous link against one class When every complete training row of a feature carries the same label, the logistic regression cannot be fitted and the chain crashed. Use a constant link with probability one for that label instead. --- glmdisc/_fitSEM.py | 42 ++++++++++++++++++++++++++------- tests/tests_SEM/test_SEM_fit.py | 26 ++++++++++++++++++++ 2 files changed, 60 insertions(+), 8 deletions(-) diff --git a/glmdisc/_fitSEM.py b/glmdisc/_fitSEM.py index 766710e..6e59b6c 100644 --- a/glmdisc/_fitSEM.py +++ b/glmdisc/_fitSEM.py @@ -56,6 +56,35 @@ def _calculate_criterion(self, emap, model_emap, current_encoder_emap): return 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 + """ + + def fit(self, X, y, sample_weight=None): + self.classes_ = np.unique(y) + return self + + 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 _init_fit_sem(self, edisc): # Initialization for following the performance of the discretization current_best = 0 @@ -83,11 +112,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 @@ -152,9 +177,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]))) diff --git a/tests/tests_SEM/test_SEM_fit.py b/tests/tests_SEM/test_SEM_fit.py index 7bf7366..771e177 100644 --- a/tests/tests_SEM/test_SEM_fit.py +++ b/tests/tests_SEM/test_SEM_fit.py @@ -187,6 +187,32 @@ def test_calculate_criterion_bic_penalty_uses_train_size(): 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_nan(): n = 100 d = 2 From 1c9c8fe45a0dfa16fda8f808124ab46c9fbd02de Mon Sep 17 00:00:00 2001 From: adimajo Date: Sun, 1 Nov 2026 13:57:32 +0100 Subject: [PATCH 4/5] Rewrite the tail of y_p as a named helper p * (2y - 1) - y + 1 is p(y=1) where y=1 and 1 - p elsewhere. Say so with np.where in one helper shared by the continuous and categorical branches. --- glmdisc/_fitSEM.py | 22 ++++++++++++++++------ tests/tests_SEM/test_SEM_fit.py | 6 ++++++ 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/glmdisc/_fitSEM.py b/glmdisc/_fitSEM.py index 6e59b6c..25fecf0 100644 --- a/glmdisc/_fitSEM.py +++ b/glmdisc/_fitSEM.py @@ -85,6 +85,16 @@ def _fit_link(x, y): 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): # Initialization for following the performance of the discretization current_best = 0 @@ -188,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( @@ -238,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 771e177..5e2648d 100644 --- a/tests/tests_SEM/test_SEM_fit.py +++ b/tests/tests_SEM/test_SEM_fit.py @@ -213,6 +213,12 @@ def test_fit_sem_collapsed_continuous_feature(): 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 From 71c9f89919cf70b0601807d884463b9963813442 Mon Sep 17 00:00:00 2001 From: adimajo Date: Sun, 1 Nov 2026 14:22:23 +0100 Subject: [PATCH 5/5] Score the validation likelihood on validation rows With validation=True the aic and bic criteria logged a validation likelihood but computed it on the training rows, so the best iteration was picked on training data. Use the validation rows. --- glmdisc/_fitSEM.py | 6 +++--- tests/tests_SEM/test_SEM_fit.py | 16 ++++++++++++++++ 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/glmdisc/_fitSEM.py b/glmdisc/_fitSEM.py index 25fecf0..da02a23 100644 --- a/glmdisc/_fitSEM.py +++ b/glmdisc/_fitSEM.py @@ -29,11 +29,12 @@ 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 @@ -45,7 +46,6 @@ def _calculate_criterion(self, emap, model_emap, current_encoder_emap): performance = -(log(len(self.train_rows)) * model_emap.coef_.shape[1] - 2 * loglik) logger.info("Current BIC on train set: " + str(- performance)) else: - rows = self.validation_rows if self.validation else self.train_rows performance = sk.metrics.roc_auc_score( self.labels[rows], model_emap.predict_proba( X=current_encoder_emap.transform( diff --git a/tests/tests_SEM/test_SEM_fit.py b/tests/tests_SEM/test_SEM_fit.py index 5e2648d..f0ddcc3 100644 --- a/tests/tests_SEM/test_SEM_fit.py +++ b/tests/tests_SEM/test_SEM_fit.py @@ -187,6 +187,22 @@ def test_calculate_criterion_bic_penalty_uses_train_size(): 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))