Repository: TRR318/scikit-psl (language: python)

Issue:
Extend PSL to be capable of working with continuos datasets
in v0.2.0 the psl can only work with binary datasets. There is a binarization class to allow for continuos datasets, but doing the binarization in a preprocessing step is probably not optimal.

future extensions could be:
- allow for categorical inputs (producing set-style decision rules)
- allow multiclass prediction
Extend PSL to be capable of working with continuos datasets
in v0.2.0 the psl can only work with binary datasets. There is a binarization class to allow for continuos datasets, but doing the binarization in a preprocessing step is probably not optimal.

future extensions could be:
- allow for categorical inputs (producing set-style decision rules)
- allow multiclass prediction


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3ab51e1..2e8f66f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,17 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## 0.3.0 - 2023-08-10
+
+### Added
+
+- PSL classifier can now run with continuous data and optimally (wrt. expected entropy) select thresholds to binarize the data
+
+### Changed
+
+- Significantly improved optimum calculation for MinEntropyBinarizer (the same optimization algorithm is shared with the psls internal binarization algorithm)
+
+
 ## 0.2.0 - 2023-08-10
 
 ### Added
diff --git a/README.md b/README.md
index eb7ede1..d7ed14f 100644
--- a/README.md
+++ b/README.md
@@ -20,30 +20,35 @@ pip install scikit-psl
 
 # Usage
 
