diff --git a/nbs/losses.ipynb b/nbs/losses.ipynb index ac0675c..058456a 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,189 @@ " )" ] }, + { + "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": [ + "#| export\n", + "def scaled_quantile_loss(\n", + " df: DFType,\n", + " models: Dict[str, 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=list(models.keys()),\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 +2178,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 +2192,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 +2210,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 +2224,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 +2316,192 @@ " 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": [ + "#| export\n", + "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=list(models.keys()),\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 +2574,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 +2595,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 +2706,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 +2726,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 +2861,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 +2888,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..a1b205a 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'), @@ -85,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 660c303..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 @@ -52,6 +52,74 @@ 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, + 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 + 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 +364,15 @@ 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( + 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( @@ -416,29 +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) - 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( + 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( @@ -530,7 +570,67 @@ def gen_expr(model): res = _pl_agg_expr(df, list(models.items()), id_col, gen_expr) return res -# %% ../nbs/losses.ipynb 69 +# %% ../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, models: Dict[str, List[str]], @@ -589,7 +689,72 @@ def mqloss( res = ufp.assign_columns(res, model, model_res[model]) return res -# %% ../nbs/losses.ipynb 75 +# %% ../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, models: List[str], @@ -648,7 +813,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], @@ -698,7 +863,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]],