Skip to content

Fix the SEM criterion and link fit #11

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
adimajo merged 5 commits into from
Nov 1, 2026
Merged

Fix the SEM criterion and link fit #11

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b0696c8
Chain the criterion branches in SEM fit
adimajo Nov 1, 2026
8779af3
Use the training size in the BIC penalty
adimajo Nov 1, 2026
dccd572
Guard the continuous link against one class
adimajo Nov 1, 2026
1c9c8fe
Rewrite the tail of y_p as a named helper
adimajo Nov 1, 2026
71c9f89
Score the validation likelihood on validation rows
adimajo Nov 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 63 additions & 35 deletions glmdisc/_fitSEM.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -160,22 +187,23 @@ 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])))

# On calcule y | q_{-j} , q_j
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(
Expand Down Expand Up @@ -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]))))

Expand Down
68 changes: 67 additions & 1 deletion tests/tests_SEM/test_SEM_fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down
Loading