added tutorial

tutorials/Risk_Matrix_Score.ipynb

@@ -0,0 +1,129 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# The Risk Matrix Score\n",
+    "\n",
+    "The Risk Matrix Score provides a consistent way of scoring forecasts and warnings within a risk matrix framework (Taggart & Wilke 2024).\n",
+    "\n",
+    "In this tutorial we will illustrate the risk matrix framework using a slimmed down version of the synthetic case study from Taggart & Wilke (2024).\n",
+    "Consider a simple warning service for hazardous outdoor heat with warning levels \"Nil\", \"Yellow Warning\", \"Orange Warning\" and \"Red Warning\". The warning level that is issued depends on a combination of certainty and severity, as depicted in the matrix below.\n",
+    "\n",
+    "![risk_matrix](images/risk_matrix.png)\n",
+    "\n",
+    "The certainty categories are defined by probability thresholds (0.1, 0.3, 0.5) and the three nested severity categories (MOD+, SEV+ and EXT) are defined by three temperature thresholds (35, 37, 40) in degrees Celsius. A forecaster makes a risk assessment by selecting one certainty category for each of the nested severity categories, according to their probability distribution of daily maximum temperature. This is equivalent to selecting one cell in each of the three right-most columns in the risk matrix. The warning level is determined by the high colour of the selected cells. For example, if a forecaster selects \"likely\" for MOD+ (yellow), \"possible\" for SEV+ (yellow) and \"possible\" for EXT (orange) then this generates an Orange Warning.\n",
+    "\n",
+    "## Risk matrix scores from a synthetic experiment\n",
+    "\n",
+    "We now set a synthetic experiment which will generate forecasts from three different forecasters and corresponding observations.\n",
+    "\n",
+    "Let $N(\\mu,\\sigma^2)$ denote the normal distribution with mean $\\mu$ and variance $\\sigma^2$. Suppose that each observed daily maximum\n",
+    "temperature $y$ is generated from random variables $y_1$, $y_2$ and $y_3$ using the formula $y = y_1 +y_2 +y_3$, where each random\n",
+    "variable $y_i$ is independent of the others, $y_1 \\sim N(20,10^2)$, $y_2 \\sim N(0,5^2)$ and $y_3 \\sim N(0,2^2)$. The random variables $y_1$, $y_2$\n",
+    "and $y_3$ can be thought to represent variability from seasonal, synoptic scale and mesoscale weather processes respectively. The\n",
+    "climatic distribution of maximum temperature is $N(20,10^2+5^2+2^2)$. The climatic probability of exceeding 35°C, 37°C and\n",
+    "40°C on any given randomly selected day is therefore 0.093, 0.067 and 0.039.\n",
+    "\n",
+    "We now consider three forecasters named NeverWarnNate, SeasonalSam and SynopticSally. NeverWarnNate only knows about the climatic distribution $N(20,10^2+5^2+2^2)$ and uses this to fill out the risk matrix. SeasonalSam has access to seasonal information $y_1$ but not synoptic or mesoscale information forecasts. He uses the probability distribution $N(y1,5^2 +2^2)$ to fill out the risk matrix. SynopticSally has access to seasonal information $y_1$ and synoptic information $y_2$ but not mesoscale information. She uses the probability distribution $N(y_1 +y_2,2^2)$ to fill out the risk matrix.\n",
+    "\n",
+    "This experiment is repeated independently 1000 times.\n",
+    "\n",
+    "We will now implement this in python.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "from scipy.stats import norm\n",
+    "\n",
+    "# the number of forecast cases\n",
+    "n_fcst_cases = 1000\n",
+    "\n",
+    "# temperature thresholds\n",
+    "temp_thresholds = [35, 37, 40]\n",
+    "\n",
+    "# probability thresholds\n",
+    "prob_thresholds = [0.1, 0.3, 0.5]\n",
+    "\n",
+    "# severity coords\n",
+    "sev_coords = [\"MOD+\", \"SEV+\", \"EXT\"]\n",
+    "\n",
+    "# generate y1, y2 and y3\n",
+    "y1 = xr.DataArray(\n",
+    "    data=norm.rvs(loc=20, scale=10, size=n_fcst_cases),\n",
+    "    dims=['fcst_case'],\n",
+    "    coords={'fcst_case': np.arange(n_fcst_cases)}\n",
+    ")\n",
+    "y2 = xr.DataArray(\n",
+    "    data=norm.rvs(loc=0, scale=5, size=n_fcst_cases),\n",
+    "    dims=['fcst_case'],\n",
+    "    coords={'fcst_case': np.arange(n_fcst_cases)}\n",
+    ")\n",
+    "y3 = xr.DataArray(\n",
+    "    data=norm.rvs(loc=0, scale=2, size=n_fcst_cases),\n",
+    "    dims=['fcst_case'],\n",
+    "    coords={'fcst_case': np.arange(n_fcst_cases)}\n",
+    ")\n",
+    "\n",
+    "# calculate the observations\n",
+    "obs = y1 + y2 + y3\n",
+    "\n",
+    "# generate broadcast arrays for the random variables\n",
+    "da_temp_thresholds = xr.DataArray(\n",
+    "    data=temp_thresholds,\n",
+    "    dims=[\"severity\"],\n",
+    "    coords={\"severity\": sev_coords},\n",
+    ")\n",
+    "y1, y2, y3, da_temp_thresholds = xr.broadcast(y1, y2, y3, da_temp_thresholds)\n",
+    "\n",
+    "# calculate SynopticSally probability of exceedance forecasts\n",
+    "sally = xr.full_like(y1, np.nan)\n",
+    "sally.values = 1 - norm.cdf(da_temp_thresholds, loc=y1 + y2, scale=2)\n",
+    "\n",
+    "# calculate SeasonalSam's probability of exceedance forecasts\n",
+    "sam = xr.full_like(y1, np.nan)\n",
+    "sam.values = 1 - norm.cdf(da_temp_thresholds, loc=y1, scale=np.sqrt(5 ** 2 + 2 ** 2))\n",
+    "\n",
+    "# calculate NeverWarnNate's probability of exceedance forecasts\n",
+    "nate = xr.full_like(y1, 20)\n",
+    "nate.values = 1 - norm.cdf(da_temp_thresholds, loc=nate, scale=np.sqrt(10 ** 2 + 5 ** 2 + 2 ** 2))\n",
+    "\n",
+    "# we'll combine the forecasts into one data set\n",
+    "fcsts = xr.Dataset({\n",
+    "    \"SynopticSally\": sally,\n",
+    "    \"SeasonalSam\": sam,\n",
+    "    \"NeverWarnNate\": nate,\n",
+    "})"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pymc_env",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.6"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}