update lorenz curve

Author: glemaitre

content/python_files/feature_engineering.py

@@ -629,6 +629,12 @@ def splitter(X, y, index_generator):
 
 def lorenz_curve(observed_value, predicted_value, n_samples=1_000):
     """Compute the Lorenz curve for a given true and predicted values."""
+
+    def gini_index(cum_proportion_population, cum_proportion_y_true):
+        from sklearn.metrics import auc
+
+        return 1 - 2 * auc(cum_proportion_population, cum_proportion_y_true)
+
     observed_value = np.asarray(observed_value)
     predicted_value = np.asarray(predicted_value)
 
@@ -642,98 +648,113 @@ def lorenz_curve(observed_value, predicted_value, n_samples=1_000):
     cum_proportion_y_true = np.cumsum(observed_value_sorted)
     cum_proportion_y_true /= cum_proportion_y_true[-1]
 
+    gini_model = gini_index(cum_proportion_population, cum_proportion_y_true)
+
     cum_proportion_population_interpolated = np.linspace(0, 1, n_samples)
     cum_proportion_y_true_interpolated = np.interp(
         cum_proportion_population_interpolated,
         cum_proportion_population,
         cum_proportion_y_true,
     )
 
-    return cum_proportion_population_interpolated, cum_proportion_y_true_interpolated
+    return pl.DataFrame(
+        {
+            "cum_population": cum_proportion_population_interpolated,
+            "cum_observed": cum_proportion_y_true_interpolated,
+        }
+    ).with_columns(
+        pl.lit(gini_model).alias("gini_index"),
+    )
 
 
-def plot_lorenz_curve(observed_value, predicted_value, n_samples=1_000):
+def plot_lorenz_curve(cv_predictions, n_samples=1_000):
     """Plot the Lorenz curve for a given true and predicted values."""
 
-    def gini_index(cum_proportion_population, cum_proportion_y_true):
-        from sklearn.metrics import auc
+    results = []
+    for cv_idx, predictions in enumerate(cv_predictions):
+        results.append(
+            lorenz_curve(
+                observed_value=predictions["load_mw"],
+                predicted_value=predictions["predicted_load_mw"],
+                n_samples=n_samples,
+            ).with_columns(
+                pl.lit(cv_idx).alias("cv_idx"),
+                pl.lit("model").alias("model"),
+            )
+        )
 
-        return 1 - 2 * auc(cum_proportion_population, cum_proportion_y_true)
+        results.append(
+            lorenz_curve(
+                observed_value=predictions["load_mw"],
+                predicted_value=predictions["load_mw"],
+                n_samples=n_samples,
+            ).with_columns(
+                pl.lit(cv_idx).alias("cv_idx"),
+                pl.lit("oracle").alias("model"),
+            )
+        )
 
-    cum_population_model, cum_observed_model = lorenz_curve(
-        observed_value, predicted_value, n_samples
-    )
-    gini_model = gini_index(cum_population_model, cum_observed_model)
+    results = pl.concat(results)
 
-    cum_population_oracle, cum_observed_oracle = lorenz_curve(
-        observed_value, observed_value, n_samples
+    gini_stats = results.group_by("model").agg(
+        [
+            pl.col("gini_index")
+            .mean()
+            .map_elements(lambda x: f"{x:.4f}", return_dtype=pl.String)
+            .alias("gini_mean"),
+            pl.col("gini_index")
+            .std()
+            .map_elements(lambda x: f"{x:.4f}", return_dtype=pl.String)
+            .alias("gini_std_dev"),
+        ]
     )
-    gini_oracle = gini_index(cum_population_oracle, cum_observed_oracle)
 
-    model_chart = (
-        altair.Chart(
-            pl.DataFrame(
-                {
-                    "cum_population": cum_population_model,
-                    "cum_observed": cum_observed_model,
-                    "model": f"Model (Gini index: {gini_model:.4f})",
-                }
-            )
-        )
-        .mark_line(strokeDash=[6, 3, 2, 3], tooltip=True)
-        .encode(
-            x=altair.X(
-                "cum_population:Q",
-                title="Fraction of observations sorted by predicted label",
-            ),
-            y=altair.Y("cum_observed:Q", title="Cumulative observed load proportion"),
-            color=altair.Color("model:N", legend=altair.Legend(title="Models")),
+    results = results.join(gini_stats, on="model").with_columns(
+        pl.format("{} ({} +/- {})", "model", "gini_mean", "gini_std_dev").alias(
+            "model_label"
         )
     )
 
-    oracle_chart = (
+    diagonal_chart = (
         altair.Chart(
             pl.DataFrame(
                 {
-                    "cum_population": cum_population_oracle,
-                    "cum_observed": cum_observed_oracle,
-                    "model": f"Oracle (Gini index: {gini_oracle:.4f})",
+                    "cum_population": [0, 1],
+                    "cum_observed": [0, 1],
+                    "model_label": "Non-informative model",
                 }
             )
         )
-        .mark_line(strokeDash=[4, 4], tooltip=True)
+        .mark_line(strokeDash=[4, 4], opacity=0.5, tooltip=True)
         .encode(
             x=altair.X(
                 "cum_population:Q",
                 title="Fraction of observations sorted by predicted label",
             ),
             y=altair.Y("cum_observed:Q", title="Cumulative observed load proportion"),
-            color=altair.Color("model:N", legend=altair.Legend(title="Models")),
+            color=altair.Color("model_label:N", legend=altair.Legend(title="Models")),
         )
     )
 
-    diagonal_chart = (
-        altair.Chart(
-            pl.DataFrame(
-                {
-                    "cum_population": [0, 1],
-                    "cum_observed": [0, 1],
-                    "model": "Non-informative model"
-                }
-            )
-        )
-        .mark_line(strokeDash=[4, 4], opacity=0.5, tooltip=True)
+    model_chart = (
+        altair.Chart(results)
+        .mark_line(opacity=0.3, tooltip=True)
         .encode(
             x=altair.X(
                 "cum_population:Q",
-                title="Fraction of observations sorted by predicted label"
+                title="Fraction of observations sorted by predicted label",
             ),
             y=altair.Y("cum_observed:Q", title="Cumulative observed load proportion"),
-            color=altair.Color("model:N", legend=altair.Legend(title="Models"))
+            color=altair.Color("model_label:N", legend=altair.Legend(title="Models")),
+            detail="cv_idx:N",
         )
     )
 
-    return model_chart + oracle_chart + diagonal_chart
+    return model_chart + diagonal_chart
+
+
+plot_lorenz_curve(cv_predictions, n_samples=500).interactive()
+
 
 # %%
 def plot_reliability_diagram(cv_predictions, n_bins=10):
@@ -1181,4 +1202,4 @@ def scoring(regressor, X, y):
 
     display(chart)
 
-# %%
+# %%