From b88b2c4144ca4b87430c224e19d28ce5c20ec028 Mon Sep 17 00:00:00 2001 From: Marco Zanotti Date: Fri, 22 Nov 2024 14:41:36 +0100 Subject: [PATCH 1/8] add MSSE and RMSSE losses --- utilsforecast/losses.py | 86 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 85 insertions(+), 1 deletion(-) diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index c46d199..444cac2 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -1,7 +1,7 @@ # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/losses.ipynb. # %% auto 0 -__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'quantile_loss', 'mqloss', 'coverage', 'calibration', +__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', 'mqloss', 'coverage', 'calibration', 'scaled_crps'] # %% ../nbs/losses.ipynb 3 @@ -374,6 +374,90 @@ def gen_expr(model, baseline) -> pl_Expr: res = res.select([id_col, *exprs]) return res + +def msse( + df: DFType, + models: List[str], + seasonality: int, + train_df: DFType, + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """Mean Squared Scaled Error (MSSE) + + MSSE measures the relative prediction + accuracy of a forecasting method by comparinng the mean squared errors + of the prediction and the observed value against the mean + squared errors of the seasonal naive model. + + Parameters + ---------- + df : pandas or polars DataFrame + Input dataframe with id, actuals and predictions. + models : list of str + Columns that identify the models predictions. + seasonality : int + Main frequency of the time series; + Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. + train_df : pandas or polars DataFrame + Training dataframe with id and actual values. Must be sorted by time. + id_col : str (default='unique_id') + Column that identifies each serie. + target_col : str (default='y') + Column that contains the target. + + Returns + ------- + pandas or polars Dataframe + dataframe with one row per id and one column per model. + + References + ---------- + [1] https://otexts.com/fpp3/accuracy.html + """ + mean_sq_err = mse(df, models, id_col, target_col) + if isinstance(train_df, pd.DataFrame): + mean_sq_err = mean_sq_err.set_index(id_col) + # assume train_df is sorted + lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) + scale = train_df[target_col].sub(lagged).pow(2) + scale = scale.groupby(train_df[id_col], observed=True).mean() + res = mean_sq_err.div(_zero_to_nan(scale), axis=0).fillna(0) + res.index.name = id_col + res = res.reset_index() + else: + # assume train_df is sorted + lagged = pl.col(target_col).shift(seasonality).over(id_col) + scale_expr = pl.col(target_col).sub(lagged).pow(2).alias("scale") + scale = train_df.select([id_col, scale_expr]) + scale = ufp.group_by(scale, id_col).mean() + scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) + + def gen_expr(model): + return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) + + full_df = mean_sq_err.join(scale, on=id_col, how="left") + res = _pl_agg_expr(full_df, models, id_col, gen_expr) + return res + + +def rmsse( + df: DFType, + models: List[str], + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """Root Mean Squared Scaled Error (RMSE) + + RMSSE is the square root of the MSSE""" + res = msse(df, models, id_col, target_col) + if isinstance(res, pd.DataFrame): + res[models] = res[models].pow(0.5) + else: + res = res.with_columns(*[pl.col(c).pow(0.5) for c in models]) + return res + + # %% ../nbs/losses.ipynb 55 def quantile_loss( df: DFType, From 7c611c1323a6462edb9ffb78dcb7b4ec629b2b65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Morales?= Date: Fri, 22 Nov 2024 10:56:40 -0600 Subject: [PATCH 2/8] update notebook --- nbs/losses.ipynb | 341 ++++++++++++++++++++++++++++++++++++--- utilsforecast/_modidx.py | 2 + utilsforecast/losses.py | 40 +++-- 3 files changed, 350 insertions(+), 33 deletions(-) diff --git a/nbs/losses.ipynb b/nbs/losses.ipynb index ed71bc7..04ee830 100644 --- a/nbs/losses.ipynb +++ b/nbs/losses.ipynb @@ -613,6 +613,8 @@ "text/markdown": [ "---\n", "\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L146){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "\n", "#### bias\n", "\n", "> bias (df:~DFType, models:List[str], id_col:str='unique_id',\n", @@ -633,6 +635,8 @@ "text/plain": [ "---\n", "\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L146){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "\n", "#### bias\n", "\n", "> bias (df:~DFType, models:List[str], id_col:str='unique_id',\n", @@ -771,7 +775,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L154){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L181){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mape\n", "\n", @@ -798,7 +802,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L154){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L181){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mape\n", "\n", @@ -909,7 +913,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L192){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L219){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### smape\n", "\n", @@ -938,7 +942,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L192){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L219){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### smape\n", "\n", @@ -1099,7 +1103,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L229){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L256){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mase\n", "\n", @@ -1128,7 +1132,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L229){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L256){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mase\n", "\n", @@ -1264,7 +1268,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L297){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L324){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### rmae\n", "\n", @@ -1289,7 +1293,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L297){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L324){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### rmae\n", "\n", @@ -1335,6 +1339,307 @@ ")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Mean Squared Scaled Error\n", + "\n", + "$$\n", + "\\mathrm{MSSE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}, \\mathbf{\\hat{y}}^{season}_{\\tau}) = \n", + "\\frac{1}{H} \\sum^{t+H}_{\\tau=t+1} \\frac{(y_{\\tau}-\\hat{y}_{\\tau})^2}{\\mathrm{MSE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}^{season}_{\\tau})}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| export\n", + "def msse(\n", + " df: DFType,\n", + " models: List[str],\n", + " seasonality: int,\n", + " train_df: DFType,\n", + " id_col: str = \"unique_id\",\n", + " target_col: str = \"y\",\n", + ") -> DFType:\n", + " \"\"\"Mean Squared Scaled Error (MSSE)\n", + "\n", + " MSSE measures the relative prediction\n", + " accuracy of a forecasting method by comparinng the mean squared errors\n", + " of the prediction and the observed value against the mean\n", + " squared errors of the seasonal naive model.\n", + "\n", + " Parameters\n", + " ----------\n", + " df : pandas or polars DataFrame\n", + " Input dataframe with id, actuals and predictions.\n", + " models : list of str\n", + " Columns that identify the models predictions.\n", + " seasonality : int\n", + " Main frequency of the time series;\n", + " Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1.\n", + " train_df : pandas or polars DataFrame\n", + " Training dataframe with id and actual values. Must be sorted by time.\n", + " id_col : str (default='unique_id')\n", + " Column that identifies each serie.\n", + " target_col : str (default='y')\n", + " Column that contains the target.\n", + "\n", + " Returns\n", + " -------\n", + " pandas or polars Dataframe\n", + " dataframe with one row per id and one column per model.\n", + "\n", + " References\n", + " ----------\n", + " [1] https://otexts.com/fpp3/accuracy.html\n", + " \"\"\"\n", + " mean_sq_err = mse(df=df, models=models, id_col=id_col, target_col=target_col)\n", + " if isinstance(train_df, pd.DataFrame):\n", + " mean_sq_err = mean_sq_err.set_index(id_col)\n", + " # assume train_df is sorted\n", + " lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality)\n", + " scale = train_df[target_col].sub(lagged).pow(2)\n", + " scale = scale.groupby(train_df[id_col], observed=True).mean()\n", + " res = mean_sq_err.div(_zero_to_nan(scale), axis=0).fillna(0)\n", + " res.index.name = id_col\n", + " res = res.reset_index()\n", + " else:\n", + " # assume train_df is sorted\n", + " lagged = pl.col(target_col).shift(seasonality).over(id_col)\n", + " scale_expr = pl.col(target_col).sub(lagged).pow(2).alias(\"scale\")\n", + " scale = train_df.select([id_col, scale_expr])\n", + " scale = ufp.group_by(scale, id_col).mean()\n", + " scale = scale.with_columns(_zero_to_nan(pl.col(\"scale\")))\n", + "\n", + " def gen_expr(model):\n", + " return pl.col(model).truediv(pl.col(\"scale\")).fill_nan(0).alias(model)\n", + "\n", + " full_df = mean_sq_err.join(scale, on=id_col, how=\"left\")\n", + " res = _pl_agg_expr(full_df, models, id_col, gen_expr)\n", + " return res" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "---\n", + "\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L378){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "\n", + "#### msse\n", + "\n", + "> msse (df:~DFType, models:List[str], seasonality:int, train_df:~DFType,\n", + "> id_col:str='unique_id', target_col:str='y')\n", + "\n", + "*Mean Squared Scaled Error (MSSE)\n", + "\n", + "MSSE measures the relative prediction\n", + "accuracy of a forecasting method by comparinng the mean squared errors\n", + "of the prediction and the observed value against the mean\n", + "squared errors of the seasonal naive model.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, actuals and predictions. |\n", + "| models | List | | Columns that identify the models predictions. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ], + "text/plain": [ + "---\n", + "\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L378){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "\n", + "#### msse\n", + "\n", + "> msse (df:~DFType, models:List[str], seasonality:int, train_df:~DFType,\n", + "> id_col:str='unique_id', target_col:str='y')\n", + "\n", + "*Mean Squared Scaled Error (MSSE)\n", + "\n", + "MSSE measures the relative prediction\n", + "accuracy of a forecasting method by comparinng the mean squared errors\n", + "of the prediction and the observed value against the mean\n", + "squared errors of the seasonal naive model.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, actuals and predictions. |\n", + "| models | List | | Columns that identify the models predictions. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show_doc(msse, title_level=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| polars\n", + "pd_vs_pl(\n", + " msse(series, models, 7, series),\n", + " msse(series_pl, models, 7, series_pl),\n", + " models,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Root Mean Squared Scaled Error\n", + "\n", + "$$\n", + "\\mathrm{RMSSE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}_{\\tau}, \\mathbf{\\hat{y}}^{season}_{\\tau}) = \n", + "\\sqrt{\\frac{1}{H} \\sum^{t+H}_{\\tau=t+1} \\frac{(y_{\\tau}-\\hat{y}_{\\tau})^2}{\\mathrm{MSE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}^{season}_{\\tau})}}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| export\n", + "def rmsse(\n", + " df: DFType,\n", + " models: List[str],\n", + " seasonality: int,\n", + " train_df: DFType,\n", + " id_col: str = \"unique_id\",\n", + " target_col: str = \"y\",\n", + ") -> DFType:\n", + " res = msse(\n", + " df,\n", + " models=models,\n", + " seasonality=seasonality,\n", + " train_df=train_df,\n", + " id_col=id_col,\n", + " target_col=target_col,\n", + " )\n", + " if isinstance(res, pd.DataFrame):\n", + " res[models] = res[models].pow(0.5)\n", + " else:\n", + " res = res.with_columns(pl.col(models).pow(0.5))\n", + " return res\n", + "\n", + "rmsse.__doc__ = msse.__doc__.replace( # type: ignore[union-attr]\n", + " 'Mean Squared Scaled Error (MSSE)', 'Root Mean Squared Scaled Error (RMSSE)'\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "---\n", + "\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L444){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "\n", + "#### rmsse\n", + "\n", + "> rmsse (df:~DFType, models:List[str], seasonality:int, train_df:~DFType,\n", + "> id_col:str='unique_id', target_col:str='y')\n", + "\n", + "*Root Mean Squared Scaled Error (RMSSE)\n", + "\n", + "MSSE measures the relative prediction\n", + "accuracy of a forecasting method by comparinng the mean squared errors\n", + "of the prediction and the observed value against the mean\n", + "squared errors of the seasonal naive model.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, actuals and predictions. |\n", + "| models | List | | Columns that identify the models predictions. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ], + "text/plain": [ + "---\n", + "\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L444){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "\n", + "#### rmsse\n", + "\n", + "> rmsse (df:~DFType, models:List[str], seasonality:int, train_df:~DFType,\n", + "> id_col:str='unique_id', target_col:str='y')\n", + "\n", + "*Root Mean Squared Scaled Error (RMSSE)\n", + "\n", + "MSSE measures the relative prediction\n", + "accuracy of a forecasting method by comparinng the mean squared errors\n", + "of the prediction and the observed value against the mean\n", + "squared errors of the seasonal naive model.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, actuals and predictions. |\n", + "| models | List | | Columns that identify the models predictions. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show_doc(rmsse, title_level=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| polars\n", + "pd_vs_pl(\n", + " rmsse(series, models, 7, series),\n", + " rmsse(series_pl, models, 7, series_pl),\n", + " models,\n", + ")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1441,7 +1746,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L351){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L462){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### quantile_loss\n", "\n", @@ -1467,7 +1772,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L351){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L462){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### quantile_loss\n", "\n", @@ -1650,7 +1955,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L413){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L524){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mqloss\n", "\n", @@ -1682,7 +1987,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L413){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L524){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mqloss\n", "\n", @@ -1860,7 +2165,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L472){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L583){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### coverage\n", "\n", @@ -1881,7 +2186,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L472){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L583){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### coverage\n", "\n", @@ -1992,7 +2297,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L531){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L642){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### calibration\n", "\n", @@ -2012,7 +2317,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L531){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L642){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### calibration\n", "\n", @@ -2147,7 +2452,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L581){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L692){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### scaled_crps\n", "\n", @@ -2174,7 +2479,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L581){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L692){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### scaled_crps\n", "\n", diff --git a/utilsforecast/_modidx.py b/utilsforecast/_modidx.py index 452a512..736e430 100644 --- a/utilsforecast/_modidx.py +++ b/utilsforecast/_modidx.py @@ -79,9 +79,11 @@ 'utilsforecast.losses.mase': ('losses.html#mase', 'utilsforecast/losses.py'), 'utilsforecast.losses.mqloss': ('losses.html#mqloss', 'utilsforecast/losses.py'), 'utilsforecast.losses.mse': ('losses.html#mse', 'utilsforecast/losses.py'), + 'utilsforecast.losses.msse': ('losses.html#msse', 'utilsforecast/losses.py'), 'utilsforecast.losses.quantile_loss': ('losses.html#quantile_loss', 'utilsforecast/losses.py'), 'utilsforecast.losses.rmae': ('losses.html#rmae', 'utilsforecast/losses.py'), 'utilsforecast.losses.rmse': ('losses.html#rmse', 'utilsforecast/losses.py'), + 'utilsforecast.losses.rmsse': ('losses.html#rmsse', 'utilsforecast/losses.py'), 'utilsforecast.losses.scaled_crps': ('losses.html#scaled_crps', 'utilsforecast/losses.py'), 'utilsforecast.losses.smape': ('losses.html#smape', 'utilsforecast/losses.py')}, 'utilsforecast.plotting': { 'utilsforecast.plotting._filter_series': ( 'plotting.html#_filter_series', diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index 444cac2..fce47ca 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -1,8 +1,8 @@ # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/losses.ipynb. # %% auto 0 -__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', 'mqloss', 'coverage', 'calibration', - 'scaled_crps'] +__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', 'mqloss', 'coverage', + 'calibration', 'scaled_crps'] # %% ../nbs/losses.ipynb 3 from typing import Callable, Dict, List, Optional, Tuple, Union @@ -374,7 +374,7 @@ def gen_expr(model, baseline) -> pl_Expr: res = res.select([id_col, *exprs]) return res - +# %% ../nbs/losses.ipynb 53 def msse( df: DFType, models: List[str], @@ -415,7 +415,7 @@ def msse( ---------- [1] https://otexts.com/fpp3/accuracy.html """ - mean_sq_err = mse(df, models, id_col, target_col) + mean_sq_err = mse(df=df, models=models, id_col=id_col, target_col=target_col) if isinstance(train_df, pd.DataFrame): mean_sq_err = mean_sq_err.set_index(id_col) # assume train_df is sorted @@ -440,25 +440,35 @@ def gen_expr(model): res = _pl_agg_expr(full_df, models, id_col, gen_expr) return res - +# %% ../nbs/losses.ipynb 57 def rmsse( df: DFType, models: List[str], + seasonality: int, + train_df: DFType, id_col: str = "unique_id", target_col: str = "y", ) -> DFType: - """Root Mean Squared Scaled Error (RMSE) - - RMSSE is the square root of the MSSE""" - res = msse(df, models, id_col, target_col) + res = msse( + df, + models=models, + seasonality=seasonality, + train_df=train_df, + id_col=id_col, + target_col=target_col, + ) if isinstance(res, pd.DataFrame): res[models] = res[models].pow(0.5) else: - res = res.with_columns(*[pl.col(c).pow(0.5) for c in models]) + res = res.with_columns(pl.col(models).pow(0.5)) return res -# %% ../nbs/losses.ipynb 55 +rmsse.__doc__ = msse.__doc__.replace( # type: ignore[union-attr] + "Mean Squared Scaled Error (MSSE)", "Root Mean Squared Scaled Error (RMSSE)" +) + +# %% ../nbs/losses.ipynb 63 def quantile_loss( df: DFType, models: Dict[str, str], @@ -520,7 +530,7 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res -# %% ../nbs/losses.ipynb 61 +# %% ../nbs/losses.ipynb 69 def mqloss( df: DFType, models: Dict[str, List[str]], @@ -579,7 +589,7 @@ def mqloss( res = ufp.assign_columns(res, model, model_res[model]) return res -# %% ../nbs/losses.ipynb 67 +# %% ../nbs/losses.ipynb 75 def coverage( df: DFType, models: List[str], @@ -638,7 +648,7 @@ def gen_expr(model): res = _pl_agg_expr(df, models, id_col, gen_expr) return res -# %% ../nbs/losses.ipynb 71 +# %% ../nbs/losses.ipynb 79 def calibration( df: DFType, models: Dict[str, str], @@ -688,7 +698,7 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res -# %% ../nbs/losses.ipynb 75 +# %% ../nbs/losses.ipynb 83 def scaled_crps( df: DFType, models: Dict[str, List[str]], From 173db66b1754d1044df0e490940f4873fcb78cde Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Morales?= Date: Fri, 22 Nov 2024 11:07:28 -0600 Subject: [PATCH 3/8] use spawn for macos --- .github/workflows/ci.yaml | 4 +-- .pre-commit-config.yaml | 3 ++- action_files/nbdev_test | 52 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 56 insertions(+), 3 deletions(-) create mode 100755 action_files/nbdev_test diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 2dca57f..c35fd24 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -51,7 +51,7 @@ jobs: run: pip install uv && uv pip install --system ".[dev]" - name: Run tests - run: nbdev_test --do_print --timing --flags 'datasets matplotlib polars pyarrow scipy' + run: ./action_files/nbdev_test --do_print --timing --flags 'datasets matplotlib polars pyarrow scipy' minimal-tests: runs-on: ${{ matrix.os }} @@ -74,4 +74,4 @@ jobs: - name: Run tests shell: bash - run: nbdev_test --do_print --timing --skip_file_re 'plotting' + run: ./action_files/nbdev_test --do_print --timing --skip_file_re 'plotting' diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 64d80ec..599a5ec 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -10,9 +10,10 @@ repos: rev: v0.2.1 hooks: - id: ruff + files: utilsforecast - repo: https://github.com/pre-commit/mirrors-mypy rev: v1.8.0 hooks: - id: mypy args: [--ignore-missing-imports] - exclude: 'setup.py' + files: utilsforecast diff --git a/action_files/nbdev_test b/action_files/nbdev_test new file mode 100755 index 0000000..9cdda4b --- /dev/null +++ b/action_files/nbdev_test @@ -0,0 +1,52 @@ +#!/usr/bin/env python3 +import time,os,sys,traceback,contextlib, inspect +from fastcore.basics import * +from fastcore.imports import * +from fastcore.foundation import * +from fastcore.parallel import * +from fastcore.script import * +from fastcore.meta import delegates + +from nbdev.config import * +from nbdev.doclinks import * +from nbdev.process import NBProcessor, nb_lang +from nbdev.frontmatter import FrontmatterProc +from nbdev.test import _keep_file, test_nb + +from execnb.nbio import * +from execnb.shell import * + + +@call_parse +@delegates(nbglob_cli) +def nbdev_test( + path:str=None, # A notebook name or glob to test + flags:str='', # Space separated list of test flags to run that are normally ignored + n_workers:int=None, # Number of workers + timing:bool=False, # Time each notebook to see which are slow + do_print:bool=False, # Print start and end of each notebook + pause:float=0.01, # Pause time (in seconds) between notebooks to avoid race conditions + ignore_fname:str='.notest', # Filename that will result in siblings being ignored + **kwargs): + "Test in parallel notebooks matching `path`, passing along `flags`" + skip_flags = get_config().tst_flags.split() + force_flags = flags.split() + files = nbglob(path, as_path=True, **kwargs) + files = [f.absolute() for f in sorted(files) if _keep_file(f, ignore_fname)] + if len(files)==0: return print('No files were eligible for testing') + + if n_workers is None: n_workers = 0 if len(files)==1 else min(num_cpus(), 8) + kw = {'method': 'spawn'} + wd_pth = get_config().nbs_path + with working_directory(wd_pth if (wd_pth and wd_pth.exists()) else os.getcwd()): + results = parallel(test_nb, files, skip_flags=skip_flags, force_flags=force_flags, n_workers=n_workers, + basepath=get_config().config_path, pause=pause, do_print=do_print, **kw) + passed,times = zip(*results) + if all(passed): print("Success.") + else: + _fence = '='*50 + failed = '\n\t'.join(f.name for p,f in zip(passed,files) if not p) + sys.stderr.write(f"\nnbdev Tests Failed On The Following Notebooks:\n{_fence}\n\t{failed}\n") + sys.exit(1) + if timing: + for i,t in sorted(enumerate(times), key=lambda o:o[1], reverse=True): print(f"{files[i].name}: {int(t)} secs") From b25097b8cf1b7091ab3435de0e439344fa2fb9d4 Mon Sep 17 00:00:00 2001 From: Marco Zanotti Date: Wed, 29 Jan 2025 10:54:41 +0100 Subject: [PATCH 4/8] adding Scaled QL and Scaled MQL to losses.py --- utilsforecast/losses.py | 149 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 149 insertions(+) diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index 660c303..635b3ac 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -530,6 +530,78 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res +def scaled_quantile_loss( + df: DFType, + models: List[str], + q: float = 0.5, + seasonality: int, + train_df: DFType, + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """Scaled Quantile Loss (SQL) + + SQL measures the deviation of a quantile forecast scaled by + the mean absolute errors of the seasonal naive model. + By weighting the absolute deviation in a non symmetric way, the + loss pays more attention to under or over estimation. + A common value for q is 0.5 for the deviation from the median. + This was the official measure used in the M5 Uncertainty competition + with seasonality = 1. + + Parameters + ---------- + df : pandas or polars DataFrame + Input dataframe with id, times, actuals and predictions. + models : dict from str to str + Mapping from model name to the model predictions for the specified quantile. + q : float (default=0.5) + Quantile for the predictions' comparison. + seasonality : int + Main frequency of the time series; + Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. + train_df : pandas or polars DataFrame + Training dataframe with id and actual values. Must be sorted by time. + id_col : str (default='unique_id') + Column that identifies each serie. + target_col : str (default='y') + Column that contains the target. + + Returns + ------- + pandas or polars Dataframe + dataframe with one row per id and one column per model. + + References + ---------- + [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 + """ + q_loss = quantile_loss(df=df, models=models, q=q, id_col=id_col, target_col=target_col) + if isinstance(train_df, pd.DataFrame): + q_loss = q_loss.set_index(id_col) + # assume train_df is sorted + lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) + scale = train_df[target_col].sub(lagged).abs() + scale = scale.groupby(train_df[id_col], observed=True).mean() + res = q_loss.div(_zero_to_nan(scale), axis=0).fillna(0) + res.index.name = id_col + res = res.reset_index() + else: + # assume train_df is sorted + lagged = pl.col(target_col).shift(seasonality).over(id_col) + scale_expr = pl.col(target_col).sub(lagged).abs().alias("scale") + scale = train_df.select([id_col, scale_expr]) + scale = ufp.group_by(scale, id_col).mean() + scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) + + def gen_expr(model): + return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) + + full_df = q_loss.join(scale, on=id_col, how="left") + res = _pl_agg_expr(full_df, models, id_col, gen_expr) + return res + + # %% ../nbs/losses.ipynb 69 def mqloss( df: DFType, @@ -589,6 +661,83 @@ def mqloss( res = ufp.assign_columns(res, model, model_res[model]) return res +def scaled_mqloss( + df: DFType, + models: Dict[str, List[str]], + quantiles: np.ndarray, + seasonality: int, + train_df: DFType, + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """Scaled Multi-Quantile loss (SMQL) + + SMQL calculates the average multi-quantile Loss for + a given set of quantiles, based on the absolute + difference between predicted quantiles and observed values + scaled by the mean absolute errors of the seasonal naive model. + + The limit behavior of MQL allows to measure the accuracy + of a full predictive distribution with + the continuous ranked probability score (CRPS). This can be achieved + through a numerical integration technique, that discretizes the quantiles + and treats the CRPS integral with a left Riemann approximation, averaging over + uniformly distanced quantiles. + This was the official measure used in the M5 Uncertainty competition + with seasonality = 1. + + Parameters + ---------- + df : pandas or polars DataFrame + Input dataframe with id, times, actuals and predictions. + models : dict from str to list of str + Mapping from model name to the model predictions for each quantile. + quantiles : numpy array + Quantiles to compare against. + seasonality : int + Main frequency of the time series; + Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. + train_df : pandas or polars DataFrame + Training dataframe with id and actual values. Must be sorted by time. + id_col : str (default='unique_id') + Column that identifies each serie. + target_col : str (default='y') + Column that contains the target. + + Returns + ------- + pandas or polars Dataframe + dataframe with one row per id and one column per model. + + References + ---------- + [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 + """ + mq_loss = mqloss(df=df, models=models, quantiles=quantiles, id_col=id_col, target_col=target_col) + if isinstance(train_df, pd.DataFrame): + mq_loss = mq_loss.set_index(id_col) + # assume train_df is sorted + lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) + scale = train_df[target_col].sub(lagged).abs() + scale = scale.groupby(train_df[id_col], observed=True).mean() + res = mq_loss.div(_zero_to_nan(scale), axis=0).fillna(0) + res.index.name = id_col + res = res.reset_index() + else: + # assume train_df is sorted + lagged = pl.col(target_col).shift(seasonality).over(id_col) + scale_expr = pl.col(target_col).sub(lagged).abs().alias("scale") + scale = train_df.select([id_col, scale_expr]) + scale = ufp.group_by(scale, id_col).mean() + scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) + + def gen_expr(model): + return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) + + full_df = mq_loss.join(scale, on=id_col, how="left") + res = _pl_agg_expr(full_df, models, id_col, gen_expr) + return res + # %% ../nbs/losses.ipynb 75 def coverage( df: DFType, From b08e3a0b62641e5953d5f39c42798a4f1049f485 Mon Sep 17 00:00:00 2001 From: Marco Zanotti Date: Wed, 29 Jan 2025 10:56:16 +0100 Subject: [PATCH 5/8] adding Scaled QL and Scaled MQL to losses.py --- utilsforecast/losses.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index 635b3ac..3bd087e 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -1,8 +1,8 @@ # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/losses.ipynb. # %% auto 0 -__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', 'mqloss', 'coverage', - 'calibration', 'scaled_crps'] +__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', + 'quantile_loss', 'scaled_quantile_loss', 'mqloss', 'scaled_mqloss', 'coverage', 'calibration', 'scaled_crps'] # %% ../nbs/losses.ipynb 3 from typing import Callable, Dict, List, Optional, Tuple, Union From 40e5449282a1960d1a8d93de22fb4204f11a8957 Mon Sep 17 00:00:00 2001 From: Marco Zanotti Date: Fri, 21 Feb 2025 09:32:41 +0100 Subject: [PATCH 6/8] add _scale_loss function to be used in mase, msse, scaled_quantile_loss, and scaled_mql --- utilsforecast/losses.py | 166 ++++++++++++++++++---------------------- 1 file changed, 73 insertions(+), 93 deletions(-) diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index 3bd087e..28839fd 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -52,6 +52,75 @@ def _pl_agg_expr( df = df.select([id_col, *exprs]) return ufp.group_by(df, id_col, maintain_order=True).mean() +def _scale_loss( + loss_df: DFType, + scale_type: str = "absolute_error", + models: List[str], + seasonality: int, + train_df: DFType, + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """ + Parameters + ---------- + loss_df : pandas or polars DataFrame + Input dataframe with id, actuals, predictions and losses results. + scale_type : str (default='absolute_error') + Type of scaling. Possible values are 'absolute_error' or 'squared_error'. + models : list of str + Columns that identify the models predictions. + seasonality : int + Main frequency of the time series; + Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. + train_df : pandas or polars DataFrame + Training dataframe with id and actual values. Must be sorted by time. + id_col : str (default='unique_id') + Column that identifies each serie. + target_col : str (default='y') + Column that contains the target. + + Returns + ------- + pandas or polars Dataframe + dataframe with one row per id and one column per model. + + References + ---------- + [1] https://robjhyndman.com/papers/mase.pdf + """ + + if isinstance(train_df, pd.DataFrame): + loss_df = loss_df.set_index(id_col) + # assume train_df is sorted + lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) + if scale_type == "absolute_error": + scale = train_df[target_col].sub(lagged).abs() + else: + scale = train_df[target_col].sub(lagged).pow(2) + scale = scale.groupby(train_df[id_col], observed=True).mean() + res = loss_df.div(_zero_to_nan(scale), axis=0).fillna(0) + res.index.name = id_col + res = res.reset_index() + else: + # assume train_df is sorted + lagged = pl.col(target_col).shift(seasonality).over(id_col) + if scale_type == "absolute_error": + scale_expr = pl.col(target_col).sub(lagged).abs().alias("scale") + else: + scale_expr = pl.col(target_col).sub(lagged).pow(2).alias("scale") + scale = train_df.select([id_col, scale_expr]) + scale = ufp.group_by(scale, id_col).mean() + scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) + + def gen_expr(model): + return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) + + full_df = loss_df.join(scale, on=id_col, how="left") + res = _pl_agg_expr(full_df, models, id_col, gen_expr) + + return res + # %% ../nbs/losses.ipynb 13 @_base_docstring def mae( @@ -296,29 +365,7 @@ def mase( [1] https://robjhyndman.com/papers/mase.pdf """ mean_abs_err = mae(df, models, id_col, target_col) - if isinstance(train_df, pd.DataFrame): - mean_abs_err = mean_abs_err.set_index(id_col) - # assume train_df is sorted - lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) - scale = train_df[target_col].sub(lagged).abs() - scale = scale.groupby(train_df[id_col], observed=True).mean() - res = mean_abs_err.div(_zero_to_nan(scale), axis=0).fillna(0) - res.index.name = id_col - res = res.reset_index() - else: - # assume train_df is sorted - lagged = pl.col(target_col).shift(seasonality).over(id_col) - scale_expr = pl.col(target_col).sub(lagged).abs().alias("scale") - scale = train_df.select([id_col, scale_expr]) - scale = ufp.group_by(scale, id_col).mean() - scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) - - def gen_expr(model): - return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) - - full_df = mean_abs_err.join(scale, on=id_col, how="left") - res = _pl_agg_expr(full_df, models, id_col, gen_expr) - return res + return _scale_loss(mean_abs_err, "absolute_error", models, seasonality, train_df, id_col, target_col) # %% ../nbs/losses.ipynb 49 def rmae( @@ -416,29 +463,7 @@ def msse( [1] https://otexts.com/fpp3/accuracy.html """ mean_sq_err = mse(df=df, models=models, id_col=id_col, target_col=target_col) - if isinstance(train_df, pd.DataFrame): - mean_sq_err = mean_sq_err.set_index(id_col) - # assume train_df is sorted - lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) - scale = train_df[target_col].sub(lagged).pow(2) - scale = scale.groupby(train_df[id_col], observed=True).mean() - res = mean_sq_err.div(_zero_to_nan(scale), axis=0).fillna(0) - res.index.name = id_col - res = res.reset_index() - else: - # assume train_df is sorted - lagged = pl.col(target_col).shift(seasonality).over(id_col) - scale_expr = pl.col(target_col).sub(lagged).pow(2).alias("scale") - scale = train_df.select([id_col, scale_expr]) - scale = ufp.group_by(scale, id_col).mean() - scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) - - def gen_expr(model): - return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) - - full_df = mean_sq_err.join(scale, on=id_col, how="left") - res = _pl_agg_expr(full_df, models, id_col, gen_expr) - return res + return _scale_loss(mean_sq_err, "squared_error", models, seasonality, train_df, id_col, target_col) # %% ../nbs/losses.ipynb 57 def rmsse( @@ -577,30 +602,7 @@ def scaled_quantile_loss( [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 """ q_loss = quantile_loss(df=df, models=models, q=q, id_col=id_col, target_col=target_col) - if isinstance(train_df, pd.DataFrame): - q_loss = q_loss.set_index(id_col) - # assume train_df is sorted - lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) - scale = train_df[target_col].sub(lagged).abs() - scale = scale.groupby(train_df[id_col], observed=True).mean() - res = q_loss.div(_zero_to_nan(scale), axis=0).fillna(0) - res.index.name = id_col - res = res.reset_index() - else: - # assume train_df is sorted - lagged = pl.col(target_col).shift(seasonality).over(id_col) - scale_expr = pl.col(target_col).sub(lagged).abs().alias("scale") - scale = train_df.select([id_col, scale_expr]) - scale = ufp.group_by(scale, id_col).mean() - scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) - - def gen_expr(model): - return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) - - full_df = q_loss.join(scale, on=id_col, how="left") - res = _pl_agg_expr(full_df, models, id_col, gen_expr) - return res - + return _scale_loss(q_loss, "absolute_error", models, seasonality, train_df, id_col, target_col) # %% ../nbs/losses.ipynb 69 def mqloss( @@ -714,29 +716,7 @@ def scaled_mqloss( [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 """ mq_loss = mqloss(df=df, models=models, quantiles=quantiles, id_col=id_col, target_col=target_col) - if isinstance(train_df, pd.DataFrame): - mq_loss = mq_loss.set_index(id_col) - # assume train_df is sorted - lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality) - scale = train_df[target_col].sub(lagged).abs() - scale = scale.groupby(train_df[id_col], observed=True).mean() - res = mq_loss.div(_zero_to_nan(scale), axis=0).fillna(0) - res.index.name = id_col - res = res.reset_index() - else: - # assume train_df is sorted - lagged = pl.col(target_col).shift(seasonality).over(id_col) - scale_expr = pl.col(target_col).sub(lagged).abs().alias("scale") - scale = train_df.select([id_col, scale_expr]) - scale = ufp.group_by(scale, id_col).mean() - scale = scale.with_columns(_zero_to_nan(pl.col("scale"))) - - def gen_expr(model): - return pl.col(model).truediv(pl.col("scale")).fill_nan(0).alias(model) - - full_df = mq_loss.join(scale, on=id_col, how="left") - res = _pl_agg_expr(full_df, models, id_col, gen_expr) - return res + return _scale_loss(mq_loss, "absolute_error", models, seasonality, train_df, id_col, target_col) # %% ../nbs/losses.ipynb 75 def coverage( From 6ff0097da8845267af4ffd44ebf6a45ab3e22207 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Morales?= Date: Fri, 21 Feb 2025 16:37:14 -0600 Subject: [PATCH 7/8] update nb --- nbs/losses.ipynb | 525 ++++++++++++++++++++++++++++++++++----- utilsforecast/_modidx.py | 1 + utilsforecast/losses.py | 143 ++--------- 3 files changed, 494 insertions(+), 175 deletions(-) diff --git a/nbs/losses.ipynb b/nbs/losses.ipynb index ac0675c..d5a0e01 100644 --- a/nbs/losses.ipynb +++ b/nbs/losses.ipynb @@ -169,7 +169,75 @@ ") -> pl_DataFrame:\n", " exprs = [gen_expr(model) for model in models]\n", " df = df.select([id_col, *exprs])\n", - " return ufp.group_by(df, id_col, maintain_order=True).mean()" + " return ufp.group_by(df, id_col, maintain_order=True).mean()\n", + "\n", + "\n", + "def _scale_loss(\n", + " loss_df: DFType,\n", + " scale_type: str,\n", + " models: List[str],\n", + " seasonality: int,\n", + " train_df: DFType,\n", + " id_col: str = \"unique_id\",\n", + " target_col: str = \"y\",\n", + ") -> DFType:\n", + " \"\"\"\n", + " Parameters\n", + " ----------\n", + " loss_df : pandas or polars DataFrame\n", + " Input dataframe with id, actuals, predictions and losses results.\n", + " scale_type : str\n", + " Type of scaling. Possible values are 'absolute_error' or 'squared_error'.\n", + " models : list of str\n", + " Columns that identify the models predictions.\n", + " seasonality : int\n", + " Main frequency of the time series;\n", + " Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1.\n", + " train_df : pandas or polars DataFrame\n", + " Training dataframe with id and actual values. Must be sorted by time.\n", + " id_col : str (default='unique_id')\n", + " Column that identifies each serie.\n", + " target_col : str (default='y')\n", + " Column that contains the target.\n", + " Returns\n", + " -------\n", + " pandas or polars Dataframe\n", + " dataframe with one row per id and one column per model.\n", + " References\n", + " ----------\n", + " [1] https://robjhyndman.com/papers/mase.pdf\n", + " \"\"\"\n", + "\n", + " if isinstance(train_df, pd.DataFrame):\n", + " loss_df = loss_df.set_index(id_col)\n", + " # assume train_df is sorted\n", + " lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality)\n", + " if scale_type == \"absolute_error\":\n", + " scale = train_df[target_col].sub(lagged).abs()\n", + " else:\n", + " scale = train_df[target_col].sub(lagged).pow(2)\n", + " scale = scale.groupby(train_df[id_col], observed=True).mean()\n", + " res = loss_df.div(_zero_to_nan(scale), axis=0).fillna(0)\n", + " res.index.name = id_col\n", + " res = res.reset_index()\n", + " else:\n", + " # assume train_df is sorted\n", + " lagged = pl.col(target_col).shift(seasonality).over(id_col)\n", + " if scale_type == \"absolute_error\":\n", + " scale_expr = pl.col(target_col).sub(lagged).abs().alias(\"scale\")\n", + " else:\n", + " scale_expr = pl.col(target_col).sub(lagged).pow(2).alias(\"scale\")\n", + " scale = train_df.select([id_col, scale_expr])\n", + " scale = ufp.group_by(scale, id_col).mean()\n", + " scale = scale.with_columns(_zero_to_nan(pl.col(\"scale\")))\n", + "\n", + " def gen_expr(model):\n", + " return pl.col(model).truediv(pl.col(\"scale\")).fill_nan(0).alias(model)\n", + "\n", + " full_df = loss_df.join(scale, on=id_col, how=\"left\")\n", + " res = _pl_agg_expr(full_df, models, id_col, gen_expr)\n", + "\n", + " return res" ] }, { @@ -1068,29 +1136,15 @@ " [1] https://robjhyndman.com/papers/mase.pdf \n", " \"\"\"\n", " mean_abs_err = mae(df, models, id_col, target_col)\n", - " if isinstance(train_df, pd.DataFrame):\n", - " mean_abs_err = mean_abs_err.set_index(id_col)\n", - " # assume train_df is sorted\n", - " lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality)\n", - " scale = train_df[target_col].sub(lagged).abs()\n", - " scale = scale.groupby(train_df[id_col], observed=True).mean()\n", - " res = mean_abs_err.div(_zero_to_nan(scale), axis=0).fillna(0)\n", - " res.index.name = id_col\n", - " res = res.reset_index()\n", - " else:\n", - " # assume train_df is sorted\n", - " lagged = pl.col(target_col).shift(seasonality).over(id_col)\n", - " scale_expr = pl.col(target_col).sub(lagged).abs().alias('scale')\n", - " scale = train_df.select([id_col, scale_expr])\n", - " scale = ufp.group_by(scale, id_col).mean()\n", - " scale = scale.with_columns(_zero_to_nan(pl.col('scale')))\n", - "\n", - " def gen_expr(model):\n", - " return pl.col(model).truediv(pl.col('scale')).fill_nan(0).alias(model)\n", - "\n", - " full_df = mean_abs_err.join(scale, on=id_col, how='left')\n", - " res = _pl_agg_expr(full_df, models, id_col, gen_expr)\n", - " return res" + " return _scale_loss(\n", + " loss_df=mean_abs_err,\n", + " scale_type=\"absolute_error\",\n", + " models=models,\n", + " seasonality=seasonality,\n", + " train_df=train_df,\n", + " id_col=id_col,\n", + " target_col=target_col\n", + " )" ] }, { @@ -1399,29 +1453,15 @@ " [1] https://otexts.com/fpp3/accuracy.html\n", " \"\"\"\n", " mean_sq_err = mse(df=df, models=models, id_col=id_col, target_col=target_col)\n", - " if isinstance(train_df, pd.DataFrame):\n", - " mean_sq_err = mean_sq_err.set_index(id_col)\n", - " # assume train_df is sorted\n", - " lagged = train_df.groupby(id_col, observed=True)[target_col].shift(seasonality)\n", - " scale = train_df[target_col].sub(lagged).pow(2)\n", - " scale = scale.groupby(train_df[id_col], observed=True).mean()\n", - " res = mean_sq_err.div(_zero_to_nan(scale), axis=0).fillna(0)\n", - " res.index.name = id_col\n", - " res = res.reset_index()\n", - " else:\n", - " # assume train_df is sorted\n", - " lagged = pl.col(target_col).shift(seasonality).over(id_col)\n", - " scale_expr = pl.col(target_col).sub(lagged).pow(2).alias(\"scale\")\n", - " scale = train_df.select([id_col, scale_expr])\n", - " scale = ufp.group_by(scale, id_col).mean()\n", - " scale = scale.with_columns(_zero_to_nan(pl.col(\"scale\")))\n", - "\n", - " def gen_expr(model):\n", - " return pl.col(model).truediv(pl.col(\"scale\")).fill_nan(0).alias(model)\n", - "\n", - " full_df = mean_sq_err.join(scale, on=id_col, how=\"left\")\n", - " res = _pl_agg_expr(full_df, models, id_col, gen_expr)\n", - " return res" + " return _scale_loss(\n", + " loss_df=mean_sq_err,\n", + " scale_type=\"squared_error\",\n", + " models=models,\n", + " seasonality=seasonality,\n", + " train_df=train_df,\n", + " id_col=id_col,\n", + " target_col=target_col\n", + " )" ] }, { @@ -1746,7 +1786,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L462){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L472){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### quantile_loss\n", "\n", @@ -1772,7 +1812,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L462){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L472){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### quantile_loss\n", "\n", @@ -1859,6 +1899,188 @@ " )" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Scaled Quantile Loss\n", + "\n", + "$$\n", + "\\mathrm{SQL}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}^{(q)}_{\\tau}) = \n", + "\\frac{1}{H} \\sum^{t+H}_{\\tau=t+1} \n", + "\\frac{(1-q)\\,( \\hat{y}^{(q)}_{\\tau} - y_{\\tau} )_{+} \n", + "+ q\\,( y_{\\tau} - \\hat{y}^{(q)}_{\\tau} )_{+}}{\\mathrm{MAE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}^{season}_{\\tau})}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def scaled_quantile_loss(\n", + " df: DFType,\n", + " models: List[str],\n", + " seasonality: int,\n", + " train_df: DFType, \n", + " q: float = 0.5,\n", + " id_col: str = \"unique_id\",\n", + " target_col: str = \"y\",\n", + ") -> DFType:\n", + " \"\"\"Scaled Quantile Loss (SQL)\n", + "\n", + " SQL measures the deviation of a quantile forecast scaled by\n", + " the mean absolute errors of the seasonal naive model.\n", + " By weighting the absolute deviation in a non symmetric way, the\n", + " loss pays more attention to under or over estimation.\n", + " A common value for q is 0.5 for the deviation from the median.\n", + " This was the official measure used in the M5 Uncertainty competition\n", + " with seasonality = 1.\n", + "\n", + " Parameters\n", + " ----------\n", + " df : pandas or polars DataFrame\n", + " Input dataframe with id, times, actuals and predictions.\n", + " models : dict from str to str\n", + " Mapping from model name to the model predictions for the specified quantile.\n", + " seasonality : int\n", + " Main frequency of the time series;\n", + " Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1.\n", + " train_df : pandas or polars DataFrame\n", + " Training dataframe with id and actual values. Must be sorted by time. \n", + " q : float (default=0.5)\n", + " Quantile for the predictions' comparison.\n", + " id_col : str (default='unique_id')\n", + " Column that identifies each serie.\n", + " target_col : str (default='y')\n", + " Column that contains the target.\n", + "\n", + " Returns\n", + " -------\n", + " pandas or polars Dataframe\n", + " dataframe with one row per id and one column per model.\n", + "\n", + " References\n", + " ----------\n", + " [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722\n", + " \"\"\" \n", + " q_loss = quantile_loss(df=df, models=models, q=q, id_col=id_col, target_col=target_col)\n", + " return _scale_loss(\n", + " loss_df=q_loss,\n", + " scale_type=\"absolute_error\",\n", + " models=models,\n", + " seasonality=seasonality,\n", + " train_df=train_df,\n", + " id_col=id_col,\n", + " target_col=target_col,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "---\n", + "\n", + "### scaled_quantile_loss\n", + "\n", + "> scaled_quantile_loss (df:~DFType, models:List[str], seasonality:int,\n", + "> train_df:~DFType, q:float=0.5,\n", + "> id_col:str='unique_id', target_col:str='y')\n", + "\n", + "*Scaled Quantile Loss (SQL)\n", + "\n", + "SQL measures the deviation of a quantile forecast scaled by\n", + "the mean absolute errors of the seasonal naive model.\n", + "By weighting the absolute deviation in a non symmetric way, the\n", + "loss pays more attention to under or over estimation.\n", + "A common value for q is 0.5 for the deviation from the median.\n", + "This was the official measure used in the M5 Uncertainty competition\n", + "with seasonality = 1.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, times, actuals and predictions. |\n", + "| models | List | | Mapping from model name to the model predictions for the specified quantile. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| q | float | 0.5 | Quantile for the predictions' comparison. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ], + "text/plain": [ + "---\n", + "\n", + "### scaled_quantile_loss\n", + "\n", + "> scaled_quantile_loss (df:~DFType, models:List[str], seasonality:int,\n", + "> train_df:~DFType, q:float=0.5,\n", + "> id_col:str='unique_id', target_col:str='y')\n", + "\n", + "*Scaled Quantile Loss (SQL)\n", + "\n", + "SQL measures the deviation of a quantile forecast scaled by\n", + "the mean absolute errors of the seasonal naive model.\n", + "By weighting the absolute deviation in a non symmetric way, the\n", + "loss pays more attention to under or over estimation.\n", + "A common value for q is 0.5 for the deviation from the median.\n", + "This was the official measure used in the M5 Uncertainty competition\n", + "with seasonality = 1.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, times, actuals and predictions. |\n", + "| models | List | | Mapping from model name to the model predictions for the specified quantile. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| q | float | 0.5 | Quantile for the predictions' comparison. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show_doc(scaled_quantile_loss)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "#| polars\n", + "q_models = {\n", + " 0.1: {\n", + " 'model0': 'model0-lo-80',\n", + " 'model1': 'model1-lo-80',\n", + " },\n", + " 0.9: {\n", + " 'model0': 'model0-hi-80',\n", + " 'model1': 'model1-hi-80',\n", + " },\n", + "}\n", + "\n", + "for q in quantiles:\n", + " pd_vs_pl(\n", + " scaled_quantile_loss(series, q_models[q], seasonality=1, train_df=series, q=q),\n", + " scaled_quantile_loss(series_pl, q_models[q], seasonality=1, train_df=series_pl, q=q),\n", + " models,\n", + " )" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1955,7 +2177,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L524){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L534){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mqloss\n", "\n", @@ -1969,7 +2191,7 @@ "difference between predicted quantiles and observed values.\n", "\n", "The limit behavior of MQL allows to measure the accuracy \n", - "of a full predictive distribution $\\mathbf{\\hat{F}}_{\\tau}$ with \n", + "of a full predictive distribution with \n", "the continuous ranked probability score (CRPS). This can be achieved \n", "through a numerical integration technique, that discretizes the quantiles \n", "and treats the CRPS integral with a left Riemann approximation, averaging over \n", @@ -1987,7 +2209,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L524){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L534){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### mqloss\n", "\n", @@ -2001,7 +2223,7 @@ "difference between predicted quantiles and observed values.\n", "\n", "The limit behavior of MQL allows to measure the accuracy \n", - "of a full predictive distribution $\\mathbf{\\hat{F}}_{\\tau}$ with \n", + "of a full predictive distribution with \n", "the continuous ranked probability score (CRPS). This can be achieved \n", "through a numerical integration technique, that discretizes the quantiles \n", "and treats the CRPS integral with a left Riemann approximation, averaging over \n", @@ -2093,6 +2315,191 @@ " assert null_vals.item() == 0" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Scaled Multi-Quantile Loss\n", + "\n", + "$$\n", + "\\mathrm{MQL}(\\mathbf{y}_{\\tau},\n", + "[\\mathbf{\\hat{y}}^{(q_{1})}_{\\tau}, ... ,\\hat{y}^{(q_{n})}_{\\tau}]) = \n", + "\\frac{1}{n} \\sum_{q_{i}} \\frac{\\mathrm{QL}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}^{(q_{i})}_{\\tau})}{\\mathrm{MAE}(\\mathbf{y}_{\\tau}, \\mathbf{\\hat{y}}^{season}_{\\tau})}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def scaled_mqloss(\n", + " df: DFType,\n", + " models: Dict[str, List[str]],\n", + " quantiles: np.ndarray,\n", + " seasonality: int,\n", + " train_df: DFType,\n", + " id_col: str = \"unique_id\",\n", + " target_col: str = \"y\",\n", + ") -> DFType:\n", + " \"\"\"Scaled Multi-Quantile loss (SMQL)\n", + "\n", + " SMQL calculates the average multi-quantile Loss for\n", + " a given set of quantiles, based on the absolute\n", + " difference between predicted quantiles and observed values \n", + " scaled by the mean absolute errors of the seasonal naive model. \n", + " The limit behavior of MQL allows to measure the accuracy\n", + " of a full predictive distribution with\n", + " the continuous ranked probability score (CRPS). This can be achieved\n", + " through a numerical integration technique, that discretizes the quantiles\n", + " and treats the CRPS integral with a left Riemann approximation, averaging over\n", + " uniformly distanced quantiles.\n", + " This was the official measure used in the M5 Uncertainty competition\n", + " with seasonality = 1.\n", + "\n", + " Parameters\n", + " ----------\n", + " df : pandas or polars DataFrame\n", + " Input dataframe with id, times, actuals and predictions.\n", + " models : dict from str to list of str\n", + " Mapping from model name to the model predictions for each quantile.\n", + " quantiles : numpy array\n", + " Quantiles to compare against.\n", + " seasonality : int\n", + " Main frequency of the time series;\n", + " Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1.\n", + " train_df : pandas or polars DataFrame\n", + " Training dataframe with id and actual values. Must be sorted by time.\n", + " id_col : str (default='unique_id')\n", + " Column that identifies each serie.\n", + " target_col : str (default='y')\n", + " Column that contains the target.\n", + "\n", + " Returns\n", + " -------\n", + " pandas or polars Dataframe\n", + " dataframe with one row per id and one column per model.\n", + "\n", + " References\n", + " ----------\n", + " [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722\n", + " \"\"\"\n", + " mq_loss = mqloss(df=df, models=models, quantiles=quantiles, id_col=id_col, target_col=target_col)\n", + " return _scale_loss(\n", + " loss_df=mq_loss,\n", + " scale_type=\"absolute_error\",\n", + " models=models,\n", + " seasonality=seasonality,\n", + " train_df=train_df,\n", + " id_col=id_col,\n", + " target_col=target_col,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "data": { + "text/markdown": [ + "---\n", + "\n", + "### scaled_mqloss\n", + "\n", + "> scaled_mqloss (df:~DFType, models:Dict[str,List[str]],\n", + "> quantiles:numpy.ndarray, seasonality:int,\n", + "> train_df:~DFType, id_col:str='unique_id',\n", + "> target_col:str='y')\n", + "\n", + "*Scaled Multi-Quantile loss (SMQL)\n", + "\n", + "SMQL calculates the average multi-quantile Loss for\n", + "a given set of quantiles, based on the absolute\n", + "difference between predicted quantiles and observed values \n", + "scaled by the mean absolute errors of the seasonal naive model. \n", + "The limit behavior of MQL allows to measure the accuracy\n", + "of a full predictive distribution with\n", + "the continuous ranked probability score (CRPS). This can be achieved\n", + "through a numerical integration technique, that discretizes the quantiles\n", + "and treats the CRPS integral with a left Riemann approximation, averaging over\n", + "uniformly distanced quantiles.\n", + "This was the official measure used in the M5 Uncertainty competition\n", + "with seasonality = 1.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, times, actuals and predictions. |\n", + "| models | Dict | | Mapping from model name to the model predictions for each quantile. |\n", + "| quantiles | ndarray | | Quantiles to compare against. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ], + "text/plain": [ + "---\n", + "\n", + "### scaled_mqloss\n", + "\n", + "> scaled_mqloss (df:~DFType, models:Dict[str,List[str]],\n", + "> quantiles:numpy.ndarray, seasonality:int,\n", + "> train_df:~DFType, id_col:str='unique_id',\n", + "> target_col:str='y')\n", + "\n", + "*Scaled Multi-Quantile loss (SMQL)\n", + "\n", + "SMQL calculates the average multi-quantile Loss for\n", + "a given set of quantiles, based on the absolute\n", + "difference between predicted quantiles and observed values \n", + "scaled by the mean absolute errors of the seasonal naive model. \n", + "The limit behavior of MQL allows to measure the accuracy\n", + "of a full predictive distribution with\n", + "the continuous ranked probability score (CRPS). This can be achieved\n", + "through a numerical integration technique, that discretizes the quantiles\n", + "and treats the CRPS integral with a left Riemann approximation, averaging over\n", + "uniformly distanced quantiles.\n", + "This was the official measure used in the M5 Uncertainty competition\n", + "with seasonality = 1.*\n", + "\n", + "| | **Type** | **Default** | **Details** |\n", + "| -- | -------- | ----------- | ----------- |\n", + "| df | DFType | | Input dataframe with id, times, actuals and predictions. |\n", + "| models | Dict | | Mapping from model name to the model predictions for each quantile. |\n", + "| quantiles | ndarray | | Quantiles to compare against. |\n", + "| seasonality | int | | Main frequency of the time series;
Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. |\n", + "| train_df | DFType | | Training dataframe with id and actual values. Must be sorted by time. |\n", + "| id_col | str | unique_id | Column that identifies each serie. |\n", + "| target_col | str | y | Column that contains the target. |\n", + "| **Returns** | **DFType** | | **dataframe with one row per id and one column per model.** |" + ] + }, + "execution_count": null, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show_doc(scaled_mqloss)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| polars\n", + "pd_vs_pl(\n", + " scaled_mqloss(series, mq_models, quantiles=quantiles, seasonality=1, train_df=series),\n", + " scaled_mqloss(series_pl, mq_models, quantiles=quantiles, seasonality=1, train_df=series_pl),\n", + " models,\n", + ")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -2165,7 +2572,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L583){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L593){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### coverage\n", "\n", @@ -2186,7 +2593,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L583){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L593){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### coverage\n", "\n", @@ -2297,7 +2704,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L642){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L652){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### calibration\n", "\n", @@ -2317,7 +2724,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L642){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L652){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### calibration\n", "\n", @@ -2452,7 +2859,7 @@ "text/markdown": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L692){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L702){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### scaled_crps\n", "\n", @@ -2479,7 +2886,7 @@ "text/plain": [ "---\n", "\n", - "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L692){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", + "[source](https://github.com/Nixtla/utilsforecast/blob/main/utilsforecast/losses.py#L702){target=\"_blank\" style=\"float:right; font-size:smaller\"}\n", "\n", "#### scaled_crps\n", "\n", diff --git a/utilsforecast/_modidx.py b/utilsforecast/_modidx.py index 736e430..e1b1e0f 100644 --- a/utilsforecast/_modidx.py +++ b/utilsforecast/_modidx.py @@ -70,6 +70,7 @@ 'utilsforecast/grouped_array.py')}, 'utilsforecast.losses': { 'utilsforecast.losses._base_docstring': ('losses.html#_base_docstring', 'utilsforecast/losses.py'), 'utilsforecast.losses._pl_agg_expr': ('losses.html#_pl_agg_expr', 'utilsforecast/losses.py'), + 'utilsforecast.losses._scale_loss': ('losses.html#_scale_loss', 'utilsforecast/losses.py'), 'utilsforecast.losses._zero_to_nan': ('losses.html#_zero_to_nan', 'utilsforecast/losses.py'), 'utilsforecast.losses.bias': ('losses.html#bias', 'utilsforecast/losses.py'), 'utilsforecast.losses.calibration': ('losses.html#calibration', 'utilsforecast/losses.py'), diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index 28839fd..8538d77 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -1,8 +1,8 @@ # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/losses.ipynb. # %% auto 0 -__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', - 'quantile_loss', 'scaled_quantile_loss', 'mqloss', 'scaled_mqloss', 'coverage', 'calibration', 'scaled_crps'] +__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', 'mqloss', 'coverage', + 'calibration', 'scaled_crps'] # %% ../nbs/losses.ipynb 3 from typing import Callable, Dict, List, Optional, Tuple, Union @@ -52,9 +52,10 @@ def _pl_agg_expr( df = df.select([id_col, *exprs]) return ufp.group_by(df, id_col, maintain_order=True).mean() + def _scale_loss( loss_df: DFType, - scale_type: str = "absolute_error", + scale_type: str, models: List[str], seasonality: int, train_df: DFType, @@ -66,7 +67,7 @@ def _scale_loss( ---------- loss_df : pandas or polars DataFrame Input dataframe with id, actuals, predictions and losses results. - scale_type : str (default='absolute_error') + scale_type : str Type of scaling. Possible values are 'absolute_error' or 'squared_error'. models : list of str Columns that identify the models predictions. @@ -79,12 +80,10 @@ def _scale_loss( Column that identifies each serie. target_col : str (default='y') Column that contains the target. - Returns ------- pandas or polars Dataframe dataframe with one row per id and one column per model. - References ---------- [1] https://robjhyndman.com/papers/mase.pdf @@ -365,7 +364,15 @@ def mase( [1] https://robjhyndman.com/papers/mase.pdf """ mean_abs_err = mae(df, models, id_col, target_col) - return _scale_loss(mean_abs_err, "absolute_error", models, seasonality, train_df, id_col, target_col) + return _scale_loss( + loss_df=mean_abs_err, + scale_type="absolute_error", + models=models, + seasonality=seasonality, + train_df=train_df, + id_col=id_col, + target_col=target_col, + ) # %% ../nbs/losses.ipynb 49 def rmae( @@ -463,7 +470,15 @@ def msse( [1] https://otexts.com/fpp3/accuracy.html """ mean_sq_err = mse(df=df, models=models, id_col=id_col, target_col=target_col) - return _scale_loss(mean_sq_err, "squared_error", models, seasonality, train_df, id_col, target_col) + return _scale_loss( + loss_df=mean_sq_err, + scale_type="squared_error", + models=models, + seasonality=seasonality, + train_df=train_df, + id_col=id_col, + target_col=target_col, + ) # %% ../nbs/losses.ipynb 57 def rmsse( @@ -555,56 +570,7 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res -def scaled_quantile_loss( - df: DFType, - models: List[str], - q: float = 0.5, - seasonality: int, - train_df: DFType, - id_col: str = "unique_id", - target_col: str = "y", -) -> DFType: - """Scaled Quantile Loss (SQL) - - SQL measures the deviation of a quantile forecast scaled by - the mean absolute errors of the seasonal naive model. - By weighting the absolute deviation in a non symmetric way, the - loss pays more attention to under or over estimation. - A common value for q is 0.5 for the deviation from the median. - This was the official measure used in the M5 Uncertainty competition - with seasonality = 1. - - Parameters - ---------- - df : pandas or polars DataFrame - Input dataframe with id, times, actuals and predictions. - models : dict from str to str - Mapping from model name to the model predictions for the specified quantile. - q : float (default=0.5) - Quantile for the predictions' comparison. - seasonality : int - Main frequency of the time series; - Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. - train_df : pandas or polars DataFrame - Training dataframe with id and actual values. Must be sorted by time. - id_col : str (default='unique_id') - Column that identifies each serie. - target_col : str (default='y') - Column that contains the target. - - Returns - ------- - pandas or polars Dataframe - dataframe with one row per id and one column per model. - - References - ---------- - [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 - """ - q_loss = quantile_loss(df=df, models=models, q=q, id_col=id_col, target_col=target_col) - return _scale_loss(q_loss, "absolute_error", models, seasonality, train_df, id_col, target_col) - -# %% ../nbs/losses.ipynb 69 +# %% ../nbs/losses.ipynb 73 def mqloss( df: DFType, models: Dict[str, List[str]], @@ -663,62 +629,7 @@ def mqloss( res = ufp.assign_columns(res, model, model_res[model]) return res -def scaled_mqloss( - df: DFType, - models: Dict[str, List[str]], - quantiles: np.ndarray, - seasonality: int, - train_df: DFType, - id_col: str = "unique_id", - target_col: str = "y", -) -> DFType: - """Scaled Multi-Quantile loss (SMQL) - - SMQL calculates the average multi-quantile Loss for - a given set of quantiles, based on the absolute - difference between predicted quantiles and observed values - scaled by the mean absolute errors of the seasonal naive model. - - The limit behavior of MQL allows to measure the accuracy - of a full predictive distribution with - the continuous ranked probability score (CRPS). This can be achieved - through a numerical integration technique, that discretizes the quantiles - and treats the CRPS integral with a left Riemann approximation, averaging over - uniformly distanced quantiles. - This was the official measure used in the M5 Uncertainty competition - with seasonality = 1. - - Parameters - ---------- - df : pandas or polars DataFrame - Input dataframe with id, times, actuals and predictions. - models : dict from str to list of str - Mapping from model name to the model predictions for each quantile. - quantiles : numpy array - Quantiles to compare against. - seasonality : int - Main frequency of the time series; - Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. - train_df : pandas or polars DataFrame - Training dataframe with id and actual values. Must be sorted by time. - id_col : str (default='unique_id') - Column that identifies each serie. - target_col : str (default='y') - Column that contains the target. - - Returns - ------- - pandas or polars Dataframe - dataframe with one row per id and one column per model. - - References - ---------- - [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 - """ - mq_loss = mqloss(df=df, models=models, quantiles=quantiles, id_col=id_col, target_col=target_col) - return _scale_loss(mq_loss, "absolute_error", models, seasonality, train_df, id_col, target_col) - -# %% ../nbs/losses.ipynb 75 +# %% ../nbs/losses.ipynb 83 def coverage( df: DFType, models: List[str], @@ -777,7 +688,7 @@ def gen_expr(model): res = _pl_agg_expr(df, models, id_col, gen_expr) return res -# %% ../nbs/losses.ipynb 79 +# %% ../nbs/losses.ipynb 87 def calibration( df: DFType, models: Dict[str, str], @@ -827,7 +738,7 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res -# %% ../nbs/losses.ipynb 83 +# %% ../nbs/losses.ipynb 91 def scaled_crps( df: DFType, models: Dict[str, List[str]], From fcc7327aec0c4611ece86b44ec75c4cbf6ca3924 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20Morales?= Date: Fri, 21 Feb 2025 16:42:30 -0600 Subject: [PATCH 8/8] export --- nbs/losses.ipynb | 8 ++- utilsforecast/_modidx.py | 3 + utilsforecast/losses.py | 129 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 135 insertions(+), 5 deletions(-) diff --git a/nbs/losses.ipynb b/nbs/losses.ipynb index d5a0e01..058456a 100644 --- a/nbs/losses.ipynb +++ b/nbs/losses.ipynb @@ -1919,9 +1919,10 @@ "metadata": {}, "outputs": [], "source": [ + "#| export\n", "def scaled_quantile_loss(\n", " df: DFType,\n", - " models: List[str],\n", + " models: Dict[str, str],\n", " seasonality: int,\n", " train_df: DFType, \n", " q: float = 0.5,\n", @@ -1969,7 +1970,7 @@ " return _scale_loss(\n", " loss_df=q_loss,\n", " scale_type=\"absolute_error\",\n", - " models=models,\n", + " models=list(models.keys()),\n", " seasonality=seasonality,\n", " train_df=train_df,\n", " id_col=id_col,\n", @@ -2334,6 +2335,7 @@ "metadata": {}, "outputs": [], "source": [ + "#| export\n", "def scaled_mqloss(\n", " df: DFType,\n", " models: Dict[str, List[str]],\n", @@ -2389,7 +2391,7 @@ " return _scale_loss(\n", " loss_df=mq_loss,\n", " scale_type=\"absolute_error\",\n", - " models=models,\n", + " models=list(models.keys()),\n", " seasonality=seasonality,\n", " train_df=train_df,\n", " id_col=id_col,\n", diff --git a/utilsforecast/_modidx.py b/utilsforecast/_modidx.py index e1b1e0f..a1b205a 100644 --- a/utilsforecast/_modidx.py +++ b/utilsforecast/_modidx.py @@ -86,6 +86,9 @@ 'utilsforecast.losses.rmse': ('losses.html#rmse', 'utilsforecast/losses.py'), 'utilsforecast.losses.rmsse': ('losses.html#rmsse', 'utilsforecast/losses.py'), 'utilsforecast.losses.scaled_crps': ('losses.html#scaled_crps', 'utilsforecast/losses.py'), + 'utilsforecast.losses.scaled_mqloss': ('losses.html#scaled_mqloss', 'utilsforecast/losses.py'), + 'utilsforecast.losses.scaled_quantile_loss': ( 'losses.html#scaled_quantile_loss', + 'utilsforecast/losses.py'), 'utilsforecast.losses.smape': ('losses.html#smape', 'utilsforecast/losses.py')}, 'utilsforecast.plotting': { 'utilsforecast.plotting._filter_series': ( 'plotting.html#_filter_series', 'utilsforecast/plotting.py'), diff --git a/utilsforecast/losses.py b/utilsforecast/losses.py index 8538d77..e98662b 100644 --- a/utilsforecast/losses.py +++ b/utilsforecast/losses.py @@ -1,8 +1,8 @@ # AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/losses.ipynb. # %% auto 0 -__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', 'mqloss', 'coverage', - 'calibration', 'scaled_crps'] +__all__ = ['mae', 'mse', 'rmse', 'bias', 'mape', 'smape', 'mase', 'rmae', 'msse', 'rmsse', 'quantile_loss', + 'scaled_quantile_loss', 'mqloss', 'scaled_mqloss', 'coverage', 'calibration', 'scaled_crps'] # %% ../nbs/losses.ipynb 3 from typing import Callable, Dict, List, Optional, Tuple, Union @@ -570,6 +570,66 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res +# %% ../nbs/losses.ipynb 68 +def scaled_quantile_loss( + df: DFType, + models: Dict[str, str], + seasonality: int, + train_df: DFType, + q: float = 0.5, + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """Scaled Quantile Loss (SQL) + + SQL measures the deviation of a quantile forecast scaled by + the mean absolute errors of the seasonal naive model. + By weighting the absolute deviation in a non symmetric way, the + loss pays more attention to under or over estimation. + A common value for q is 0.5 for the deviation from the median. + This was the official measure used in the M5 Uncertainty competition + with seasonality = 1. + + Parameters + ---------- + df : pandas or polars DataFrame + Input dataframe with id, times, actuals and predictions. + models : dict from str to str + Mapping from model name to the model predictions for the specified quantile. + seasonality : int + Main frequency of the time series; + Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. + train_df : pandas or polars DataFrame + Training dataframe with id and actual values. Must be sorted by time. + q : float (default=0.5) + Quantile for the predictions' comparison. + id_col : str (default='unique_id') + Column that identifies each serie. + target_col : str (default='y') + Column that contains the target. + + Returns + ------- + pandas or polars Dataframe + dataframe with one row per id and one column per model. + + References + ---------- + [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 + """ + q_loss = quantile_loss( + df=df, models=models, q=q, id_col=id_col, target_col=target_col + ) + return _scale_loss( + loss_df=q_loss, + scale_type="absolute_error", + models=list(models.keys()), + seasonality=seasonality, + train_df=train_df, + id_col=id_col, + target_col=target_col, + ) + # %% ../nbs/losses.ipynb 73 def mqloss( df: DFType, @@ -629,6 +689,71 @@ def mqloss( res = ufp.assign_columns(res, model, model_res[model]) return res +# %% ../nbs/losses.ipynb 79 +def scaled_mqloss( + df: DFType, + models: Dict[str, List[str]], + quantiles: np.ndarray, + seasonality: int, + train_df: DFType, + id_col: str = "unique_id", + target_col: str = "y", +) -> DFType: + """Scaled Multi-Quantile loss (SMQL) + + SMQL calculates the average multi-quantile Loss for + a given set of quantiles, based on the absolute + difference between predicted quantiles and observed values + scaled by the mean absolute errors of the seasonal naive model. + The limit behavior of MQL allows to measure the accuracy + of a full predictive distribution with + the continuous ranked probability score (CRPS). This can be achieved + through a numerical integration technique, that discretizes the quantiles + and treats the CRPS integral with a left Riemann approximation, averaging over + uniformly distanced quantiles. + This was the official measure used in the M5 Uncertainty competition + with seasonality = 1. + + Parameters + ---------- + df : pandas or polars DataFrame + Input dataframe with id, times, actuals and predictions. + models : dict from str to list of str + Mapping from model name to the model predictions for each quantile. + quantiles : numpy array + Quantiles to compare against. + seasonality : int + Main frequency of the time series; + Hourly 24, Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1. + train_df : pandas or polars DataFrame + Training dataframe with id and actual values. Must be sorted by time. + id_col : str (default='unique_id') + Column that identifies each serie. + target_col : str (default='y') + Column that contains the target. + + Returns + ------- + pandas or polars Dataframe + dataframe with one row per id and one column per model. + + References + ---------- + [1] https://www.sciencedirect.com/science/article/pii/S0169207021001722 + """ + mq_loss = mqloss( + df=df, models=models, quantiles=quantiles, id_col=id_col, target_col=target_col + ) + return _scale_loss( + loss_df=mq_loss, + scale_type="absolute_error", + models=list(models.keys()), + seasonality=seasonality, + train_df=train_df, + id_col=id_col, + target_col=target_col, + ) + # %% ../nbs/losses.ipynb 83 def coverage( df: DFType,