01-objective-benchmark-demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2e372c41-beb6-49c3-a00e-8bd938fa354a",
   "metadata": {},
   "source": [
    "# How to select an objective function using information theory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fbc99a9-9a6c-48ff-bfcf-4c5c30b3a04d",
   "metadata": {},
   "source": [
    " [![GitHub tag (latest by date)](https://img.shields.io/github/v/tag/hytest-org/workflow-hodson-2022-objective-benchmark)](https://github.com)\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/hytest-org/workflow-hodson-2022-objective-benchmark/blob/main/01-objective-benchmark-demo.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e35c615c-4dfa-425c-bc6a-6a1ed0c4482a",
   "metadata": {},
   "source": [
    "## Abstract\n",
    "Science tests competing theories or models by evaluating the similarity of their predictions against observational experience. \n",
    "Thus, how we measure similarity fundamentally determines what we learn.\n",
    "In machine learning and scientific modeling, similarity metrics are used as objective functions.\n",
    "A classic example being mean squared error, which is the optimal measure of similarity when errors are normally distributed and independent and identically distributed (iid). \n",
    "In many cases, however, the error distribution is neither normal nor iid, so it is left to the scientist to determine an appropriate objective.\n",
    "Here, we review how information theory can guide that selection, then demonstrate the approach with a simple hydrologic model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8923947-ed9f-4e24-9ff7-dc14f04fd32d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Introduction\n",
    "Science seeks to create useful representations of reality in the form of hypotheses, models, or theories.\n",
    "What sets science apart from other pursuits is that it rigorously tests those representations against observational experience,\n",
    "favoring those that best fit the evidence.\n",
    "An analogous process occurs when calibrating a numerical model or evaluating among competing models.\n",
    "To select the ``best'' model, experiment by varying the model while keeping the test data and objective fixed.\n",
    "If mean squared error (MSE) is the objective, compute the MSE between the test data and the model predictions, then select the model with the lowest MSE.\n",
    "But why choose MSE and not another objective function?\n",
    "The answer: MSE is the optimal measure when errors are normally distributed and iid.\n",
    "But for many problems, the true error distribution is complex or unknown.\n",
    "Rather than simply assuming some de facto objective function, \n",
    "compare them against the evidence.\n",
    "This paper demonstrates how.\n",
    "\n",
    "To select the best objective, the experiment is essentially the same except the objective is varied while the model and data are held fixed.\n",
    "Now, select the objective indicating the greatest similarity between data and model.\n",
    "Different objective functions have different scales, so they are normalized\n",
    "such that each integrates to one, thereby representing them as probability distributions.\n",
    "The normalized form of MSE is the normal distribution, for example (Hodson, 2022).\n",
    "When used to evaluate model fit, the probability distribution is called a likelihood function\n",
    "and its output the likelihood.\n",
    "So, to select among objectives, compare their likelihoods, and favor the most likely.\n",
    "Taking the natural logarithm of the likelihood, denoted as $\\ell$,\n",
    "does not change the models ranks\n",
    "but simplifies the math by converting products to sums:\n",
    "likelihoods multiply, so log likelihoods add.\n",
    "\n",
    "\n",
    "So, to benchmark objectives, compare their log likelihoods.\n",
    "The maximum likelihood estimators for the log likelihoods of several objectives are implemented below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ef72a63-4417-4cb0-beea-260bb1b80c42",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# compute likelihood\n",
    "import numpy as np\n",
    "from scipy.stats import pearsonr\n",
    "import scipy.stats\n",
    "\n",
    "\n",
    "def normal_ll(y, y_hat, transform=None, gradient=1):\n",
    "    '''Log likelihood for the normal distribution with change of variable\n",
    "    \n",
    "    The normal distribution is the formal likelihood for the mean squared error (MSE).\n",
    "    \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "        Observations.\n",
    "    y_hat : array_like\n",
    "        Predictions.\n",
    "    transform : function\n",
    "        Change of variable transformation.\n",
    "    gradient : function\n",
    "        Gradient of the transform function.\n",
    "        \n",
    "    Proof\n",
    "    -----\n",
    "    https://www.statlect.com/probability-distributions/normal-distribution\n",
    "    '''\n",
    "    if transform is not None:\n",
    "        y = transform(y)\n",
    "        y_hat = transform(y_hat)\n",
    "        \n",
    "    e = y - y_hat\n",
    "    n = len(e)\n",
    "    sigma = e.std()\n",
    "    log_gradient = np.sum(np.log(np.abs(gradient)))\n",
    "    ll = -n * np.log(sigma) - n/2*np.log(2*np.pi) - 1/(2*sigma**2) * (e**2).sum() + log_gradient\n",
    "    return ll\n",
    "\n",
    "\n",
    "def laplace_ll(y, y_hat, transform=None, gradient=1):\n",
    "    '''Log likelihood for Laplace distribution with change of variable\n",
    "    \n",
    "    The laplace distribution is the formal likelihood for the mean absolute\n",
    "    error (MAE).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "        Observations.\n",
    "    y_hat : array_like\n",
    "        Predictions.\n",
    "    transform : function\n",
    "        Change of variable transformation.\n",
    "    gradient : function\n",
    "        Gradient of the transform function.\n",
    "    '''\n",
    "    if transform is not None:\n",
    "        y = transform(y)\n",
    "        y_hat = transform(y_hat)\n",
    "        \n",
    "    e = (y - y_hat).abs()\n",
    "    n = len(e)\n",
    "    b = e.mean()\n",
    "    log_gradient = np.sum(np.log(np.abs(gradient)))\n",
    "    ll = -n * np.log(2*b) - 1/b * e.sum() + log_gradient\n",
    "    return ll.sum()\n",
    "                                   \n",
    "\n",
    "def msre_ll(y, y_hat):\n",
    "    '''Log likelihood for mean squared square-root error\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    return normal_ll(y, y_hat, transform=lambda x: np.sqrt(x), gradient=-1/(2*np.sqrt(y)))\n",
    "\n",
    "\n",
    "def mare_ll(y, y_hat):\n",
    "    '''Log likelihood for mean absolute square-root error\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    return laplace_ll(y, y_hat, transform=lambda x: np.sqrt(x), gradient=-1/(2*np.sqrt(y)))\n",
    "\n",
    "\n",
    "def lognormal_ll(y, y_hat):\n",
    "    '''Lognormal log likelihood\n",
    "    \n",
    "    The lognormal distribution is the formal likelihood for the mean squared\n",
    "    log error (MSLE).\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    return normal_ll(y, y_hat, transform=lambda x: np.log(x), gradient=1/y)\n",
    "\n",
    "\n",
    "def mspe_ll(y, y_hat):\n",
    "    '''Log likelhood for mean squared percentage error\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    \n",
    "    '''\n",
    "    return normal_ll(y, y_hat, transform=lambda x: x/y, gradient=-1/(y**2)) \n",
    "\n",
    "\n",
    "def nse_ll(y, y_hat, group='gage_id'):\n",
    "    '''Log likelihood for normalized squared error (NSE)\n",
    "    \n",
    "    NSE is equivalent to the Nash–Sutcliffe model efficiency coefficient.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    sigma_o = y.groupby('gage_id').transform(lambda x: x.std())\n",
    "    return normal_ll(y, y_hat, transform=lambda x: x/sigma_o, gradient=1/sigma_o)\n",
    "\n",
    "\n",
    "def loglaplace_ll(y, y_hat):\n",
    "    '''Log likelihood for log Laplace distribution\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    return laplace_ll(y, y_hat, transform=lambda x: np.log(x), gradient=1/y)\n",
    "\n",
    "\n",
    "def uniform_ll(y, y_hat):\n",
    "    '''Log likelihood for uniform distribution.\n",
    "    \n",
    "    The uniform log likelihood minimizes the maximum error.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    e = np.abs(y - y_hat)\n",
    "    n = len(e)\n",
    "    #ll = -n * np.log(e.max()-e.min()) # standard formulation\n",
    "    ll = -n * np.log(e.max() - 0)\n",
    "    return ll\n",
    "\n",
    "\n",
    "def bernoulli_ll(y, y_hat, groupby=None):\n",
    "    '''TODO and use within zi_ll\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    pass\n",
    "\n",
    "\n",
    "\n",
    "def zi_ll(y, y_hat, ll=normal_ll, threshold=0.01, groupby=None):\n",
    "    ''' Zero-inflated log likelihood.\n",
    "    \n",
    "     Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    ll : function\n",
    "        Zero-inflated log likelihood \n",
    "    threshold : float\n",
    "        Value below which is treated as zero\n",
    "    groupby : string\n",
    "        Optional groupby term (testing)\n",
    "    '''\n",
    "    y_o = y <= threshold\n",
    "    y_hat_o = y_hat <= threshold\n",
    "    \n",
    "    if groupby is None:\n",
    "        n1 = (y_o & y_hat_o).sum() # correct zero-flow prediction\n",
    "        n2 = (y_o ^ y_hat_o).sum() # incorrect zero-flow prediction \n",
    "    else:\n",
    "        n1 = (y_o & y_hat_o).groupby(groupby).sum() # correct zero-flow prediction\n",
    "        n2 = (y_o ^ y_hat_o).groupby(groupby).sum() # incorrect zero-flow prediction\n",
    "\n",
    "    n3 = (~y_o & ~y_hat_o) # correct flow predictions\n",
    "    \n",
    "    # fraction of correctly predicted zero flows\n",
    "    rho = np.where( (n1+n2) == 0, 0, n1 / (n1 + n2))\n",
    "    n_rho = 1-rho\n",
    "    \n",
    "    # n1 * np.log(rho) + n2 * np.log(1-rho)\n",
    "    ll_zero = n1[rho!=0] * np.log(rho[rho!=0]) + n2[n_rho!=0]* np.log(n_rho[n_rho!=0])\n",
    "    \n",
    "    return ll_zero.sum() + ll(y[n3], y_hat[n3])\n",
    "\n",
    "\n",
    "def zilognormal_ll(y, y_hat):\n",
    "    '''Log likelihood for zero-inflated lognormal.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "       \n",
    "    return zi_ll(y, y_hat, ll=lognormal_ll, threshold=0.01)\n",
    "\n",
    "\n",
    "def ziloglaplace_ll(y, y_hat):\n",
    "    '''Log likelihood for zero-inflated laplace.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    y : array_like\n",
    "    y_hat : array_like\n",
    "    '''\n",
    "    return zi_ll(y, y_hat, ll=loglaplace_ll, threshold=0.01)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "165d83af-4345-4235-9ff8-b0a417a418b6",
   "metadata": {},
   "source": [
    "## Weights\n",
    "Given a set of $m$ models,\n",
    "the \"weight\" of evidence for each model $w_i$ is\n",
    "$$\n",
    "w_i = \\frac{ x^{-\\hat H_i} }{ \\sum^{m}_{i} x^{-\\hat H_i}  }\n",
    "$$\n",
    "where the base $x$ is 2 for bits or $e$ for nats (adapted from Burnham and Anderson, 2002)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0adaa4c8-de06-489f-a0fb-4d2b01c9c3a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_weights(series, base=np.e):\n",
    "    '''Compute posterior weights\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    series : array_like\n",
    "        Log likelihoods\n",
    "    base: float\n",
    "        Base of the logarithm used to compute log likelihood\n",
    "    '''\n",
    "    s = base**series\n",
    "    return s/s.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c6eb8b6-99f9-419b-a1b0-29e80b756ac5",
   "metadata": {},
   "source": [
    "## Benchmark demonstration\n",
    "To demonstrate, we benchmark the entropies of several objective functions that might be considered for a streamflow model.\n",
    "The test data are streamflow observations from 1,385 streamgages in the conterminous United States \\citep{Russell_2020};\n",
    "roughly 14 million observations.\n",
    "As streamflow can be zero or negative, which is undefined for some objective functions,\n",
    "flows below 0.0028 m$^3$ s$^{-1}$ (0.01 ft$^3$ s$^{-1}$) were thresholded and treated as the ``zero-flow'' state in the comparison.\n",
    "Different thresholds may yield slightly different results,\n",
    "particularly among logged objectives because of their greater sensitivity near zero.\n",
    "\n",
    "The model is simple: \n",
    "predict streamflow at a location by scaling the nearest concurrent observation by the ratio of the two drainage areas.\n",
    "So when predicting flow in a large river using observations from a smaller one, scale up the observations accordingly.\n",
    "By nature, the predictions are out of sample, so neither cross validation nor bias adjustment is necessary. \n",
    "We chose this example because it represents the case\n",
    "in which the model is physically correct, but its boundary conditions are uncertain (a common problem in Earth science).\n",
    "An alternate experiment takes a general model, \n",
    "like a neural network or a physical simulation,\n",
    "then calibrates it to each objective.\n",
    "The former experiment tests how well the objectives represent uncertainty about the model input,\n",
    "whereas the latter tests how well they represent uncertainty about the model structure as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcf028c0-c5fe-4033-9f66-e5af7f0ba7db",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data from s3 (run once)\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import fsspec\n",
    "\n",
    "fs_read = fsspec.filesystem('s3', anon=True, client_kwargs={'endpoint_url': \"https://usgs.osn.mghpcc.org/\"})\n",
    "\n",
    "with fs_read.open('s3://hytest/nwis_gages2_streamflow_estimates/gages2_nndar.parquet') as f:\n",
    "    df = pd.read_parquet(f)\n",
    "    \n",
    "# save local copy\n",
    "df.to_parquet('gages2_nndar.parquet')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "858c72e9-c5d9-4a1a-ab24-6be7e6c37d9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read local copy\n",
    "df = pd.read_parquet('gages2_nndar.parquet')\n",
    "df[df < 0.01] = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d17876b7-aec1-48e6-86d5-94c6ce788768",
   "metadata": {},
   "outputs": [],
   "source": [
    "# step 1: create a table of objective functions\n",
    "objectives = {\n",
    "    'U' : {'name':'uniformly distributed error', 'f':uniform_ll, 'k':1},\n",
    "    'MSE' : {'name':'mean squared error', 'f':normal_ll, 'k':1},\n",
    "    'NSE' : {'name':'normalized squared error', 'f':nse_ll, 'k':1},\n",
    "    'MAE' : {'name': 'mean absolute error', 'f':laplace_ll, 'k':1},\n",
    "    'MSPE' : {'name': 'mean squared percent error', 'f':mspe_ll, 'k':1},\n",
    "    'MSLE' : {'name':'mean squared log error*', 'f':lognormal_ll, 'k':1},\n",
    "    'MALE' : {'name':'mean absolute log error*', 'f':loglaplace_ll, 'k':1},\n",
    "    'ZMSLE' : {'name':'zero-inflated MSLE', 'f':zilognormal_ll, 'k':2},\n",
    "    'ZMALE' : {'name':'zero-inflated MALE', 'f':ziloglaplace_ll, 'k':2},\n",
    "    'MARE' : {'name':'mean absolute square root error', 'f':mare_ll, 'k':1},\n",
    "}\n",
    "\n",
    "obj_df = pd.DataFrame.from_dict(objectives, orient='index')\n",
    "\n",
    "# step 2: compute the information in each objective function\n",
    "for index, row in obj_df.iterrows():\n",
    "    # nats is the negative log likelihood or the info in the error\n",
    "    obj_df.loc[index, 'bits'] = (row.k - row.f(df.obs, df.NNDAR))/len(df)/np.log(2)\n",
    "\n",
    "# step 3: compute weights\n",
    "obj_df['weight'] = compute_weights(-obj_df.bits, base=2) #this was a negative\n",
    "\n",
    "# step 4: format output table\n",
    "\n",
    "table = obj_df[['name', 'k', 'bits', 'weight']].sort_values('weight').round(2)#.rename(columns=names)\n",
    "\n",
    "table['rank'] = len(table) - np.argsort(table['weight'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f19532f6-cf72-4adf-96fa-d98aa34d5f21",
   "metadata": {},
   "source": [
    "Results are shown in the table below.\n",
    "In the experiment, the data and model were fixed;\n",
    "therefore, so was the information in the error.\n",
    "All that changed was how we measured it.\n",
    "Relative to ZMALE, the excess bits in the other objective functions are noise.\n",
    "So, MSE measures at least 40 percent noise,\n",
    "and NSE at least 38 percent.\n",
    "Consider stochastic gradient descent, where noisier gradients require more iterations to reach the solution.\n",
    "In that case, each iteration completes faster, so the solution may be reached quicker overall.\n",
    "A poorly chosen objective incurs a similar penalty but potentially without benefit.\n",
    "In general, noisier objectives convey less information and so require more iterations during calibration,\n",
    "more data to reach the solution,\n",
    "and produce models that require more storage space\n",
    "(better model, better data compression)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a920adc-2381-4d72-bc8b-54024a412374",
   "metadata": {},
   "outputs": [],
   "source": [
    "#print(table.to_latex())\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "227556b9-d9e1-48f7-af8b-0028f0c830e2",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "Despite their ubiquitous use as a basis for learning,\n",
    "objective functions are rarely benchmarked, except by Bayesians.\n",
    "We do not advocate for one over another\n",
    "---the choice varies from problem to problem---\n",
    "only that benchmarking objectives is good practice that will yield better models,\n",
    "both in terms of their accuracy and in how they are evaluated:\n",
    "objective functions quantify model performance,\n",
    "as well as uncertainty and potential compression.\n",
    "Perhaps given this simple demonstration, more scientists will test different measures of similarity by benchmarking their objectives.\n",
    "After all, how well machines---and scientists---learn and think depends a lot on how well they measure similarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "def8d189-0817-4e9e-aa06-c5ce4928da5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute noise in MSE and NSE\n",
    "mse_bits = table.loc['MSE', 'bits']\n",
    "nse_bits = table.loc['NSE', 'bits']\n",
    "zmale_bits = table.loc['ZMALE', 'bits']\n",
    "mse_noise = round(100*(mse_bits - zmale_bits) / mse_bits)\n",
    "nse_noise = round(100*(nse_bits - zmale_bits) / nse_bits)\n",
    "print(f'MSE is at least {mse_noise} percent noise and NSE at least {nse_noise} percent.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c4d8534-2034-49a6-8083-8e765642b23f",
   "metadata": {},
   "source": [
    "## Data availability\n",
    "The streamflow data are from Russell et al. (2020) and are available at https://doi.org/10.5066/P9XT4WSP.\n",
    "This demonstration notebook is available at at https://code.usgs.gov/wma/hytest/workflow-hodson-2022-objective-benchmark. \n",
    "\n",
    "\n",
    "## References \n",
    "Burnham, K.P. and Anderson, D.R. (2002). Model selection and multimodel inference: A Practical Information-Theoretic Approach. 2nd Edition, Springer-Verlag, New York.\n",
    "\n",
    "Hodson, T.O. (2022). Root-mean-square error (RMSE) or mean absolute error (MAE): when to use them or not, Geosci. Model Dev., 15, 5481–5487. https://doi.org/10.5194/gmd-15-5481-2022\n",
    "\n",
    "Russell, A.M., Over, T.M., and Farmer, W.H. (2020). Cross-validation results for five statistical methods of daily streamflow estimation at 1,385 reference streamgages in the conterminous United States, Water Years 1981-2017: U.S. Geological Survey data release. https://doi.org/10.5066/P9XT4WSP"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}