diff --git a/glmdisc/__init__.py b/glmdisc/__init__.py
index b65229e..7accb6e 100644
--- a/glmdisc/__init__.py
+++ b/glmdisc/__init__.py
@@ -215,8 +215,14 @@ def __init__(self, algorithm="SEM", test=True, validation=True, criterion="bic",
             raise ValueError(msg)
 
         # m_start doit être pas déconnant
-        if not 2 <= m_start <= 50:
-            msg = 'Please set 2 <= m_start <= 50'
+        if not isinstance(m_start, (int, np.integer)) or not 2 <= m_start <= 50:
+            msg = 'Please set an integer 2 <= m_start <= 50'
+            logger.error(msg)
+            raise ValueError(msg)
+
+        # burn_in must be a non-negative integer
+        if not isinstance(burn_in, (int, np.integer)) or isinstance(burn_in, bool) or burn_in < 0:
+            msg = 'Please set an integer burn_in >= 0'
             logger.error(msg)
             raise ValueError(msg)
 
@@ -266,6 +272,8 @@ def _check_is_fitted(self):
         """
         if self.algorithm == "SEM":
             try:
+                if self.best_reglog is None:
+                    raise sk.exceptions.NotFittedError("This Glmdisc instance is not fitted yet.")
                 sk.utils.validation.check_is_fitted(self.best_reglog)
                 for link in self.best_link:
                     if isinstance(link, sk.linear_model.LogisticRegression):
diff --git a/glmdisc/_bestFormula.py b/glmdisc/_bestFormula.py
index 130736b..bb8ea1c 100644
--- a/glmdisc/_bestFormula.py
+++ b/glmdisc/_bestFormula.py
@@ -16,9 +16,14 @@ def _best_formula_cont_sem(self, emap_best):
     for j in range(self.d_cont):
         best_disc.append([])
         for k in np.unique(emap_best[:, j]):
-            best_disc[j].append(
-                np.nanmin(self.predictors_cont[emap_best[:, j] == k, j]))
-        del best_disc[j][np.where(best_disc[j] == np.nanmin(best_disc[j]))[0][0]]
+            lower_bound = self.predictors_cont[emap_best[:, j] == k, j]
+            lower_bound = lower_bound[~np.isnan(lower_bound)]
+            # the level reserved for missing values has no lower bound
+            if lower_bound.size > 0:
+                best_disc[j].append(lower_bound.min())
+        # the smallest lower bound is not a cut-point
+        if best_disc[j]:
+            del best_disc[j][np.argmin(best_disc[j])]
         if len(best_disc[j]) < 1:
             logger.info("No cut-points found for continuous variable " + str(j) + ".")
         else:
diff --git a/glmdisc/_discretize.py b/glmdisc/_discretize.py
index 86ca867..c6e9c05 100644
--- a/glmdisc/_discretize.py
+++ b/glmdisc/_discretize.py
@@ -36,7 +36,7 @@ def _check_args_discretize_nn(self, predictors_cont, predictors_qual):
         raise ValueError(msg)
     if d_2 != self.d_qual:
         msg = ('Shape of ' + str(d_2) +
-               ' for predictors_cont does not match training set of size ' + str(self.d_qual) + '.')
+               ' for predictors_qual does not match training set of size ' + str(self.d_qual) + '.')
         logger.error(msg)
         raise ValueError(msg)
 
@@ -70,7 +70,7 @@ def _check_args_discretize_sem(self, predictors_cont, predictors_qual):
         raise ValueError(msg)
     if d_2 != sum(d_2bis) or d_2 != self.d_qual:
         msg = ('Shape of ' + str(d_2) +
-               ' for predictors_cont does not match provided link function '
+               ' for predictors_qual does not match provided link function '
                'of size ' + str(sum(d_2bis)) + ' and/or training set of size ' + str(self.d_qual) + '.')
         logger.error(msg)
         raise ValueError(msg)
diff --git a/glmdisc/_fit.py b/glmdisc/_fit.py
index 3f36c7c..5240086 100644
--- a/glmdisc/_fit.py
+++ b/glmdisc/_fit.py
@@ -5,7 +5,6 @@
 import numpy as np
 import sklearn as sk
 from loguru import logger
-from scipy import stats
 
 from glmdisc._fitNN import _fit_nn
 from glmdisc._fitSEM import _fit_sem
@@ -96,7 +95,7 @@ def _init_disc(self, continu_complete_case):
     initial quantization and transformation of categorical features to integers
     """
     self.affectations = [None] * (self.d_cont + self.d_qual)
-    edisc = np.random.choice(list(range(self.m_start)), size=(self.n, self.d_cont + self.d_qual))
+    edisc = np.random.randint(0, self.m_start, size=(self.n, self.d_cont + self.d_qual))
 
     for j in range(self.d_cont):
         edisc[np.invert(continu_complete_case[:, j]), j] = self.m_start
@@ -105,16 +104,9 @@ def _init_disc(self, continu_complete_case):
 
     for j in range(self.d_qual):
         self.affectations[j + self.d_cont] = sk.preprocessing.LabelEncoder().fit(self.predictors_qual[:, j])
-        if (self.m_start > stats.describe(sk.preprocessing.LabelEncoder().fit(self.predictors_qual[:, j]).transform(
-                self.predictors_qual[:, j])).minmax[1] + 1):
-            edisc[:, j + self.d_cont] = np.random.choice(list(range(
-                stats.describe(sk.preprocessing.LabelEncoder().fit(
-                    self.predictors_qual[:, j]).transform(
-                        self.predictors_qual[:, j])).minmax[1])),
-                size=self.n)
-        else:
-            edisc[:, j + self.d_cont] = np.random.choice(list(range(self.m_start)),
-                                                         size=self.n)
+        # a feature cannot start with more groups than it has levels
+        m_init = min(self.m_start, len(self.affectations[j + self.d_cont].classes_))
+        edisc[:, j + self.d_cont] = np.random.randint(0, m_init, size=self.n)
 
         predictors_trans[:, j] = (self.affectations[j + self.d_cont].transform(
             self.predictors_qual[:, j])).astype(int)
diff --git a/glmdisc/_fitNN.py b/glmdisc/_fitNN.py
index 8ed753f..e058359 100644
--- a/glmdisc/_fitNN.py
+++ b/glmdisc/_fitNN.py
@@ -12,7 +12,6 @@
 import tensorflow as tf
 import tensorflow.keras.optimizers
 from loguru import logger
-from scipy import stats
 from tensorflow.keras import Model
 from tensorflow.keras.callbacks import Callback, ReduceLROnPlateau
 from tensorflow.keras.layers import Dense
@@ -152,8 +151,7 @@ def _initialize_neural_net(self):
     m_qual = [self.m_start] * (self.d_cont + self.d_qual)
 
     for j in range(self.d_qual):
-        num_levels = stats.describe(sk.preprocessing.LabelEncoder().fit(self.predictors_qual[:, j]).transform(
-            self.predictors_qual[:, j])).minmax[1]
+        num_levels = len(self.affectations[j + self.d_cont].classes_) - 1
         if self.m_start > num_levels:
             m_qual[j] = num_levels
 
diff --git a/glmdisc/_generateData.py b/glmdisc/_generateData.py
index 06bc84e..619e680 100644
--- a/glmdisc/_generateData.py
+++ b/glmdisc/_generateData.py
@@ -36,19 +36,20 @@ def generate_data(n, d, theta=None, plot=False):
     for i in range(d):
         xd[:, i] = pd.cut(x[:, i], bins=cuts, labels=[0, 1, 2])
 
-    log_odd = np.array([0] * n)
+    log_odd = np.zeros(n)
     for i in range(n):
         for j in range(d):
             log_odd[i] += theta[int(xd[i, j]), j]
 
     p = 1 / (1 + np.exp(- log_odd))
-    y = np.random.binomial(1, p)
+    y = np.random.binomial(1, p, size=n)
 
     if plot and d > 1:
-        for unique_pair_of_feature_as_tuple in combinations(range(d), 2):
-            plt.scatter(x[:, unique_pair_of_feature_as_tuple[0]],
-                        x[:, unique_pair_of_feature_as_tuple[1]],
-                        c=y,
-                        s=[1 for _ in range(len(x))])
+        for j1, j2 in combinations(range(d), 2):
+            plt.figure()
+            plt.scatter(x[:, j1], x[:, j2], c=y, s=1)
+            plt.xlabel("x" + str(j1))
+            plt.ylabel("x" + str(j2))
+            plt.show(block=False)
 
     return [x, y, theta]
diff --git a/glmdisc/_gini_utils.py b/glmdisc/_gini_utils.py
index 78bd017..f1de898 100644
--- a/glmdisc/_gini_utils.py
+++ b/glmdisc/_gini_utils.py
@@ -24,7 +24,7 @@ def compute_midrank(x):
     J = np.argsort(x)
     Z = x[J]
     N = len(x)
-    T = np.zeros(N, dtype=np.float)
+    T = np.zeros(N, dtype=float)
     i = 0
     while i < N:
         j = i
@@ -32,7 +32,7 @@ def compute_midrank(x):
             j += 1
         T[i:j] = 0.5 * (i + j - 1)
         i = j
-    T2 = np.empty(N, dtype=np.float)
+    T2 = np.empty(N, dtype=float)
     # Note(kazeevn) +1 is due to Python using 0-based indexing
     # instead of 1-based in the AUC formula in the paper
     T2[J] = T + 1
@@ -50,7 +50,7 @@ def compute_midrank_weight(x, sample_weight):
     Z = x[J]
     cumulative_weight = np.cumsum(sample_weight[J])
     N = len(x)
-    T = np.zeros(N, dtype=np.float)
+    T = np.zeros(N, dtype=float)
     i = 0
     while i < N:
         j = i
@@ -58,7 +58,7 @@ def compute_midrank_weight(x, sample_weight):
             j += 1
         T[i:j] = cumulative_weight[i:j].mean()
         i = j
-    T2 = np.empty(N, dtype=np.float)
+    T2 = np.empty(N, dtype=float)
     T2[J] = T
     return T2
 
@@ -99,9 +99,9 @@ def fast_de_long_weights(predictions_sorted_transposed, label_1_count, sample_we
     negative_examples = predictions_sorted_transposed[:, m:]
     k = predictions_sorted_transposed.shape[0]
 
-    tx = np.empty([k, m], dtype=np.float)
-    ty = np.empty([k, n], dtype=np.float)
-    tz = np.empty([k, m + n], dtype=np.float)
+    tx = np.empty([k, m], dtype=float)
+    ty = np.empty([k, n], dtype=float)
+    tz = np.empty([k, m + n], dtype=float)
     for r in range(k):
         tx[r, :] = compute_midrank_weight(positive_examples[r, :], sample_weight[:m])
         ty[r, :] = compute_midrank_weight(negative_examples[r, :], sample_weight[m:])
@@ -149,9 +149,9 @@ def fast_de_long_no_weights(predictions_sorted_transposed, label_1_count):
     negative_examples = predictions_sorted_transposed[:, m:]
     k = predictions_sorted_transposed.shape[0]
 
-    tx = np.empty([k, m], dtype=np.float)
-    ty = np.empty([k, n], dtype=np.float)
-    tz = np.empty([k, m + n], dtype=np.float)
+    tx = np.empty([k, m], dtype=float)
+    ty = np.empty([k, n], dtype=float)
+    tz = np.empty([k, m + n], dtype=float)
     for r in range(k):
         tx[r, :] = compute_midrank(positive_examples[r, :])
         ty[r, :] = compute_midrank(negative_examples[r, :])
@@ -201,7 +201,7 @@ def gini(ground_truth, predictions, sample_weight=None, level=.95):
        sample_weight: np.array of floats of the sample weights
        level: float (between 0 and 1) describing the level required for the confidence interval
     """
-    if (not isinstance(level, float)) or (not level < 1 and level > 0):
+    if not isinstance(level, float) or not 0 < level < 1:
         raise ValueError("level should be a single float between 0 and 1")
     auc, auc_cov = delong_roc_variance(ground_truth, predictions, sample_weight)
     auc_std = np.sqrt(auc_cov)
diff --git a/glmdisc/_plot.py b/glmdisc/_plot.py
index 39e4b4d..f2a939d 100644
--- a/glmdisc/_plot.py
+++ b/glmdisc/_plot.py
@@ -50,24 +50,31 @@ def _plot_disc(self, predictors_cont_number, predictors_qual_number):
 
 
 def _plot_logodd(self):
-    # Gérer les manquants dans le GAM
-    lignes_completes = np.invert(
-        np.isnan(self.predictors_cont).sum(axis=1).astype(bool))
+    # Handle missing values in the GAM: keep the complete rows of the continuous features
+    if self.d_cont > 0:
+        lignes_completes = np.invert(np.isnan(self.predictors_cont).any(axis=1))
+    else:
+        lignes_completes = np.ones(self.n, dtype=bool)
+
+    features = []
+    if self.d_cont > 0:
+        features.append(pd.DataFrame(self.predictors_cont[lignes_completes, :]).apply(
+            lambda x: x.astype('float')))
+    if self.d_qual > 0:
+        features.append(pd.DataFrame(self.predictors_qual[lignes_completes, :]).apply(
+            lambda x: x.astype('category')))
 
-    # Fit du GAM sur tout le monde
+    # Fit the GAM on everyone
     gam = LogisticGAM(dtype=['numerical' for _ in range(self.d_cont)] + ['categorical' for _ in range(
-        self.d_qual)]).fit(
-        pd.concat([pd.DataFrame(self.predictors_cont[lignes_completes, :]).apply(
-            lambda x: x.astype('float')),
-            pd.DataFrame(self.predictors_qual[lignes_completes, :]).apply(
-                lambda x: x.astype('category'))], axis=1), self.labels[lignes_completes])
-
-    # Quelles que soient les valeurs de predictors_cont_number et
-    # predictors_qual_number, on plot tout pour l'instant
+        self.d_qual)]).fit(pd.concat(features, axis=1), self.labels[lignes_completes])
+
+    # Whatever the values of predictors_cont_number and
+    # predictors_qual_number, everything is plotted for now
     plt.figure()
-    fig, axs = plt.subplots(1, self.d_cont + self.d_qual)
+    # squeeze=False keeps axs a 2D array even with a single feature
+    fig, axs = plt.subplots(1, self.d_cont + self.d_qual, squeeze=False)
     plt.rcParams['figure.figsize'] = (28, 8)
-    for i, ax in enumerate(axs):
+    for i, ax in enumerate(axs[0]):
         try:
             XX = gam.generate_X_grid(term=i)
             ax.plot(XX[:, i], gam.partial_dependence(term=i, X=XX))
diff --git a/tests/test_gini.py b/tests/test_gini.py
new file mode 100644
index 0000000..64dbd80
--- /dev/null
+++ b/tests/test_gini.py
@@ -0,0 +1,50 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+import pytest
+from sklearn.metrics import roc_auc_score
+
+import glmdisc
+
+
+def _scores(n=200, seed=0):
+    rng = np.random.RandomState(seed)
+    y = rng.binomial(1, 0.4, size=n)
+    p = np.clip(0.3 * y + rng.uniform(size=n) * 0.7, 0, 1)
+    return y, p
+
+
+def test_gini_matches_auc():
+    y, p = _scores()
+    value, ci = glmdisc.gini(y, p)
+    assert value == pytest.approx(2 * roc_auc_score(y, p) - 1)
+    assert ci[0] < value < ci[1] <= 1
+
+
+def test_gini_with_ties():
+    y, p = _scores()
+    p = np.round(p, 1)
+    value, _ = glmdisc.gini(y, p)
+    assert value == pytest.approx(2 * roc_auc_score(y, p) - 1)
+
+
+def test_gini_with_sample_weight():
+    y, p = _scores()
+    w = np.random.RandomState(1).uniform(0.5, 2, size=len(y))
+    value, _ = glmdisc.gini(y, p, sample_weight=w)
+    assert value == pytest.approx(2 * roc_auc_score(y, p, sample_weight=w) - 1)
+
+
+@pytest.mark.parametrize("level", [0., -.5, 1., 1.5, 1, "0.95", None])
+def test_gini_invalid_level(level):
+    y, p = _scores()
+    with pytest.raises(ValueError):
+        glmdisc.gini(y, p, level=level)
+
+
+def test_gini_level_widens_interval():
+    y, p = _scores()
+    _, ci_90 = glmdisc.gini(y, p, level=.9)
+    _, ci_99 = glmdisc.gini(y, p, level=.99)
+    assert ci_99[0] < ci_90[0]
+    assert ci_99[1] >= ci_90[1]
diff --git a/tests/test_init.py b/tests/test_init.py
index ff574ab..90abe26 100644
--- a/tests/test_init.py
+++ b/tests/test_init.py
@@ -64,6 +64,17 @@ def test_m_start():
     with pytest.raises(ValueError):
         glmdisc.Glmdisc(m_start=1)
 
+    with pytest.raises(ValueError):
+        glmdisc.Glmdisc(m_start=2.5)
+
+
+def test_burn_in():
+    assert glmdisc.Glmdisc(burn_in=0).burn_in == 0
+
+    for burn_in in [-1, 2.5, "5", None, True]:
+        with pytest.raises(ValueError):
+            glmdisc.Glmdisc(burn_in=burn_in)
+
 
 def test_gini_penalized(caplog):
     glmdisc.Glmdisc(validation=False,
@@ -83,3 +94,11 @@ def test_validation_criterion(caplog):
 
     assert caplog.records[0].message == ('No need to penalize the log-likelihood when a validation set is used. '
                                          'Using log-likelihood instead.')
+
+
+def test_not_fitted_sem():
+    model = glmdisc.Glmdisc()
+    with pytest.raises(glmdisc.NotFittedError):
+        model._check_is_fitted()
+    with pytest.raises(glmdisc.NotFittedError):
+        model.predict(np.zeros((3, 1)), None)
diff --git a/tests/tests_NN/test_NN_fit.py b/tests/tests_NN/test_NN_fit.py
index 1d8523a..113154a 100644
--- a/tests/tests_NN/test_NN_fit.py
+++ b/tests/tests_NN/test_NN_fit.py
@@ -5,6 +5,8 @@
 import pandas as pd
 import random
 import glmdisc
+import sklearn as sk
+from glmdisc._fitNN import _initialize_neural_net
 import tensorflow
 from tensorflow.keras.optimizers import Adagrad
 from tensorflow.keras.callbacks import EarlyStopping
@@ -209,3 +211,20 @@ def test_fit():
     model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=20)
     model = glmdisc.Glmdisc(algorithm="NN", criterion="gini", test=False, validation=False)
     model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=20)
+
+
+def test_initialize_neural_net_reuses_label_encoders(monkeypatch):
+    np.random.seed(0)
+    x, y, theta = glmdisc.Glmdisc.generate_data(100, 2)
+    xd = np.digitize(x, [0.333, 0.666]).astype(float)
+    model = glmdisc.Glmdisc(algorithm="NN", validation=False, test=False)
+    model.predictors_cont = x
+    model.predictors_qual = xd
+    model.labels = y
+    model._init_disc(model._calculate_shape())
+
+    def fail_fit(self, y):
+        raise AssertionError("LabelEncoder fitted again")
+
+    monkeypatch.setattr(sk.preprocessing.LabelEncoder, "fit", fail_fit)
+    _initialize_neural_net(model)
diff --git a/tests/tests_SEM/test_SEM_best_formula.py b/tests/tests_SEM/test_SEM_best_formula.py
index 52e1508..46b02f9 100644
--- a/tests/tests_SEM/test_SEM_best_formula.py
+++ b/tests/tests_SEM/test_SEM_best_formula.py
@@ -3,6 +3,7 @@
 import numpy as np
 import pandas as pd
 import glmdisc
+from glmdisc._bestFormula import _best_formula_cont_sem
 
 
 def test_best_formula(caplog):
@@ -45,3 +46,27 @@ def test_best_formula_no_cont():
     model = glmdisc.Glmdisc(validation=False, test=False)
     model.fit(predictors_cont=None, predictors_qual=xd, labels=y, iter=11)
     model.best_formula()
+
+
+class _Stub:
+    d_cont = 2
+
+
+def test_best_formula_cont_sem_missing_values():
+    stub = _Stub()
+    stub.predictors_cont = np.array([[1., np.nan],
+                                     [2., np.nan],
+                                     [3., np.nan],
+                                     [4., np.nan]])
+    emap_best = np.array([[0, 0], [0, 0], [1, 0], [1, 0]])
+    assert _best_formula_cont_sem(stub, emap_best) == [[3.], []]
+
+
+def test_best_formula_cont_sem_missing_level():
+    stub = _Stub()
+    stub.predictors_cont = np.array([[1., 5.],
+                                     [2., 6.],
+                                     [3., np.nan],
+                                     [np.nan, 7.]])
+    emap_best = np.array([[0, 0], [0, 0], [1, 1], [2, 1]])
+    assert _best_formula_cont_sem(stub, emap_best) == [[3.], [7.]]
diff --git a/tests/tests_SEM/test_SEM_discretize.py b/tests/tests_SEM/test_SEM_discretize.py
index 0fac42a..a941a98 100644
--- a/tests/tests_SEM/test_SEM_discretize.py
+++ b/tests/tests_SEM/test_SEM_discretize.py
@@ -5,6 +5,7 @@
 import pytest
 import random
 import glmdisc
+from glmdisc._discretize import _check_args_discretize_nn
 
 
 def test_discretize_new():
@@ -67,3 +68,34 @@ def test_discretize_wrong():
         model.discretize(predictors_cont=None, predictors_qual=xd)
     with pytest.raises(ValueError):
         model.discretize(predictors_cont=x, predictors_qual=xd)
+
+
+def test_discretize_wrong_messages(caplog):
+    n = 100
+    d = 2
+    x, y, theta = glmdisc.Glmdisc.generate_data(n, d)
+    model = glmdisc.Glmdisc(validation=False, test=False)
+    random.seed(1)
+    np.random.seed(1)
+    model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11)
+    with pytest.raises(ValueError):
+        model.discretize(predictors_cont=x, predictors_qual=x.astype(str))
+    assert "Shape of 2 for predictors_qual" in caplog.records[-1].message
+    with pytest.raises(ValueError):
+        model.discretize(predictors_cont=x[:, :1], predictors_qual=None)
+    assert "Shape of 1 for predictors_cont" in caplog.records[-1].message
+
+
+class _Stub:
+    d_cont = 1
+    d_qual = 1
+
+
+def test_check_args_discretize_nn_messages(caplog):
+    x = np.zeros((5, 1))
+    with pytest.raises(ValueError):
+        _check_args_discretize_nn(_Stub(), x, np.zeros((5, 3)))
+    assert "Shape of 3 for predictors_qual" in caplog.records[-1].message
+    with pytest.raises(ValueError):
+        _check_args_discretize_nn(_Stub(), np.zeros((5, 2)), np.zeros((5, 1)))
+    assert "Shape of 2 for predictors_cont" in caplog.records[-1].message
diff --git a/tests/tests_SEM/test_SEM_generate.py b/tests/tests_SEM/test_SEM_generate.py
index c270c16..da75a35 100644
--- a/tests/tests_SEM/test_SEM_generate.py
+++ b/tests/tests_SEM/test_SEM_generate.py
@@ -1,5 +1,6 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
+import matplotlib.pyplot as plt
 import pytest
 import numpy as np
 import glmdisc
@@ -34,3 +35,44 @@ def test_generate_theta_error():
 
     with pytest.raises(ValueError):
         glmdisc.Glmdisc.generate_data(n=800, d=d, theta=theta)
+
+
+def test_generate_log_odd_not_truncated(monkeypatch):
+    theta = np.array([[0.4], [0.9], [-0.6]])
+    drawn = {}
+
+    def fake_binomial(n, p, size=None):
+        drawn["p"] = p
+        return np.zeros(len(p), dtype=int)
+
+    monkeypatch.setattr(np.random, "binomial", fake_binomial)
+    np.random.seed(0)
+    x, y, theta2 = glmdisc.Glmdisc.generate_data(n=50, d=1, theta=theta)
+
+    levels = np.digitize(x[:, 0], [0.333, 0.666])
+    expected = 1 / (1 + np.exp(-theta[levels, 0]))
+    np.testing.assert_allclose(drawn["p"], expected)
+
+
+def test_generate_labels_size(monkeypatch):
+    sizes = []
+    binomial = np.random.binomial
+
+    def spy_binomial(n, p, size=None):
+        sizes.append(size)
+        return binomial(n, p, size=size)
+
+    monkeypatch.setattr(np.random, "binomial", spy_binomial)
+    x, y, theta = glmdisc.Glmdisc.generate_data(n=1, d=3)
+    assert sizes == [1]
+    assert y.shape == (1,)
+
+
+def test_generate_plot_one_figure_per_pair():
+    plt.close('all')
+    glmdisc.Glmdisc.generate_data(n=50, d=3, plot=True)
+    assert len(plt.get_fignums()) == 3
+    for num in plt.get_fignums():
+        ax = plt.figure(num).axes[0]
+        assert len(ax.collections) == 1
+    plt.close('all')
diff --git a/tests/tests_SEM/test_SEM_init_disc.py b/tests/tests_SEM/test_SEM_init_disc.py
new file mode 100644
index 0000000..27a141a
--- /dev/null
+++ b/tests/tests_SEM/test_SEM_init_disc.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import numpy as np
+from sklearn.preprocessing import LabelEncoder
+
+import glmdisc
+
+
+def _model(n=300, levels=("a", "b", "c"), m_start=20):
+    model = glmdisc.Glmdisc(m_start=m_start)
+    model.predictors_cont = np.random.uniform(size=(n, 2))
+    model.predictors_qual = np.random.choice(np.array(levels), size=(n, 2)).astype(object)
+    model.labels = np.random.binomial(1, .5, size=n)
+    return model
+
+
+def test_init_disc_fits_each_label_encoder_once(monkeypatch):
+    np.random.seed(0)
+    model = _model()
+    fits = []
+    original_fit = LabelEncoder.fit
+
+    def counting_fit(self, y):
+        fits.append(len(y))
+        return original_fit(self, y)
+
+    monkeypatch.setattr(LabelEncoder, "fit", counting_fit)
+    complete = model._calculate_shape()
+    model._init_disc(complete)
+    assert len(fits) == model.d_qual
+
+
+def test_init_disc_draws_every_level():
+    np.random.seed(0)
+    model = _model(levels=("a", "b", "c"), m_start=20)
+    edisc, predictors_trans = model._init_disc(model._calculate_shape())
+    for j in range(model.d_qual):
+        assert set(np.unique(edisc[:, model.d_cont + j])) == {0, 1, 2}
+
+
+def test_init_disc_caps_groups_at_m_start():
+    np.random.seed(0)
+    model = _model(levels=tuple("abcdefgh"), m_start=3)
+    edisc, predictors_trans = model._init_disc(model._calculate_shape())
+    for j in range(model.d_qual):
+        assert set(np.unique(edisc[:, model.d_cont + j])) == {0, 1, 2}
+
+
+def test_init_disc_single_level():
+    np.random.seed(0)
+    model = _model(levels=("a",), m_start=20)
+    edisc, predictors_trans = model._init_disc(model._calculate_shape())
+    np.testing.assert_array_equal(edisc[:, 2:], 0)
diff --git a/tests/tests_SEM/test_SEM_plot.py b/tests/tests_SEM/test_SEM_plot.py
index 564eece..59b9f68 100644
--- a/tests/tests_SEM/test_SEM_plot.py
+++ b/tests/tests_SEM/test_SEM_plot.py
@@ -56,3 +56,39 @@ def test_plot_log(caplog):
     assert caplog.records[-1].message == ("A single int (more than 0 and less than the "
                                           "number of columns in predictors_cont) must be "
                                           "provided for predictors_cont_number")
+
+
+def _digitized(x):
+    return np.digitize(x, [0.333, 0.666]).astype(float)
+
+
+def test_plot_logodd_cont_only():
+    np.random.seed(1)
+    x, y, theta = glmdisc.Glmdisc.generate_data(300, 2)
+    model = glmdisc.Glmdisc(validation=False, test=False)
+    model.fit(predictors_cont=x, predictors_qual=None, labels=y, iter=11)
+    model.plot(plot_type="logodd")
+    plt.close('all')
+
+
+def test_plot_logodd_qual_only():
+    np.random.seed(1)
+    x, y, theta = glmdisc.Glmdisc.generate_data(300, 2)
+    model = glmdisc.Glmdisc(validation=False, test=False)
+    model.fit(predictors_cont=None, predictors_qual=_digitized(x), labels=y, iter=11)
+    model.plot(plot_type="logodd")
+    plt.close('all')
+
+
+def test_plot_logodd_single_feature():
+    np.random.seed(1)
+    x, y, theta = glmdisc.Glmdisc.generate_data(300, 2)
+    model = glmdisc.Glmdisc(validation=False, test=False)
+    model.fit(predictors_cont=x[:, :1], predictors_qual=None, labels=y, iter=11)
+    model.plot(plot_type="logodd")
+    plt.close('all')
+
+    model = glmdisc.Glmdisc(validation=False, test=False)
+    model.fit(predictors_cont=None, predictors_qual=_digitized(x[:, :1]), labels=y, iter=11)
+    model.plot(plot_type="logodd")
+    plt.close('all')