-
Notifications
You must be signed in to change notification settings - Fork 473
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.
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
Uncertainty: Conformal Prediction V1.1 - extend to multiple forecast steps instead of only a single forecast step #1073
Changes from all commits
fc1936a
322b239
e084a41
d2be041
e366bb6
d6b229a
bd4b225
feef1d2
fa5e861
54458bf
ed989ff
d2863a7
889a8fa
4fbe1b5
195fe60
9c91965
d491bea
8ad5088
1a9e165
2aebf76
2a1ddf4
f4629e2
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -5,7 +5,10 @@ | |
import numpy as np | ||
import pandas as pd | ||
|
||
from neuralprophet.plot_forecast_matplotlib import plot_nonconformity_scores | ||
from neuralprophet.plot_forecast_matplotlib import plot_interval_width_per_timestep, plot_nonconformity_scores | ||
from neuralprophet.plot_forecast_plotly import ( | ||
plot_interval_width_per_timestep as plot_interval_width_per_timestep_plotly, | ||
) | ||
from neuralprophet.plot_forecast_plotly import plot_nonconformity_scores as plot_nonconformity_scores_plotly | ||
|
||
|
||
|
@@ -23,13 +26,16 @@ class Conformal: | |
Options | ||
* ``naive``: Naive or Absolute Residual | ||
* ``cqr``: Conformalized Quantile Regression | ||
n_forecasts : int | ||
optional, number of steps ahead of prediction time step to forecast | ||
quantiles : list | ||
optional, list of quantiles for quantile regression uncertainty estimate | ||
|
||
""" | ||
|
||
alpha: float | ||
method: str | ||
n_forecasts: int = 1 | ||
quantiles: Optional[List[float]] = None | ||
|
||
def predict(self, df: pd.DataFrame, df_cal: pd.DataFrame) -> pd.DataFrame: | ||
|
@@ -48,34 +54,50 @@ def predict(self, df: pd.DataFrame, df_cal: pd.DataFrame) -> pd.DataFrame: | |
test dataframe with uncertainty prediction intervals | ||
|
||
""" | ||
# conformalize | ||
self.noncon_scores = self._get_nonconformity_scores(df_cal) | ||
self.q_hat = self._get_q_hat(df_cal) | ||
df["qhat1"] = self.q_hat | ||
if self.method == "naive": | ||
df["yhat1 - qhat1"] = df["yhat1"] - self.q_hat | ||
df["yhat1 + qhat1"] = df["yhat1"] + self.q_hat | ||
elif self.method == "cqr": | ||
quantile_hi = str(max(self.quantiles) * 100) | ||
quantile_lo = str(min(self.quantiles) * 100) | ||
df[f"yhat1 {quantile_hi}% - qhat1"] = df[f"yhat1 {quantile_hi}%"] - self.q_hat | ||
df[f"yhat1 {quantile_hi}% + qhat1"] = df[f"yhat1 {quantile_hi}%"] + self.q_hat | ||
df[f"yhat1 {quantile_lo}% - qhat1"] = df[f"yhat1 {quantile_lo}%"] - self.q_hat | ||
df[f"yhat1 {quantile_lo}% + qhat1"] = df[f"yhat1 {quantile_lo}%"] + self.q_hat | ||
else: | ||
raise ValueError( | ||
f"Unknown conformal prediction method '{self.method}'. Please input either 'naive' or 'cqr'." | ||
) | ||
self.q_hats = [] | ||
for step_number in range(1, self.n_forecasts + 1): | ||
# conformalize | ||
noncon_scores = self._get_nonconformity_scores(df_cal, step_number) | ||
q_hat = self._get_q_hat(df_cal, noncon_scores) | ||
df[f"qhat{step_number}"] = q_hat | ||
if self.method == "naive": | ||
df[f"yhat{step_number} - qhat{step_number}"] = df[f"yhat{step_number}"] - q_hat | ||
df[f"yhat{step_number} + qhat{step_number}"] = df[f"yhat{step_number}"] + q_hat | ||
elif self.method == "cqr": | ||
quantile_hi = str(max(self.quantiles) * 100) | ||
quantile_lo = str(min(self.quantiles) * 100) | ||
df[f"yhat{step_number} {quantile_hi}% - qhat{step_number}"] = ( | ||
df[f"yhat{step_number} {quantile_hi}%"] - q_hat | ||
) | ||
df[f"yhat{step_number} {quantile_hi}% + qhat{step_number}"] = ( | ||
df[f"yhat{step_number} {quantile_hi}%"] + q_hat | ||
) | ||
df[f"yhat{step_number} {quantile_lo}% - qhat{step_number}"] = ( | ||
df[f"yhat{step_number} {quantile_lo}%"] - q_hat | ||
) | ||
df[f"yhat{step_number} {quantile_lo}% + qhat{step_number}"] = ( | ||
df[f"yhat{step_number} {quantile_lo}%"] + q_hat | ||
) | ||
else: | ||
raise ValueError( | ||
f"Unknown conformal prediction method '{self.method}'. Please input either 'naive' or 'cqr'." | ||
) | ||
if step_number == 1: | ||
# save nonconformity scores of the first timestep | ||
self.noncon_scores = noncon_scores | ||
self.q_hats.append(q_hat) | ||
|
||
return df | ||
|
||
def _get_nonconformity_scores(self, df_cal: pd.DataFrame) -> np.ndarray: | ||
def _get_nonconformity_scores(self, df_cal: pd.DataFrame, step_number: int) -> np.ndarray: | ||
"""Get the nonconformity scores using the given conformal prediction technique. | ||
|
||
Parameters | ||
---------- | ||
df_cal : pd.DataFrame | ||
calibration dataframe | ||
step_number : int | ||
i-th step ahead forecast | ||
|
||
Returns | ||
------- | ||
|
@@ -89,35 +111,40 @@ def _get_nonconformity_scores(self, df_cal: pd.DataFrame) -> np.ndarray: | |
quantile_lo = str(min(self.quantiles) * 100) | ||
cqr_scoring_func = ( | ||
lambda row: [None, None] | ||
if row[f"yhat1 {quantile_lo}%"] is None or row[f"yhat1 {quantile_hi}%"] is None | ||
if row[f"yhat{step_number} {quantile_lo}%"] is None or row[f"yhat{step_number} {quantile_hi}%"] is None | ||
else [ | ||
max( | ||
row[f"yhat1 {quantile_lo}%"] - row["y"], | ||
row["y"] - row[f"yhat1 {quantile_hi}%"], | ||
row[f"yhat{step_number} {quantile_lo}%"] - row["y"], | ||
row["y"] - row[f"yhat{step_number} {quantile_hi}%"], | ||
), | ||
0 if row[f"yhat1 {quantile_lo}%"] - row["y"] > row["y"] - row[f"yhat1 {quantile_hi}%"] else 1, | ||
0 | ||
if row[f"yhat{step_number} {quantile_lo}%"] - row["y"] | ||
> row["y"] - row[f"yhat{step_number} {quantile_hi}%"] | ||
else 1, | ||
] | ||
) | ||
scores_df = df_cal.apply(cqr_scoring_func, axis=1, result_type="expand") | ||
scores_df.columns = ["scores", "arg"] | ||
noncon_scores = scores_df["scores"].values | ||
else: # self.method == "naive" | ||
# Naive nonconformity scoring function | ||
noncon_scores = abs(df_cal["y"] - df_cal["yhat1"]).values | ||
noncon_scores = abs(df_cal["y"] - df_cal[f"yhat{step_number}"]).values | ||
# Remove NaN values | ||
noncon_scores = noncon_scores[~pd.isnull(noncon_scores)] | ||
# Sort | ||
noncon_scores.sort() | ||
|
||
return noncon_scores | ||
|
||
def _get_q_hat(self, df_cal: pd.DataFrame) -> float: | ||
def _get_q_hat(self, df_cal: pd.DataFrame, noncon_scores: np.ndarray) -> float: | ||
"""Get the q_hat that is derived from the nonconformity scores. | ||
|
||
Parameters | ||
---------- | ||
df_cal : pd.DataFrame | ||
calibration dataframe | ||
noncon_scores : np.ndarray | ||
nonconformity scores | ||
|
||
Returns | ||
------- | ||
|
@@ -126,8 +153,8 @@ def _get_q_hat(self, df_cal: pd.DataFrame) -> float: | |
|
||
""" | ||
# Get the q-hat index and value | ||
q_hat_idx = int(len(self.noncon_scores) * self.alpha) | ||
q_hat = self.noncon_scores[-q_hat_idx] | ||
q_hat_idx = int(len(noncon_scores) * self.alpha) | ||
q_hat = noncon_scores[-q_hat_idx] | ||
|
||
return q_hat | ||
|
||
|
@@ -146,8 +173,16 @@ def plot(self, plotting_backend: str): | |
""" | ||
method = self.method.upper() if "cqr" in self.method.lower() else self.method.title() | ||
if plotting_backend == "plotly": | ||
fig = plot_nonconformity_scores_plotly(self.noncon_scores, self.alpha, self.q_hat, method) | ||
if self.n_forecasts == 1: | ||
# includes nonconformity scores of the first timestep | ||
fig = plot_nonconformity_scores_plotly(self.noncon_scores, self.alpha, self.q_hats[0], method) | ||
else: | ||
fig = plot_interval_width_per_timestep_plotly(self.q_hats, method) | ||
elif plotting_backend == "matplotlib": | ||
fig = plot_nonconformity_scores(self.noncon_scores, self.alpha, self.q_hat, method) | ||
if self.n_forecasts == 1: | ||
# includes nonconformity scores of the first timestep | ||
fig = plot_nonconformity_scores(self.noncon_scores, self.alpha, self.q_hats[0], method) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Ideally we would have the same behavior for There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Further, there seems to be these two options only in the matplotlib case but not for plotly - However the plot should be as identical as possible independent of the plotting backend. How can we make them return the same kind of plot? There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Two different plots show whether n_forecasts == 1 or n_forecasts >1. See your 1st and 3rd questions for details. As for the matplotlib and plotly, both have the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Just added There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. It makes more sense to me now, thank you! |
||
else: | ||
fig = plot_interval_width_per_timestep(self.q_hats, method) | ||
if plotting_backend in ["matplotlib", "plotly"] and matplotlib.is_interactive(): | ||
fig.show() |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -443,7 +443,7 @@ def plot_multiforecast_component( | |
|
||
|
||
def plot_nonconformity_scores(scores, alpha, q, method): | ||
"""Plot the NeuralProphet forecast components. | ||
"""Plot the nonconformity scores as well as the one-sided interval width (q). | ||
|
||
Parameters | ||
---------- | ||
|
@@ -470,8 +470,37 @@ def plot_nonconformity_scores(scores, alpha, q, method): | |
ax.plot(confidence_levels, scores, label="score") | ||
ax.axvline(x=1 - alpha, color="g", linestyle="-", label=f"(1-alpha) = {1-alpha}", linewidth=1) | ||
ax.axhline(y=q, color="r", linestyle="-", label=f"q1 = {round(q, 2)}", linewidth=1) | ||
ax.set_title(f"{method} One-Sided Interval Width with q") | ||
ax.set_xlabel("Confidence Level") | ||
ax.set_ylabel("One-Sided Interval Width") | ||
ax.set_title(f"{method} One-Sided Interval Width with q") | ||
ax.legend() | ||
return fig | ||
|
||
|
||
def plot_interval_width_per_timestep(q_hats, method): | ||
"""Plot the nonconformity scores as well as the one-sided interval width (q). | ||
|
||
Parameters | ||
---------- | ||
q_hats : list | ||
prediction interval widths (or q) for each timestep | ||
method : str | ||
name of conformal prediction technique used | ||
|
||
Options | ||
* (default) ``naive``: Naive or Absolute Residual | ||
* ``cqr``: Conformalized Quantile Regression | ||
|
||
Returns | ||
------- | ||
matplotlib.pyplot.figure | ||
Figure showing the q-values for each timestep | ||
""" | ||
fig, ax = plt.subplots() | ||
ax.plot(range(1, len(q_hats) + 1), q_hats) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I am sorry, I am a bit confused about what is being plotted based on the docstring and code. There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This is the You can see how q linearly increases the further the timestep (until it plateaus around t+12)? This shows the greater uncertainty, the further out the forecast that it requires wider intervals for the same confidence level (e.g., 90% in this example). The reason for this plateau could be attributed to the half-day seasonality of this hospital energy load dataset (i.e., day and night usage). There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Thank you for explaining! |
||
ax.set_title(f"{method} One-Sided Interval Width with q per Timestep") | ||
ax.set_xlabel("Timestep Number") | ||
ax.set_ylabel("One-Sided Interval Width") | ||
# ax.set_xlim(left=1) | ||
ax.set_ylim(bottom=0) | ||
return fig |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why are these saved for the first step (only)?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The nonconformity scores are saved in order to be inputted into
plot_nonconformity_scores()
, where it plots:The nonconformity score is the score (blue line).
This plot is only called when n_forecasts ==1. For models with n_forecasts >1, it will call the
plot_interval_width_per_timestep()
instead (more details under your 3rd question), which doesn't require nonconformity scores. Therefore, only the first step needs to be saved.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for the explanation!