diff --git a/tutorials/Binary_Contingency_Scores.ipynb b/tutorials/Binary_Contingency_Scores.ipynb index 234accd5..8e4fe16b 100644 --- a/tutorials/Binary_Contingency_Scores.ipynb +++ b/tutorials/Binary_Contingency_Scores.ipynb @@ -15,16 +15,16 @@ { "cell_type": "markdown", "id": "7e3b4f1c-7e82-4b29-9353-3c52eae9fa1e", - "metadata": { - "jp-MarkdownHeadingCollapsed": true - }, + "metadata": {}, "source": [ - "## Scores for Binary Categorical Forecasts\n", + "# Scores for Binary Categorical Forecasts\n", "\n", "Examples of binary categorical scores include: hit rate (true positive rate), false alarm rate and accuracy.\n", "\n", "These scores are calculated based on a binary contingency table, which comprises true positives (hits), true negatives (correct negatives), false positives (false alarms) and false negatives (misses).\n", "\n", + "If you would like to know more, visit the [Methods for dichotomous (yes/no) forecasts](https://www.cawcr.gov.au/projects/verification/#Methods_for_dichotomous_forecasts) information produced by the WWRP/WGNE Joint Working Group on Forecast Verification Research.\n", + "\n", "Producing these categorical statistics is inherently a three-step process.\n", "\n", "1. For data that is either continuous or categorical, first convert it to binary data using an \"event definition\". For data that is already binary, go straight to step two.\n", @@ -41,11 +41,11 @@ " - threat score (critical success index)\n", " - the Peirce skill score (true skill statistic, Hanssen and Kuipers discriminant)\n", " - specificity (true negative rate)\n", - " - f1 score\n", + " - F1 score\n", " - equitable threat score (Gilbert's skill score)\n", " - Heidke skill score (Cohen's kappa)\n", " - odds ratio\n", - " - odds ratio skill score (Yule's q)\n", + " - odds ratio skill score (Yule's Q)\n", "\n", "\n", "`scores` provides support for all three steps, plus convenient functions for working efficiently. Most of the `scores` APIs are functional in nature, however introducing some classes (see below) makes binary categorical metrics much easier to calculate and work with. This notebook starts with simple examples and builds up to more complex ones." @@ -86,7 +86,7 @@ "\n", "Here is an example of using `scores` to convert non-binary data into binary data using the \"event operator\" class.\n", "\n", - "Consider the following two idealized tables of forecast and observed values (in mm per hour) for precipitation." + "Consider the following two idealized tables of forecast and observed values (in mm) for precipitation." ] }, { @@ -100,17 +100,17 @@ "simple_forecast = xr.DataArray(\n", " [\n", "\t\t[\n", - "\t\t\t[0.9, 0.0, 5], \n", - "\t\t\t[0.7, 1.4, 2.8],\n", - "\t\t\t[.4, 0.5, 2.3],\n", + "\t\t\t[35.9, 44.0, 49], \n", + "\t\t\t[50.7, 51.4, 52.8],\n", + "\t\t\t[38.4, 38.5, 22.3],\n", "\t\t], \n", "\t\t\t[\n", - "\t\t\t[1.9, 1.0, 1.5], \n", - "\t\t\t[1.7, 2.4, 1.1],\n", - "\t\t\t[1.4, 1.5, 3.3],\n", + "\t\t\t[25.6, 38.1, 38.2], \n", + "\t\t\t[52.4, 55.9, 51.0],\n", + "\t\t\t[29.1, 33.1, 24.5],\n", "\t\t], \n", "\t],\n", - "\tcoords=[[10, 20], [0, 1, 2], [5, 6, 7]], dims=[\"height\", \"lat\", \"lon\"])" + "\tcoords=[[10, 11], [0, 1, 2], [5, 6, 7]], dims=[\"lead_hour\", \"lat\", \"lon\"])" ] }, { @@ -120,22 +120,22 @@ "metadata": {}, "outputs": [], "source": [ - "# Idealized observations are within 0.1 or 0.2 of the forecast in all cases except one\n", - "# Can be used to find some exact matches, and some close matches\n", + "# Idealized observations are within 1.0 or 4.0 of the forecast in all cases except one\n", + "# This can be used to find some exact matches, and some close matches\n", "simple_obs = xr.DataArray(\n", " [\n", "\t\t[\n", - "\t\t\t[0.9, 0.0, 5], \n", - "\t\t\t[0.7, 1.3, 2.7],\n", - "\t\t\t[.3, 0.4, 2.2],\n", + "\t\t\t[35.8, 44.2, 52], \n", + "\t\t\t[53.7, 52.3, 53.7],\n", + "\t\t\t[38.3, 36.4, 21.4],\n", "\t\t], \n", "\t\t\t[\n", - "\t\t\t[1.7, 1.2, 1.7], \n", - "\t\t\t[1.7, 2.2, 3.9],\n", - "\t\t\t[1.6, 1.2, 9.9],\n", + "\t\t\t[25.7, 38.2, 38.7], \n", + "\t\t\t[51.7, 56.2, 49.9],\n", + "\t\t\t[29.6, 30.2, 16.9],\n", "\t\t], \n", "\t],\n", - "\tcoords=[[10, 20], [0, 1, 2], [5, 6, 7]], dims=[\"height\", \"lat\", \"lon\"])" + "\tcoords=[[10, 11], [0, 1, 2], [5, 6, 7]], dims=[\"lead_hour\", \"lat\", \"lon\"])" ] }, { @@ -143,11 +143,11 @@ "id": "f3b9a4ad-3a63-45ee-90f6-bf2087aa18cf", "metadata": {}, "source": [ - "For the purposes of this example, let's say an event occurs if the precipitation value is greater than 1.3. We therefore define a threshold based event operator with a threshold of 1.3.\n", + "For the purposes of this example, let's say an event occurs if the precipitation value reaches 50. We therefore define a threshold based event operator with a threshold of 50.\n", "\n", - "Applying this operator means `scores` converts the non-binary data into binary data. It does this by creating binary event data - precipitation data with a value less than or equal to 1.3 is \"false\" and precipitation data with a value greater than 1.3 is \"true\".\n", + "Applying this operator means `scores` converts the non-binary data into binary data. It does this by creating binary event data - precipitation data with a value less than 50 is \"false\" and precipitation data with a value greater than or equal to 50 is \"true\".\n", "\n", - "The following code creates a \"greater than\" operator - using `scores` built-in ThresholdEventOperator - which creates events with a threshold of \"> 1.3\". It can be repeatedly called with different event thresholds, for reasons that will be explained later.\n" + "The following code creates a \"greater than or equal to\" operator - using `scores` built-in ThresholdEventOperator - which creates events with a threshold of \">= 50\". It can be repeatedly called with different event thresholds, for reasons that will be explained later.\n" ] }, { @@ -157,21 +157,23 @@ "metadata": {}, "outputs": [], "source": [ - "# An event here is defined as precipitation with a value greater than 1.3\n", + "# An event here is defined as precipitation with a value greater than 50\n", "\n", "# The EventThresholdOperator can take a variety of operators from the python \"operator\" module, or a user-defined function\n", - "# The default is operator.gt, which is the same as \">\" but in functional form.\n", + "# The default is operator.ge, which is the same as \">=\" but in functional form.\n", + "# Here we pass it in explicitly so it's very clear what will happen\n", "\n", - "event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=1.3)\n", + "event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=50.0, default_op_fn=operator.ge)\n", + "# event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=50.0) # This is the same thing\n", "\n", "# https://docs.python.org/3/library/operator.html provides an overview of what standard operator types can be supplied\n", "# You may wish to exeriment with the following alternatives now\n", "\n", "# Greater-than (not equal to) operator\n", - "# event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=1.3, default_op_fn=operator.gt)\n", + "# event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=50, default_op_fn=operator.gt)\n", "\n", "# Less-than operator\n", - "# event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=1.3, default_op_fn=operator.lt)" + "# event_operator = scores.categorical.ThresholdEventOperator(default_event_threshold=50, default_op_fn=operator.lt)" ] }, { @@ -184,18 +186,18 @@ "name": "stdout", "output_type": "stream", "text": [ - " Size: 144B\n", - "array([[[0., 0., 1.],\n", - " [0., 1., 1.],\n", - " [0., 0., 1.]],\n", + " Size: 144B\n", + "array([[[0., 0., 0.],\n", + " [1., 1., 1.],\n", + " [0., 0., 0.]],\n", "\n", - " [[1., 0., 1.],\n", - " [1., 1., 0.],\n", - " [1., 1., 1.]]])\n", + " [[0., 0., 0.],\n", + " [1., 1., 1.],\n", + " [0., 0., 0.]]])\n", "Coordinates:\n", - " * height (height) int64 16B 10 20\n", - " * lat (lat) int64 24B 0 1 2\n", - " * lon (lon) int64 24B 5 6 7\n" + " * lead_hour (lead_hour) int64 16B 10 11\n", + " * lat (lat) int64 24B 0 1 2\n", + " * lon (lon) int64 24B 5 6 7\n" ] } ], @@ -207,11 +209,41 @@ "forecast_binary, observed_binary = event_operator.make_event_tables(simple_forecast, simple_obs)\n", "\n", "# It is also possible to vary both the event threshold and the event operator from the defaults during this step\n", - "# forecast_binary, observed_binary = event_operator.make_event_tables(simple_forecast, simple_obs, event_threshold=2.0, op_fn=operator.lt)\n", + "# This commented-out code uses the less-than operator with a threshold of 45\n", + "# forecast_binary, observed_binary = event_operator.make_event_tables(simple_forecast, simple_obs, event_threshold=45, op_fn=operator.lt)\n", "\n", "print(forecast_binary)" ] }, + { + "cell_type": "code", + "execution_count": 6, + "id": "d7fccf1b-f55c-442f-a265-78022d853ca8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Size: 144B\n", + "array([[[0., 0., 1.],\n", + " [1., 1., 1.],\n", + " [0., 0., 0.]],\n", + "\n", + " [[0., 0., 0.],\n", + " [1., 1., 0.],\n", + " [0., 0., 0.]]])\n", + "Coordinates:\n", + " * lead_hour (lead_hour) int64 16B 10 11\n", + " * lat (lat) int64 24B 0 1 2\n", + " * lon (lon) int64 24B 5 6 7\n" + ] + } + ], + "source": [ + "print(observed_binary)" + ] + }, { "cell_type": "markdown", "id": "dec28e21-6d2c-486c-8387-b6d164dce570", @@ -224,7 +256,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "id": "1cab2d73-ca4b-4c90-b0e2-db940c0cd198", "metadata": {}, "outputs": [ @@ -595,19 +627,19 @@ " fill: currentColor;\n", "}\n", "
<xarray.DataArray (contingency: 5)> Size: 40B\n",
-       "array([10.,  6.,  1.,  1., 18.])\n",
+       "array([ 5., 11.,  1.,  1., 18.])\n",
        "Coordinates:\n",
-       "  * contingency  (contingency) <U11 220B 'tp_count' 'tn_count' ... 'total_count'
" + " * contingency (contingency) <U11 220B 'tp_count' 'tn_count' ... 'total_count'" ], "text/plain": [ " Size: 40B\n", - "array([10., 6., 1., 1., 18.])\n", + "array([ 5., 11., 1., 1., 18.])\n", "Coordinates:\n", " * contingency (contingency) 
<xarray.DataArray (contingency: 5)> Size: 40B\n",
-       "array([ 9.,  6.,  2.,  1., 18.])\n",
+       "array([ 5., 11.,  1.,  1., 18.])\n",
        "Coordinates:\n",
-       "  * contingency  (contingency) <U11 220B 'tp_count' 'tn_count' ... 'total_count'
" + " * contingency (contingency) <U11 220B 'tp_count' 'tn_count' ... 'total_count'" ], "text/plain": [ " Size: 40B\n", - "array([ 9., 6., 2., 1., 18.])\n", + "array([ 5., 11., 1., 1., 18.])\n", "Coordinates:\n", " * contingency (contingency) 
<xarray.DataArray (height: 2, lat: 3, lon: 3)> Size: 144B\n",
-       "array([[[0., 0., 1.],\n",
-       "        [0., 1., 1.],\n",
-       "        [0., 0., 1.]],\n",
+       "
<xarray.DataArray (lead_hour: 2, lat: 3, lon: 3)> Size: 144B\n",
+       "array([[[0., 0., 0.],\n",
+       "        [1., 1., 1.],\n",
+       "        [0., 0., 0.]],\n",
        "\n",
-       "       [[1., 0., 1.],\n",
-       "        [1., 1., 0.],\n",
-       "        [1., 1., 1.]]])\n",
+       "       [[0., 0., 0.],\n",
+       "        [1., 1., 1.],\n",
+       "        [0., 0., 0.]]])\n",
        "Coordinates:\n",
-       "  * height   (height) int64 16B 10 20\n",
-       "  * lat      (lat) int64 24B 0 1 2\n",
-       "  * lon      (lon) int64 24B 5 6 7
"
+       "  * lead_hour  (lead_hour) int64 16B 10 11\n",
+       "  * lat        (lat) int64 24B 0 1 2\n",
+       "  * lon        (lon) int64 24B 5 6 7
" ], "text/plain": [ - " Size: 144B\n", - "array([[[0., 0., 1.],\n", - " [0., 1., 1.],\n", - " [0., 0., 1.]],\n", - "\n", - " [[1., 0., 1.],\n", - " [1., 1., 0.],\n", - " [1., 1., 1.]]])\n", + " Size: 144B\n", + "array([[[0., 0., 0.],\n", + " [1., 1., 1.],\n", + " [0., 0., 0.]],\n", + "\n", + " [[0., 0., 0.],\n", + " [1., 1., 1.],\n", + " [0., 0., 0.]]])\n", "Coordinates:\n", - " * height (height) int64 16B 10 20\n", - " * lat (lat) int64 24B 0 1 2\n", - " * lon (lon) int64 24B 5 6 7" + " * lead_hour (lead_hour) int64 16B 10 11\n", + " * lat (lat) int64 24B 0 1 2\n", + " * lon (lon) int64 24B 5 6 7" ] }, - "execution_count": 9, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } @@ -1456,13 +1488,13 @@ "\n", "Below are two examples of using `scores` to calculate metrics directly from the contigency manager:\n", "\n", - "- \"accuracy\": (true positives + false negatives)/(total number of samples)\n", - "- \"false_alarm_rate\": (number of false positives)/(total number of samples)" + "- \"accuracy\": (true positives + true negatives)/(total number of samples)\n", + "- \"false_alarm_rate\": (number of false positives)/(true negatives + false negatives)" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 11, "id": "f69c5863-57f9-4b85-a3bd-c2d81163d63a", "metadata": {}, "outputs": [ @@ -1833,14 +1865,14 @@ " fill: currentColor;\n", "}\n", "
<xarray.DataArray ()> Size: 8B\n",
-       "array(0.83333333)
" + "array(0.88888889)" ], "text/plain": [ " Size: 8B\n", - "array(0.83333333)" + "array(0.88888889)" ] }, - "execution_count": 10, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -1851,7 +1883,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 12, "id": "2b2ae506-5e5f-4934-bb7f-509ca955b36c", "metadata": {}, "outputs": [ @@ -2222,14 +2254,14 @@ " fill: currentColor;\n", "}\n", "
<xarray.DataArray ()> Size: 8B\n",
-       "array(0.25)
" + "array(0.08333333)" ], "text/plain": [ " Size: 8B\n", - "array(0.25)" + "array(0.08333333)" ] }, - "execution_count": 11, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -2254,7 +2286,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, "id": "336d9f57-8fb7-440b-bcf0-21cdb862bf8e", "metadata": {}, "outputs": [ @@ -2625,14 +2657,14 @@ " fill: currentColor;\n", "}\n", "
<xarray.DataArray ()> Size: 8B\n",
-       "array(0.83333333)
" + "array(0.88888889)" ], "text/plain": [ " Size: 8B\n", - "array(0.83333333)" + "array(0.88888889)" ] }, - "execution_count": 12, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -2658,7 +2690,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 14, "id": "b6da1026-9e82-4dd5-ba7e-f66b14d31ed4", "metadata": {}, "outputs": [ @@ -2667,27 +2699,27 @@ "output_type": "stream", "text": [ "Contingency Manager (xarray view of table):\n", - " Size: 80B\n", - "array([[3., 6.],\n", - " [5., 1.],\n", - " [1., 1.],\n", + " Size: 80B\n", + "array([[3., 2.],\n", + " [5., 6.],\n", " [0., 1.],\n", + " [1., 0.],\n", " [9., 9.]])\n", "Coordinates:\n", - " * height (height) int64 16B 10 20\n", + " * lead_hour (lead_hour) int64 16B 10 11\n", " * contingency (contingency) 
<xarray.DataArray (height: 2)> Size: 16B\n",
-       "array([0.88888889, 0.77777778])\n",
+       "
<xarray.DataArray (lead_hour: 2)> Size: 16B\n",
+       "array([0.88888889, 0.88888889])\n",
        "Coordinates:\n",
-       "  * height   (height) int64 16B 10 20
"
+       "  * lead_hour  (lead_hour) int64 16B 10 11
" ], "text/plain": [ - " Size: 16B\n", - "array([0.88888889, 0.77777778])\n", + " Size: 16B\n", + "array([0.88888889, 0.88888889])\n", "Coordinates:\n", - " * height (height) int64 16B 10 20" + " * lead_hour (lead_hour) int64 16B 10 11" ] }, - "execution_count": 14, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } @@ -3091,6 +3123,14 @@ "\n", " - Calculating multiple contingency tables for different thresholds to observe how accuracy changes according to different event thresholds" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e68d29fb-9259-463c-8fa5-a9fa7391ab5d", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -3109,7 +3149,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.8" + "version": "3.12.2" } }, "nbformat": 4,