+For examples have a look at the `examples` folder, but here is a simple example
+
+
 ```python
 from sklearn.datasets import make_classification
-from sklearn.model_selection import ShuffleSplit
+from sklearn.model_selection import train_test_split
 
 from skpsl import ProbabilisticScoringList
 
 # Generating synthetic data with continuous features and a binary target variable
-X, y = make_classification(random_state=42)
-X = (X > .5).astype(int)
-
-psl = ProbabilisticScoringList([-1, 1, 2])
-
-for train, test in ShuffleSplit(1, test_size=.2, random_state=42).split(X):
-    psl.fit(X[train], y[train])
-    print(f"Brier score: {psl.score(X[test], y[test]):.4f}")
-    #>  Brier score: 0.1924  (lower is better)
-
-    df = psl.inspect(5)
-    print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))    
-    #>  Stage  Score  T = -3  T = -2  T = -1  T = 0  T = 1  T = 2  T = 3
-    #>   0        -       -       -       -   0.54      -      -      - 
-    #>   1     2.00       -       -       -   0.18      -   0.97      - 
-    #>   2    -1.00       -       -    0.00   0.28   0.91   1.00      - 
-    #>   3    -1.00       -    0.00    0.07   0.86   0.91   1.00      - 
-    #>   4     1.00       -    0.00    0.00   0.29   0.92   1.00   1.00 
-    #>   5    -1.00    0.00    0.00    0.00   0.40   1.00   1.00   1.00
-```
+X, y = make_classification(n_informative=10, random_state=42)
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
+
+psl = ProbabilisticScoringList({-1, 1, 2})
+psl.fit(X_train, y_train)
+print(f"Brier score: {psl.score(X_test, y_test):.4f}")
+"""
+Brier score: 0.2438  (lower is better)
+"""
+
+df = psl.inspect(5)
+print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))
+"""
+ Stage Threshold  Score  T = -2  T = -1  T = 0  T = 1  T = 2  T = 3  T = 4  T = 5
+  0            -     -       -       -   0.51      -      -      -      -      - 
+  1     >-2.4245  2.00       -       -   0.00      -   0.63      -      -      - 
+  2     >-0.9625 -1.00       -    0.00   0.00   0.48   1.00      -      -      - 
+  3      >0.4368 -1.00    0.00    0.00   0.12   0.79   1.00      -      -      - 
+  4     >-0.9133  1.00    0.00    0.00   0.12   0.12   0.93   1.00      -      - 
+  5      >2.4648  2.00    0.00    0.00   0.07   0.07   0.92   1.00   1.00   1.00 
+"""
+```
\ No newline at end of file
diff --git a/examples/fitting psl on partially continuous data.py b/examples/fitting psl on partially continuous data.py
new file mode 100644
index 0000000..f1d9d27
--- /dev/null
+++ b/examples/fitting psl on partially continuous data.py
@@ -0,0 +1,18 @@
+from sklearn.datasets import make_classification
+from sklearn.model_selection import train_test_split
+
+from skpsl import ProbabilisticScoringList
+
+if __name__ == '__main__':
+    X, y = make_classification(n_informative=10, random_state=42)
+    # make
+    X[:, [9, 10]] = X[:, [9, 10]] > 0
+    X[:, [2, 6, 7]] = X[:, [2, 6, 7]] > -5
+    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
+
+    psl = ProbabilisticScoringList({-1, 1, 2})
+    psl.fit(X_train, y_train)
+    print(f"Brier score: {psl.score(X_test, y_test):.4f}")
+
+    df = psl.inspect(6)
+    print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))
diff --git a/examples/simple.py b/examples/simple.py
new file mode 100644
index 0000000..3814356
--- /dev/null
+++ b/examples/simple.py
@@ -0,0 +1,27 @@
+from sklearn.datasets import make_classification
+from sklearn.model_selection import train_test_split
+
+from skpsl import ProbabilisticScoringList
+
+if __name__ == '__main__':
+    X, y = make_classification(n_informative=10, random_state=42)
+    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
+
+    psl = ProbabilisticScoringList({-1, 1, 2})
+    psl.fit(X_train, y_train)
+    print(f"Brier score: {psl.score(X_test, y_test):.4f}")
+    """
+    Brier score: 0.2438  (lower is better)
+    """
+
+    df = psl.inspect(5)
+    print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))
+    """
+     Stage Threshold  Score  T = -2  T = -1  T = 0  T = 1  T = 2  T = 3  T = 4  T = 5
+      0            -     -       -       -   0.51      -      -      -      -      - 
+      1     >-2.4245  2.00       -       -   0.00      -   0.63      -      -      - 
+      2     >-0.9625 -1.00       -    0.00   0.00   0.48   1.00      -      -      - 
+      3      >0.4368 -1.00    0.00    0.00   0.12   0.79   1.00      -      -      - 
+      4     >-0.9133  1.00    0.00    0.00   0.12   0.12   0.93   1.00      -      - 
+      5      >2.4648  2.00    0.00    0.00   0.07   0.07   0.92   1.00   1.00   1.00 
+    """
diff --git a/examples/using binarizer pipeline.py b/examples/using binarizer pipeline.py
new file mode 100644
index 0000000..6c207ac
--- /dev/null
+++ b/examples/using binarizer pipeline.py
@@ -0,0 +1,25 @@
+from sklearn.datasets import make_classification
+from sklearn.model_selection import train_test_split
+from sklearn.pipeline import Pipeline
+
+from skpsl import MinEntropyBinarizer, ProbabilisticScoringList
+
+if __name__ == '__main__':
+    X, y = make_classification(n_informative=10, random_state=42)
+    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
+
+    # The MinEntropyBinarizer calculates an optimal threshold beforehand.
+    # It will calculate a threshold for binarizing the continuous data
+    # The optimization criterion is an entropy based impurity measure, similar to the ones used for decision trees.
+    # The PSL than operates on the binary data
+    pipe = Pipeline([("binarizer", MinEntropyBinarizer()),
+                     ("psl", ProbabilisticScoringList({-1, 1, 2}))])
+    pipe.fit(X_train, y_train)
+    print(f"Brier score: {pipe.score(X_test, y_test):.4f}")
+    # >  Brier score: 0.2184  (lower is better)
+
+    df = pipe["binarizer"].inspect()
+    print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))
+
+    df = pipe["psl"].inspect(5)
+    print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))
diff --git a/pyproject.toml b/pyproject.toml
index 0731911..67d72aa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "scikit-psl"
-version = "0.2.0"
+version = "0.3.0"
 description = "Probabilistic Scoring List classifier"
 license = "MIT"
 authors = ["Stefan Heid <stefan.heid@upb.de>", "Jonas Hanselle <jonas.hanselle@ifi.lmu.de>"]
@@ -12,20 +12,16 @@ packages = [{ include = "skpsl" }]
 python = ">=3.9,<3.13"
 scikit-learn = "^1.3.0"
 numpy = "^1.25.2"
-scipy = "^1.11.1"
+scipy = "^1.11.2"
 joblib = "^1.3.2"
-pandas = "^2.0.3"
+pandas = "^2.1.0"
+sortedcontainers = "^2.4.0"
 
-sphinx = { version = "^7.1", optional = true }
-sphinx_rtd_theme = { version = "^1.2", optional = true  }
-
-pytest = { version = "^7.4", optional = true }
-
-
-[tool.poetry.extras]
-doc = ["sphinx", "sphinx_rtd_theme"]
-test = ["pytest"]
+[tool.poetry.group.dev.dependencies]
+sphinx = "^7.2.5"
+sphinx_rtd_theme = "^1.3"
 
+pytest = "^7.4"
 
 [build-system]
 requires = ["poetry-core"]
diff --git a/scratch/psl_describe.ipynb b/scratch/psl_describe.ipynb
index 2e98e3a..06a8ec8 100644
--- a/scratch/psl_describe.ipynb
+++ b/scratch/psl_describe.ipynb
@@ -18,7 +18,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 6,
    "id": "764414e5-261c-4a77-b14b-2235472b9baf",
    "metadata": {},
    "outputs": [
@@ -26,13 +26,14 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Brier score: 0.1924\n"
+      "Brier score: 0.6500\n"
      ]
     }
    ],
    "source": [
     "# Generating synthetic data with continuous features and a binary target variable\n",
     "X, y = make_classification(random_state=42)\n",
+    "#X = (X*10).round().astype(int)\n",
     "X = (X > .5).astype(int)\n",
     "\n",
     "for train, test in ShuffleSplit(1, test_size=.2, random_state=42).split(X):\n",
@@ -43,7 +44,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 7,
    "id": "26a044d2-d3a1-43eb-bbef-b8e7050ce568",
    "metadata": {},
    "outputs": [
@@ -70,13 +71,13 @@
        "      <th></th>\n",
        "      <th>Stage</th>\n",
        "      <th>Score</th>\n",
-       "      <th>T = -3</th>\n",
-       "      <th>T = -2</th>\n",
-       "      <th>T = -1</th>\n",
-       "      <th>T = 0</th>\n",
-       "      <th>T = 1</th>\n",
-       "      <th>T = 2</th>\n",
-       "      <th>T = 3</th>\n",
+       "      <th>T = -4.0</th>\n",
+       "      <th>T = -3.1666666666666665</th>\n",
+       "      <th>T = -2.333333333333333</th>\n",
+       "      <th>T = -1.5</th>\n",
+       "      <th>T = -0.6666666666666665</th>\n",
+       "      <th>T = 0.16666666666666696</th>\n",
+       "      <th>T = 1.0</th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
@@ -84,89 +85,242 @@
        "      <th>0</th>\n",
        "      <td>0</td>\n",
        "      <td>-</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
        "      <td>0.5375</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
        "      <td>1</td>\n",
-       "      <td>2.0</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
+       "      <td>1.0</td>\n",
        "      <td>0.1818</td>\n",
-       "      <td>-</td>\n",
+       "      <td>0.1818</td>\n",
+       "      <td>0.1818</td>\n",
+       "      <td>0.1818</td>\n",
+       "      <td>0.1818</td>\n",
+       "      <td>0.3136</td>\n",
        "      <td>0.9722</td>\n",
-       "      <td>-</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
        "      <td>2</td>\n",
        "      <td>-1.0</td>\n",
-       "      <td>-</td>\n",
-       "      <td>-</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.2759</td>\n",
-       "      <td>0.9091</td>\n",
-       "      <td>1.0</td>\n",
-       "      <td>-</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.1293</td>\n",
+       "      <td>0.4898</td>\n",
+       "      <td>1.0000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>3</th>\n",
        "      <td>3</td>\n",
        "      <td>-1.0</td>\n",
-       "      <td>-</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.069</td>\n",
-       "      <td>0.8571</td>\n",
-       "      <td>0.9091</td>\n",
-       "      <td>1.0</td>\n",
-       "      <td>-</td>\n",
+       "      <td>0.2692</td>\n",
+       "      <td>0.2692</td>\n",
+       "      <td>0.2692</td>\n",
+       "      <td>0.2692</td>\n",
+       "      <td>0.3590</td>\n",
+       "      <td>0.6154</td>\n",
+       "      <td>1.0000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>4</th>\n",
        "      <td>4</td>\n",
-       "      <td>1.0</td>\n",
+       "      <td>-1.0</td>\n",
+       "      <td>0.3333</td>\n",
+       "      <td>0.3333</td>\n",
+       "      <td>0.3333</td>\n",
+       "      <td>0.3333</td>\n",
+       "      <td>0.4373</td>\n",
+       "      <td>0.7043</td>\n",
+       "      <td>1.0000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>5</th>\n",
+       "      <td>5</td>\n",
+       "      <td>-1.0</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.1667</td>\n",
+       "      <td>0.3173</td>\n",
+       "      <td>0.5722</td>\n",
+       "      <td>0.9561</td>\n",
+       "      <td>1.0000</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "   Stage Score  T = -4.0  T = -3.1666666666666665  T = -2.333333333333333  \\\n",
+       "0      0     -    0.5375                   0.5375                  0.5375   \n",
+       "1      1   1.0    0.1818                   0.1818                  0.1818   \n",
+       "2      2  -1.0    0.0000                   0.0000                  0.0000   \n",
+       "3      3  -1.0    0.2692                   0.2692                  0.2692   \n",
+       "4      4  -1.0    0.3333                   0.3333                  0.3333   \n",
+       "5      5  -1.0    0.0000                   0.0000                  0.1667   \n",
+       "\n",
+       "   T = -1.5  T = -0.6666666666666665  T = 0.16666666666666696  T = 1.0  \n",
+       "0    0.5375                   0.5375                   0.5375   0.5375  \n",
+       "1    0.1818                   0.1818                   0.3136   0.9722  \n",
+       "2    0.0000                   0.1293                   0.4898   1.0000  \n",
+       "3    0.2692                   0.3590                   0.6154   1.0000  \n",
+       "4    0.3333                   0.4373                   0.7043   1.0000  \n",
+       "5    0.3173                   0.5722                   0.9561   1.0000  "
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "df = psl.inspect(5)\n",
+    "\n",
+    "pd.set_option(\"display.precision\", 4)\n",
+    "df.fillna(\"-\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "aa327ca6-0ec9-4ca2-8c9e-3d59a26701e4",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>Stage</th>\n",
+       "      <th>Score</th>\n",
+       "      <th>T = -4.0</th>\n",
+       "      <th>T = -3.0</th>\n",
+       "      <th>T = -2.0</th>\n",
+       "      <th>T = -1.0</th>\n",
+       "      <th>T = 0.0</th>\n",
+       "      <th>T = 1.0</th>\n",
+       "      <th>T = 2.0</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>0</td>\n",
        "      <td>-</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.2857</td>\n",
-       "      <td>0.9167</td>\n",
-       "      <td>1.0</td>\n",
-       "      <td>1.0</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "      <td>0.5375</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>1</td>\n",
+       "      <td>-1.0</td>\n",
+       "      <td>0.3667</td>\n",
+       "      <td>0.3667</td>\n",
+       "      <td>0.3667</td>\n",
+       "      <td>0.3667</td>\n",
+       "      <td>0.4783</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>2</td>\n",
+       "      <td>-1.0</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.3065</td>\n",
+       "      <td>0.7143</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>3</td>\n",
+       "      <td>-1.0</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.1364</td>\n",
+       "      <td>0.6111</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>4</td>\n",
+       "      <td>2.0</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0400</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>5</th>\n",
        "      <td>5</td>\n",
        "      <td>-1.0</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.0</td>\n",
-       "      <td>0.4000</td>\n",
-       "      <td>1.0</td>\n",
-       "      <td>1.0</td>\n",
-       "      <td>1.0</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>0.0000</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
+       "      <td>1.0000</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "text/plain": [
-       "   Stage Score T = -3 T = -2 T = -1   T = 0   T = 1   T = 2 T = 3\n",
-       "0      0     -      -      -      -  0.5375       -       -     -\n",
-       "1      1   2.0      -      -      -  0.1818       -  0.9722     -\n",
-       "2      2  -1.0      -      -    0.0  0.2759  0.9091     1.0     -\n",
-       "3      3  -1.0      -    0.0  0.069  0.8571  0.9091     1.0     -\n",
-       "4      4   1.0      -    0.0    0.0  0.2857  0.9167     1.0   1.0\n",
-       "5      5  -1.0    0.0    0.0    0.0  0.4000     1.0     1.0   1.0"
+       "   Stage Score  T = -4.0  T = -3.0  T = -2.0  T = -1.0  T = 0.0  T = 1.0  \\\n",
+       "0      0     -    0.5375    0.5375    0.5375    0.5375   0.5375   0.5375   \n",
+       "1      1  -1.0    0.3667    0.3667    0.3667    0.3667   0.4783   1.0000   \n",
+       "2      2  -1.0    0.0000    0.0000    0.0000    0.3065   0.7143   1.0000   \n",
+       "3      3  -1.0    0.0000    0.0000    0.0000    0.1364   0.6111   1.0000   \n",
+       "4      4   2.0    0.0000    0.0000    0.0000    0.0400   1.0000   1.0000   \n",
+       "5      5  -1.0    0.0000    0.0000    0.0000    0.0000   1.0000   1.0000   \n",
+       "\n",
+       "   T = 2.0  \n",
+       "0   0.5375  \n",
+       "1   1.0000  \n",
+       "2   1.0000  \n",
+       "3   1.0000  \n",
+       "4   1.0000  \n",
+       "5   1.0000  "
       ]
      },
-     "execution_count": 21,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
diff --git a/skpsl/binarizer.py b/skpsl/binarizer.py
index a45a770..7d0f9c4 100644
--- a/skpsl/binarizer.py
+++ b/skpsl/binarizer.py
@@ -1,13 +1,13 @@
 import logging
 from functools import partial
-from typing import Tuple
 
 import numpy as np
-from scipy.optimize import minimize
+import pandas as pd
 from scipy.stats import entropy
 from sklearn.base import TransformerMixin, BaseEstimator
 from sklearn.exceptions import NotFittedError
-from sklearn.preprocessing import MinMaxScaler
+
+from skpsl.data.util import binary_search_optimizer
 
 logger = logging.getLogger()
 
@@ -17,51 +17,49 @@ def __init__(self):
         self.threshs = None
 
     def fit(self, X, y=None):
-        self.threshs = np.apply_along_axis(lambda x: MinEntropyBinarizer._binarize(x, y=y)[0], 0, X)
+        def binarize(x):
+            is_data_binary = set(np.unique(x).astype(int)) == {0, 1}
+            if is_data_binary:
+                return np.nan
+            return binary_search_optimizer(partial(MinEntropyBinarizer._cut_entropy, y=y), x)
+
+        self.threshs = np.apply_along_axis(binarize, 0, X)
         return self
 
     def transform(self, X):
         if self.threshs is None:
             raise NotFittedError()
-        return np.apply_along_axis(lambda x: MinEntropyBinarizer._binarize(x[1:], thresh=x[0])[1], 0,
-                                   np.vstack([self.threshs.reshape(1, -1), X]))
+        threshs = self.threshs
+        threshs[np.isnan(threshs)] = .5
+        return (np.array(X) > threshs).astype(int)
 
-    @staticmethod
-    def _cut_entropy(x: np.ndarray, y: np.ndarray, thresh: float) -> float:
+    def inspect(self, feature_names=None) -> pd.DataFrame:
         """
-        https://www.ijcai.org/Proceedings/93-2/Papers/022.pdf
-        @param x: one-dimensional float array of the feature variable
-        @param y: one-dimensional float array of the target variable
-        @param thresh: scalar or scalar array
-        @return: combined entropy as float
+        Returns a dataframe that visualizes the internal model
         """
-        thresh = np.array(thresh).squeeze()
-        s1 = y[x < thresh]
-        s2 = y[x >= thresh]
-        s1_freqs = np.unique(s1, return_counts=True)[1] / s1.size
-        s2_freqs = np.unique(s2, return_counts=True)[1] / s2.size
-
-        return (s1.size / x.size * entropy(s1_freqs, base=2) +
-                s2.size / x.size * entropy(s2_freqs, base=2))
+        k = len(self.threshs)
+        df = pd.DataFrame(columns=["Threshold"], data=[f">{t:.4f}" for t in self.threshs])
+        if feature_names is not None:
+            df.insert(0, "Feature", feature_names[:k] + [np.nan] * (k - len(feature_names)))
+        return df
 
     @staticmethod
-    def _binarize(x: np.ndarray, y: np.ndarray = None, thresh=None) -> Tuple[float, np.ndarray]:
+    def _cut_entropy(thresh: float, x: np.array, y: np.array) -> float:
         """
-        @param x: one-dimensional float array of the feature variable
-        @param y: one-dimensional float array of the target variable
-        @return: threshold, one-dimensional integer array of binarized feature
+        https://www.ijcai.org/Proceedings/93-2/Papers/022.pdf
+
+        :param x: one-dimensional float array of the feature variable
+        :param y: one-dimensional float array of the target variable
+        :param thresh: scalar
+        :return: combined entropy as float
         """
-        uniq = np.unique(x)
-        if uniq.size < 3:
-            thresh = np.min(uniq) + (np.max(uniq) - np.min(uniq)) / 2
-        if thresh is None:
-            scaler = MinMaxScaler()
-            x_ = scaler.fit_transform(x.reshape(-1, 1)).squeeze()  # type:np.ndarray
-            thresh = minimize(partial(MinEntropyBinarizer._cut_entropy, x_, y), x.mean(), method="Nelder-Mead",
-                              options=dict(initial_simplex=[[0], [1]])).x.squeeze().item()
-            thresh = scaler.inverse_transform([[thresh]]).squeeze()
-            logger.debug(thresh)
-        return thresh, (x >= thresh).astype(int)
+        mask = x > np.array(thresh).squeeze()
+        s1, s2 = y[mask], y[~mask]
+        _, s1_freqs = np.unique(s1, return_counts=True)
+        _, s2_freqs = np.unique(s2, return_counts=True)
+
+        return (s1.size * entropy(s1_freqs, base=2) +
+                s2.size * entropy(s2_freqs, base=2)) / x.size
 
 
 if __name__ == '__main__':
@@ -69,25 +67,13 @@ def _binarize(x: np.ndarray, y: np.ndarray = None, thresh=None) -> Tuple[float,
 
     logging.basicConfig()
     logger.setLevel(logging.DEBUG)
-    # Generating synthetic data with continuous features and a binary target variable
-    X, y = make_classification(n_samples=5, n_features=2, n_informative=2, n_redundant=0, random_state=42)
-
-    print(np.hstack([X, y.reshape(-1, 1)]))
-    print(MinEntropyBinarizer().fit_transform(X, y))
-
-    x1 = np.linspace(0, 100, 10)
-    y = (x1 > 30).astype(int)
-    X = x1.reshape(-1, 1)
-    print(np.hstack([X, y.reshape(-1, 1)]))
-    print(MinEntropyBinarizer().fit_transform(X, y))
-
-    X = np.random.randint(0, 2, size=(10, 3))
-    y = np.random.randint(0, 2, size=(10,))
-    print(np.hstack([X, y.reshape(-1, 1)]))
-    print(X == MinEntropyBinarizer().fit_transform(X, y))
 
     x1 = np.linspace(0, 100, 10)
-    y = (x1 > 30).astype(int)
-    X = x1.reshape(-1, 1)
-    print(np.hstack([X, y.reshape(-1, 1)]))
-    print(MinEntropyBinarizer().fit(X, y).transform(X))
+    y_ = (x1 > 30).astype(int)
+    X_ = x1.reshape(-1, 1)
+    print(np.hstack([X_, y_.reshape(-1, 1)]))
+    print(MinEntropyBinarizer().fit(X_, y_).inspect())
+
+    X_, y_ = make_classification(n_samples=50, n_features=3, n_informative=2, n_redundant=0, random_state=42)
+    print(np.hstack([X_, y_.reshape(-1, 1)]))
+    print(MinEntropyBinarizer().fit(X_, y_).inspect(feature_names=["width", "height"]))
diff --git a/skpsl/data/dataset.py b/skpsl/data/dataset.py
index 2a88a15..efb6281 100644
--- a/skpsl/data/dataset.py
+++ b/skpsl/data/dataset.py
@@ -20,7 +20,7 @@ def lookahead_example():
         score, invscore = [_ClassifierAtK(features=f, scores=s_).fit(X_, y_, ).score(X_) for s_ in permutations(s)]
         if score > invscore:
             score_l1 = psl.score(X_, y_)
-            l2_psl = ProbabilisticScoringList(s).fit(X_, y_, l=2)
+            l2_psl = ProbabilisticScoringList(s).fit(X_, y_, lookahead=2)
             score_l2 = l2_psl.score(X_, y_)
             if score_l2 < score_l1:
                 print(
@@ -45,7 +45,7 @@ def lookahead_example():
             l1_out = l1_psl.score(X_[test_index], y_[test_index])
             l1_in = l1_psl.score(X_[train_index], y_[train_index])
 
-            l2_psl = ProbabilisticScoringList(s).fit(X_[train_index], y_[train_index], l=2)
+            l2_psl = ProbabilisticScoringList(s).fit(X_[train_index], y_[train_index], lookahead=2)
             l2_out = l2_psl.score(X_[test_index], y_[test_index])
             l2_in = l2_psl.score(X_[train_index], y_[train_index])
             if l1_in > l2_in:
diff --git a/skpsl/data/util.py b/skpsl/data/util.py
new file mode 100644
index 0000000..355640d
--- /dev/null
+++ b/skpsl/data/util.py
@@ -0,0 +1,49 @@
+import numpy as np
+from sortedcontainers import SortedSet
+
+
+def binary_search_optimizer(func: callable, data: np.array, minimize=True) -> float:
+    """
+    This algorithm employs a hierarchical logarithmic search to find the global minimum of a parametrized metric.
+
+    Parameter data is used to calculate potential threshold values used for cutting data.
+    It is passed to func without changes. The function is assumed to be pseudo convex in x.
+    This allows for the hierarchical logarithmic search to find a global optimum efficiently.
+
+    :param data: one-dimensional float array of the feature variable
+    :param func: Metric must have the signature (x, data) and return a score to be minimized
+    :return: optimal threshold
+    """
+    values = np.sort(np.unique(data))
+    # Adding the extremal values <min and >max might not be necessary, but it is not trivial to proof.
+    cuts = np.concatenate([[data.min() - 1], (values[:-1] + values[1:]) / 2, [data.max() + 1]])
+    minimal_points = set()
+    min_ = np.inf
+    thresh = None
+    evaluated = SortedSet()
+    to_evaluate = {0, cuts.size - 1}
+
+    while to_evaluate:
+        # evaluate points
+        while to_evaluate:
+            k = to_evaluate.pop()
+            evaluated.add(k)
+            entropy = (-1)**int(not minimize) * func(cuts[k], data)
+            if entropy < min_:
+                min_ = entropy
+                thresh = cuts[k]
+                minimal_points = {k}
+            elif entropy == min_:
+                minimal_points.add(k)
+
+        # calculate new points to evaluate
+        for k in minimal_points:
+            k_pos = evaluated.index(k)
+            candidates = set()
+            for offset in [-1, 1]:
+                try:
+                    candidates.add(k + (evaluated[k_pos + offset] - k) // 2)
+                except IndexError:
+                    pass
+            to_evaluate.update(candidates - evaluated)
+    return thresh
diff --git a/skpsl/probabilistic_scoring_list.py b/skpsl/probabilistic_scoring_list.py
index eb2e2d9..4bc7f77 100644
--- a/skpsl/probabilistic_scoring_list.py
+++ b/skpsl/probabilistic_scoring_list.py
@@ -1,6 +1,6 @@
 import logging
 from itertools import combinations, product, repeat
-from typing import List
+from typing import Optional
 
 import numpy as np
 import pandas as pd
@@ -11,50 +11,81 @@
 from sklearn.isotonic import IsotonicRegression
 from sklearn.metrics import brier_score_loss
 
+from skpsl.data.util import binary_search_optimizer
+
 
 class _ClassifierAtK(BaseEstimator, ClassifierMixin):
     """
     Internal class for the classifier at stage k of the probabilistic scoring list
     """
 
-    def __init__(self, features, scores):
+    def __init__(
+        self,
+        features: tuple[int],
+        scores: tuple[int],
+        initial_thresholds: tuple[Optional[float]],
+    ):
+        """
+
+        :param features: tuple of feature indices. used for selecting data from X
+        :param scores: tuple of scores, corresponding to the feature indices
+        :param initial_thresholds: tuple of thresholds to binarize the feature values
+        """
         self.features = features
         self.scores = scores
+        self.initial_thresholds = initial_thresholds
 
-        self.logger = logging.getLogger(__name__)
+        self.thresholds = list(initial_thresholds)
         self.scores_vec = np.array(scores)
-        self.probabilities = None
+        self.logger = logging.getLogger(__name__)
+        self.calibrator = None
 
     def fit(self, X, y) -> "_ClassifierAtK":
-        scores = self._scores_per_record(X)
+        for i, (f, t) in enumerate(zip(self.features, self.initial_thresholds)):
+            feature_values = X[:, f]
+            is_data_binary = set(np.unique(feature_values).astype(int)) == {0, 1}
+            self.logger.debug(f"feature {f} is non-binary, calculating threshold")
+            if t is np.nan and not is_data_binary:
+                # fit optimal threshold
+                self.thresholds[i] = binary_search_optimizer(
+                    lambda t_, _: self._expected_entropy(
+                        X=X,
+                        y=y,
+                        features=self.features[: i + 1],
+                        scores=self.scores_vec[: i + 1],
+                        thresholds=self.thresholds[:i] + [t_],
+                    ),
+                    feature_values,
+                )
 
-        # compute all possible total scores using subset-summation
-        total_scores = {0}
-        for score in self.scores_vec:
-            total_scores |= {prev_sum + score for prev_sum in total_scores}
-        total_scores = np.array(sorted(total_scores))
+        scores = self._scores_per_binarized_record(
+            X, self.features, self.scores_vec, self.thresholds
+        )
 
         # compute probabilities
-        calibrator = IsotonicRegression(y_min=0.0, y_max=1.0, increasing=True, out_of_bounds="clip")
-        calibrator.fit(scores, y)
-        self.probabilities = {T: p for T, p in zip(total_scores, calibrator.transform(total_scores))}
+        self.calibrator = IsotonicRegression(
+            y_min=0.0, y_max=1.0, increasing=True, out_of_bounds="clip"
+        )
+        self.calibrator.fit(scores, y)
 
         return self
 
     def predict(self, X):
-        if self.probabilities is None:
+        if self.calibrator is None:
             raise NotFittedError()
         return self.predict_proba(X).argmax(axis=1)
 
     def predict_proba(self, X):
-        """Predicts the probability for
         """
-        if self.probabilities is None:
+        Predicts the probability for
+        """
+        if self.calibrator is None:
             raise NotFittedError()
-        scores = self._scores_per_record(X)
-        proba_true = np.empty_like(scores, dtype=float)
-        for total_score in np.unique(scores):
-            proba_true[scores == total_score] = self.probabilities[total_score]
+        proba_true = self.calibrator.transform(
+            self._scores_per_binarized_record(
+                X, self.features, self.scores_vec, self.thresholds
+            )
+        )
         proba = np.vstack([1 - proba_true, proba_true]).T
         return proba
 
@@ -66,17 +97,47 @@ def score(self, X, y=None, sample_weight=None):
         :param sample_weight:
         :return:
         """
-        if self.probabilities is None:
+        if self.calibrator is None:
             raise NotFittedError()
-        scores = self._scores_per_record(X)
+        if not self.features:
+            return stats_entropy(self.calibrator.transform([[0]]))
+        return self._expected_entropy(
+            X,
+            self.features,
+            self.scores_vec,
+            self.thresholds,
+            calibrator=self.calibrator,
+        )
+
+    @staticmethod
+    def _expected_entropy(
+        X, features, scores: np.ndarray, thresholds, calibrator=None, y=None
+    ):
+        scores = _ClassifierAtK._scores_per_binarized_record(
+            X, features, scores, thresholds
+        )
+        if calibrator is None:
+            if y is None:
+                raise AttributeError(
+                    "_expected_entropy must not be called with both 'calibrator' and 'y' being None"
+                )
+            calibrator = IsotonicRegression(
+                y_min=0.0, y_max=1.0, increasing=True, out_of_bounds="clip"
+            )
+            calibrator.fit(scores, y)
         total_scores, score_freqs = np.unique(scores, return_counts=True)
-        score_probas = np.array([self.probabilities[ti] for ti in total_scores])
+        score_probas = np.array(calibrator.transform(total_scores))
         entropy_values = stats_entropy([score_probas, 1 - score_probas], base=2)
         return np.sum((score_freqs / scores.size) * entropy_values)
 
-    # Helper functions
-    def _scores_per_record(self, X):
-        return X[:, self.features] @ self.scores_vec
+    @staticmethod
+    def _scores_per_binarized_record(X, features, scores: np.ndarray, thresholds):
+        if not features:
+            return np.zeros(X.shape[0])
+        data = np.array(X)[:, features]
+        thresholds = np.array(thresholds)
+        thresholds[np.isnan(thresholds)] = 0.5
+        return (data > thresholds[None, :]) @ scores
 
 
 class ProbabilisticScoringList(BaseEstimator, ClassifierMixin):
@@ -85,8 +146,11 @@ class ProbabilisticScoringList(BaseEstimator, ClassifierMixin):
     A probabilistic classifier that greedily creates a PSL selecting one feature at a time
     """
 
-    def __init__(self, score_set, entropy_threshold=-1):
-        """ IMPORTANT: Shannon entropy is calculated with respect to base 2
+    def __init__(self, score_set: set, entropy_threshold: float = -1):
+        """
+
+        :param score_set: Set score values to be considered. Basically feature weights.
+        :param entropy_threshold: Shannon Entropy base 2 threshold after which to stop fitting more stages.
         """
         self.score_set = score_set
         self.entropy_threshold = entropy_threshold
@@ -94,15 +158,17 @@ def __init__(self, score_set, entropy_threshold=-1):
         self.logger = logging.getLogger(__name__)
         self.sorted_score_set = sorted(self.score_set, reverse=True)
         self._stage_clf = _ClassifierAtK
-        self.stage_clfs = []  # type: List[_ClassifierAtK]
+        self.stage_clfs = []  # type: list[_ClassifierAtK]
 
-    def fit(self, X, y, l=1, n_jobs=1, predef_features=None, predef_scores=None) -> "ProbabilisticScoringList":
+    def fit(
+        self, X, y, lookahead=1, n_jobs=1, predef_features=None, predef_scores=None
+    ) -> "ProbabilisticScoringList":
         """
         Fits a probabilistic scoring list to the given data
 
         :param X:
         :param y:
-        :param l: steps of look ahead
+        :param lookahead: steps of look ahead
         :param n_jobs: passed to joblib for parallelization
         :param predef_features:
         :param predef_scores:
@@ -113,32 +179,65 @@ def fit(self, X, y, l=1, n_jobs=1, predef_features=None, predef_scores=None) ->
         remaining_features = list(range(number_features))
 
         if predef_features is not None and len(predef_features) != number_features:
-            raise ValueError("Predefined features must be a permutation of all features!")
+            raise ValueError(
+                "Predefined features must be a permutation of all features!"
+            )
 
         stage = 0
 
-        # first 
+        # first
         expected_entropy = self._fit_and_store_clf_at_k(X, y)
 
         while remaining_features and expected_entropy > self.entropy_threshold:
             stage += 1
 
-            features_to_consider = remaining_features if predef_features is None else [predef_features[stage - 1]]
-            scores_to_consider = self.sorted_score_set if predef_scores is None else [predef_scores[stage - 1]]
-
-            entropies, f, s = zip(*Parallel(n_jobs=n_jobs)(
-                delayed(self._optimize)(self.features, f_seq, self.scores, list(s_seq), self._stage_clf, X, y)
-                for (f_seq, s_seq) in product(
-                    self._gen_lookahead(features_to_consider, l),
-                    # cartesian power of scores
-                    product(*repeat(scores_to_consider, min(l, len(features_to_consider))))
+            # noinspection PyUnresolvedReferences
+            features_to_consider = (
+                remaining_features
+                if predef_features is None
+                else [predef_features[stage - 1]]
+            )
+            scores_to_consider = (
+                self.sorted_score_set
+                if predef_scores is None
+                else [predef_scores[stage - 1]]
+            )
+
+            entropies, f, s, t = zip(
+                *Parallel(n_jobs=n_jobs)(
+                    delayed(self._optimize)(
+                        self.features,
+                        f_seq,
+                        self.scores,
+                        list(s_seq),
+                        self.thresholds,
+                        self._stage_clf,
+                        X,
+                        y,
+                    )
+                    for (f_seq, s_seq) in product(
+                        self._gen_lookahead(features_to_consider, lookahead),
+                        # cartesian power of scores
+                        product(
+                            *repeat(
+                                scores_to_consider,
+                                min(lookahead, len(features_to_consider)),
+                            )
+                        ),
+                    )
                 )
-            ))
+            )
 
             i = np.argmin(entropies)
             remaining_features.remove(f[i])
 
-            expected_entropy = self._fit_and_store_clf_at_k(X, y, self.features + [f[i]], self.scores + [s[i]])
+            expected_entropy = self._fit_and_store_clf_at_k(
+                X,
+                y,
+                self.features + [f[i]],
+                self.scores + [s[i]],
+                self.thresholds + [t[i]],
+            )
         return self
 
     def predict(self, X, k=-1):
@@ -159,7 +258,9 @@ def predict_proba(self, X, k=-1):
         :return:
         """
         if not self.stage_clfs:
-            raise NotFittedError("Please fit the probabilistic scoring classifier before usage.")
+            raise NotFittedError(
+                "Please fit the probabilistic scoring classifier before usage."
+            )
 
         return self.stage_clfs[k].predict_proba(X)
 
@@ -184,14 +285,48 @@ def inspect(self, k=None, feature_names=None) -> pd.DataFrame:
         """
         k = k or len(self.stage_clfs) - 1
 
-        pmfs = [clf.probabilities for clf in self.stage_clfs[:k + 1]]
-        all_total_scores = sorted(set.union(*[set(pmf.keys()) for pmf in pmfs]))
-        data = [[pmfs[i].get(t_, np.nan) for t_ in all_total_scores] for i in range(k + 1)]
-
-        df = pd.DataFrame(columns=[f"T = {t_}" for t_ in all_total_scores], data=data)
-        df.insert(0, "Score", [np.nan] + self.stage_clfs[k].scores)
+        scores = self.stage_clfs[k].scores
+        features = self.stage_clfs[k].features
+        thresholds = self.stage_clfs[k].thresholds
+
+        all_total_scores = [{0}]
+        for i, score in enumerate(scores):
+            all_total_scores.append(
+                all_total_scores[i]
+                | {prev_sum + score for prev_sum in all_total_scores[i]}
+            )
+
+        data = []
+        for clf, T in zip(self.stage_clfs[: k + 1], all_total_scores):
+            a = {t: np.nan for t in all_total_scores[-1]}
+            probas = clf.calibrator.transform(np.array(list(T)))
+            a.update(dict(zip(T, probas)))
+            data.append(a)
+
+        df = pd.DataFrame(data)
+        df = df[sorted(df.columns)]
+        df.columns = [f"T = {t_}" for t_ in df.columns]
+        df.insert(0, "Score", np.array([np.nan] + list(scores)))
         if feature_names is not None:
-            df.insert(0, "Feature", [np.nan] + feature_names[:k] + [np.nan] * (k - len(feature_names)))
+            if len(feature_names) != len(features):
+                raise ValueError(
+                    f"Passed feature names are of incorrect length! Passed {len(feature_names)}, expected {len(features)}."
+                )
+            feature_names = [feature_names[i] for i in features]
+            df.insert(0, "Feature", [np.nan] + feature_names[:k])
+        else:
+            df.insert(
+                0,
+                "Feature Index",
+                [np.nan] + features[:k],
+            )
+        if not all(np.isnan(t) for t in thresholds):
+            df.insert(
+                0,
+                "Threshold",
+                [np.nan]
+                + [(np.nan if np.isnan(t) else f">{t:.4f}") for t in thresholds],
+            )
         return df.reset_index(names=["Stage"])
 
     @property
@@ -202,33 +337,50 @@ def features(self):
     def scores(self):
         return self.stage_clfs[-1].scores if self.stage_clfs else []
 
-    def _fit_and_store_clf_at_k(self, X, y, f=None, s=None):
-        f, s = f or [], s or []
-        k_clf = self._stage_clf(features=f, scores=s).fit(X, y)
+    @property
+    def thresholds(self):
+        return self.stage_clfs[-1].thresholds if self.stage_clfs else []
+
+    def _fit_and_store_clf_at_k(self, X, y, f=None, s=None, t=None):
+        f, s, t = f or [], s or [], t or []
+        k_clf = self._stage_clf(features=f, scores=s, initial_thresholds=t).fit(X, y)
         self.stage_clfs.append(k_clf)
         return k_clf.score(X)
 
     @staticmethod
-    def _optimize(features, feature_extension, scores, score_extension, clfcls, X, y):
-        clf = clfcls(features=features + feature_extension, scores=scores + score_extension).fit(X, y)
-        return clf.score(X), feature_extension[0], score_extension[0]
+    def _optimize(
+        features, feature_extension, scores, score_extension, thresholds, clfcls, X, y
+    ):
+        clf = clfcls(
+            features=features + feature_extension,
+            scores=scores + score_extension,
+            initial_thresholds=thresholds + [np.nan] * len(feature_extension),
+        ).fit(X, y)
+        return (
+            clf.score(X),
+            feature_extension[0],
+            score_extension[0],
+            clf.thresholds[len(features)],
+        )
 
     @staticmethod
     def _gen_lookahead(list_, lookahead):
         # generate sequences of shortening lookaheads (because combinations returns empty list if len(list) < l)
-        combination_seqs = ([list(tup) for tup in combinations(list_, _l)] for _l in range(lookahead, 0, -1))
+        combination_seqs = (
+            [list(tup) for tup in combinations(list_, _l)]
+            for _l in range(lookahead, 0, -1)
+        )
         # get first non-empty sequence
         seqs = next((seq for seq in combination_seqs if seq))
         return seqs
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     from sklearn.datasets import make_classification
     from sklearn.model_selection import cross_val_score
 
     # Generating synthetic data with continuous features and a binary target variable
-    X, y = make_classification(random_state=42)
-    X = (X > .5).astype(int)
+    X_, y_ = make_classification(random_state=42)
 
-    clf = ProbabilisticScoringList([-1, 1, 2])
-    print("Brier score:", cross_val_score(clf, X, y, cv=5).mean())
+    clf_ = ProbabilisticScoringList({-1, 1, 2})
+    print("Brier score:", cross_val_score(clf_, X_, y_, cv=5).mean())


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/test_probabilistic_scoring_list.py b/tests/test_probabilistic_scoring_list.py
new file mode 100644
index 0000000..52fb7fe
--- /dev/null
+++ b/tests/test_probabilistic_scoring_list.py
@@ -0,0 +1,45 @@
+import numpy as np
+from sklearn.datasets import make_classification
+from sklearn.model_selection import train_test_split, cross_val_score
+from sklearn.pipeline import Pipeline
+
+from skpsl import ProbabilisticScoringList, MinEntropyBinarizer
+
+
+def test_binary_data():
+    X, y = make_classification(random_state=42)
+    X = (X > .5).astype(int)
+
+    psl = ProbabilisticScoringList({-1, 1, 2})
+
+    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
+    psl.fit(X_train, y_train)
+    assert psl.thresholds == [np.nan] * X.shape[1]
+
+
+def test_inspect():
+    # Generating synthetic data with continuous features and a binary target variable
+    X, y = make_classification(n_features=5, random_state=42)
+    X = (X > .5).astype(int)
+
+    psl = ProbabilisticScoringList({-1, 1, 2})
+
+    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2, random_state=42)
+    psl.fit(X_train, y_train)
+    df = psl.inspect(3)
+    print(df.to_string(index=False, na_rep="-", justify="center", float_format=lambda x: f"{x:.2f}"))
+
+
+def test_improvement_for_internalized_binarization():
+    X, y = make_classification(n_features=6, n_informative=4, random_state=42)
+
+    clf_ = Pipeline([("binarizer", MinEntropyBinarizer()),
+                     ("psl", ProbabilisticScoringList({-1, 1, 2}))])
+    pipe_score = cross_val_score(clf_, X, y, cv=5).mean()
+
+    clf_ = ProbabilisticScoringList({-1, 1, 2})
+    psl_score = cross_val_score(clf_, X, y, cv=5).mean()
+
+    # lower score is better
+    # the internalized binarization should perform better
+    assert psl_score < pipe_score