03-bootstrap-nse-old.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 = np.sqrt(np.var(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": 3,
   "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://renc.osn.xsede.org\"})\n",
    "\n",
    "with fs_read.open('s3://rsignellbucket2/hytest/thodson/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": 5,
   "id": "858c72e9-c5d9-4a1a-ab24-6be7e6c37d9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read local copy\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "df = pd.read_parquet('gages2_nndar.parquet')\n",
    "df[df < 0.01] = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "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},\n",
    "    'MSE' : {'name':'mean squared error', 'f':normal_ll},\n",
    "    'NSE' : {'name':'normalized squared error', 'f':nse_ll},\n",
    "    'MAE' : {'name': 'mean absolute error', 'f':laplace_ll},\n",
    "    'MSPE' : {'name': 'mean squared percent error', 'f':mspe_ll},\n",
    "    'MSLE' : {'name':'mean squared log error*', 'f':lognormal_ll},\n",
    "    'MALE' : {'name':'mean absolute log error*', 'f':loglaplace_ll},\n",
    "    'ZMSLE' : {'name':'zero-inflated MSLE', 'f':zilognormal_ll},\n",
    "    'ZMALE' : {'name':'zero-inflated MALE', 'f':ziloglaplace_ll},\n",
    "    'MARE' : {'name':'mean absolute square root error', 'f':mare_ll},\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.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)\n",
    "\n",
    "# step 4: format output table\n",
    "\n",
    "table = obj_df[['name','bits','weight']].sort_values('weight').round(2)#.rename(columns=names)\n",
    "\n",
    "table['rank'] = len(table) - np.argsort(table['weight'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "b77462ab-d4c1-4f12-9dbe-2e5a3ba12326",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>name</th>\n",
       "      <th>f</th>\n",
       "      <th>bits</th>\n",
       "      <th>weight</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>U</th>\n",
       "      <td>uniformly distributed error</td>\n",
       "      <td>&lt;function uniform_ll at 0x7f35f7066480&gt;</td>\n",
       "      <td>18.172241</td>\n",
       "      <td>0.000093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSE</th>\n",
       "      <td>mean squared error</td>\n",
       "      <td>&lt;function normal_ll at 0x7f35f7066ac0&gt;</td>\n",
       "      <td>11.616919</td>\n",
       "      <td>0.008716</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSE</th>\n",
       "      <td>normalized squared error</td>\n",
       "      <td>&lt;function nse_ll at 0x7f36146080e0&gt;</td>\n",
       "      <td>11.195196</td>\n",
       "      <td>0.011675</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MAE</th>\n",
       "      <td>mean absolute error</td>\n",
       "      <td>&lt;function laplace_ll at 0x7f35f70668e0&gt;</td>\n",
       "      <td>9.492831</td>\n",
       "      <td>0.037994</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSPE</th>\n",
       "      <td>mean squared percent error</td>\n",
       "      <td>&lt;function mspe_ll at 0x7f35f70662a0&gt;</td>\n",
       "      <td>23.537709</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSLE</th>\n",
       "      <td>mean squared log error*</td>\n",
       "      <td>&lt;function lognormal_ll at 0x7f35f70660c0&gt;</td>\n",
       "      <td>7.466534</td>\n",
       "      <td>0.154772</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MALE</th>\n",
       "      <td>mean absolute log error*</td>\n",
       "      <td>&lt;function loglaplace_ll at 0x7f35f70663e0&gt;</td>\n",
       "      <td>7.039863</td>\n",
       "      <td>0.208034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZMSLE</th>\n",
       "      <td>zero-inflated MSLE</td>\n",
       "      <td>&lt;function zilognormal_ll at 0x7f363c2b4680&gt;</td>\n",
       "      <td>7.182802</td>\n",
       "      <td>0.188411</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZMALE</th>\n",
       "      <td>zero-inflated MALE</td>\n",
       "      <td>&lt;function ziloglaplace_ll at 0x7f35f7066660&gt;</td>\n",
       "      <td>6.952254</td>\n",
       "      <td>0.221059</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MARE</th>\n",
       "      <td>mean absolute square root error</td>\n",
       "      <td>&lt;function mare_ll at 0x7f3614608040&gt;</td>\n",
       "      <td>7.337572</td>\n",
       "      <td>0.169245</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                  name  \\\n",
       "U          uniformly distributed error   \n",
       "MSE                 mean squared error   \n",
       "NSE           normalized squared error   \n",
       "MAE                mean absolute error   \n",
       "MSPE        mean squared percent error   \n",
       "MSLE           mean squared log error*   \n",
       "MALE          mean absolute log error*   \n",
       "ZMSLE               zero-inflated MSLE   \n",
       "ZMALE               zero-inflated MALE   \n",
       "MARE   mean absolute square root error   \n",
       "\n",
       "                                                  f       bits    weight  \n",
       "U           <function uniform_ll at 0x7f35f7066480>  18.172241  0.000093  \n",
       "MSE          <function normal_ll at 0x7f35f7066ac0>  11.616919  0.008716  \n",
       "NSE             <function nse_ll at 0x7f36146080e0>  11.195196  0.011675  \n",
       "MAE         <function laplace_ll at 0x7f35f70668e0>   9.492831  0.037994  \n",
       "MSPE           <function mspe_ll at 0x7f35f70662a0>  23.537709  0.000002  \n",
       "MSLE      <function lognormal_ll at 0x7f35f70660c0>   7.466534  0.154772  \n",
       "MALE     <function loglaplace_ll at 0x7f35f70663e0>   7.039863  0.208034  \n",
       "ZMSLE   <function zilognormal_ll at 0x7f363c2b4680>   7.182802  0.188411  \n",
       "ZMALE  <function ziloglaplace_ll at 0x7f35f7066660>   6.952254  0.221059  \n",
       "MARE           <function mare_ll at 0x7f3614608040>   7.337572  0.169245  "
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obj_df"
   ]
  },
  {
   "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": 7,
   "id": "9a920adc-2381-4d72-bc8b-54024a412374",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>name</th>\n",
       "      <th>bits</th>\n",
       "      <th>weight</th>\n",
       "      <th>rank</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>MSPE</th>\n",
       "      <td>mean squared percent error</td>\n",
       "      <td>23.54</td>\n",
       "      <td>0.00</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>U</th>\n",
       "      <td>uniformly distributed error</td>\n",
       "      <td>18.17</td>\n",
       "      <td>0.00</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSE</th>\n",
       "      <td>mean squared error</td>\n",
       "      <td>11.62</td>\n",
       "      <td>0.01</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NSE</th>\n",
       "      <td>normalized squared error</td>\n",
       "      <td>11.20</td>\n",
       "      <td>0.01</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MAE</th>\n",
       "      <td>mean absolute error</td>\n",
       "      <td>9.49</td>\n",
       "      <td>0.04</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MSLE</th>\n",
       "      <td>mean squared log error*</td>\n",
       "      <td>7.47</td>\n",
       "      <td>0.15</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MARE</th>\n",
       "      <td>mean absolute square root error</td>\n",
       "      <td>7.34</td>\n",
       "      <td>0.17</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZMSLE</th>\n",
       "      <td>zero-inflated MSLE</td>\n",
       "      <td>7.18</td>\n",
       "      <td>0.19</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>MALE</th>\n",
       "      <td>mean absolute log error*</td>\n",
       "      <td>7.04</td>\n",
       "      <td>0.21</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZMALE</th>\n",
       "      <td>zero-inflated MALE</td>\n",
       "      <td>6.95</td>\n",
       "      <td>0.22</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                  name   bits  weight  rank\n",
       "MSPE        mean squared percent error  23.54    0.00    10\n",
       "U          uniformly distributed error  18.17    0.00     9\n",
       "MSE                 mean squared error  11.62    0.01     8\n",
       "NSE           normalized squared error  11.20    0.01     7\n",
       "MAE                mean absolute error   9.49    0.04     6\n",
       "MSLE           mean squared log error*   7.47    0.15     5\n",
       "MARE   mean absolute square root error   7.34    0.17     4\n",
       "ZMSLE               zero-inflated MSLE   7.18    0.19     3\n",
       "MALE          mean absolute log error*   7.04    0.21     2\n",
       "ZMALE               zero-inflated MALE   6.95    0.22     1"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#print(table.to_latex())\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38243edf-a33e-44c3-a9be-d3d6fbff2a2c",
   "metadata": {},
   "source": [
    "### Bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "98c22a51-d390-402b-8354-7980cb25fdb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample=100\n",
    "n = 1000\n",
    "b = 100\n",
    "seed=12345\n",
    "block_size = n//b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4670d52f-d7d4-4304-ab4c-6521ea949f7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nse_ll(y, y_hat, sigma_o=None, 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",
    "    if sigma_o is None:\n",
    "        sigma_o = y.groupby('gage_id').transform(lambda x: x.std(ddof=0))\n",
    "        \n",
    "    return normal_ll(y, y_hat, transform=lambda x: x/sigma_o, gradient=1/sigma_o)\n",
    "\n",
    "sigma_o = df['obs'].groupby('gage_id').transform(lambda x: x.std())\n",
    "sigma_global = sigma_o.groupby('gage_id').sample(sample, replace=True, random_state=seed)\n",
    "\n",
    "def nse2_ll(y, y_hat):\n",
    "    return nse_ll(y, y_hat, sigma_o=sigma_global)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "7878b76e-b51a-4553-878c-6701c1d9b1b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#modify nse for the bootstrap\n",
    "objectives = {\n",
    "    'U' : {'name':'uniformly distributed error', 'f':uniform_ll},\n",
    "    'MSE' : {'name':'mean squared error', 'f':normal_ll},\n",
    "    'NSE' : {'name':'normalized squared error', 'f':nse2_ll},\n",
    "    'MAE' : {'name': 'mean absolute error', 'f':laplace_ll},\n",
    "    'MSPE' : {'name': 'mean squared percent error', 'f':mspe_ll},\n",
    "    'MSLE' : {'name':'mean squared log error*', 'f':lognormal_ll},\n",
    "    'MALE' : {'name':'mean absolute log error*', 'f':loglaplace_ll},\n",
    "    'ZMSLE' : {'name':'zero-inflated MSLE', 'f':zilognormal_ll},\n",
    "    'ZMALE' : {'name':'zero-inflated MALE', 'f':ziloglaplace_ll},\n",
    "    'MARE' : {'name':'mean absolute square root error', 'f':mare_ll},\n",
    "}\n",
    "\n",
    "obj_df = pd.DataFrame.from_dict(objectives, orient='index')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea68fe4d-8124-4af2-93e9-f3176fd1079a",
   "metadata": {},
   "source": [
    "## Run Bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "id": "9873ca17-81e0-4151-9449-c89687d07828",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "boot_array = np.zeros((block_size, obj_df.f.shape[0]))\n",
    "\n",
    "for b in range(0,b): # for each block\n",
    "    for i in range(0,block_size):\n",
    "        state = b*block_size + i\n",
    "        #temp_df = df.sample(df.shape[0], replace=True, random_state=state)\n",
    "        temp_df = df.groupby('gage_id').sample(sample, replace=True, random_state=state)\n",
    "        \n",
    "        for j,f in enumerate(obj_df.f):\n",
    "            #print(f)\n",
    "            boot_array[i,j]= f(temp_df.obs, temp_df.NNDAR)/len(temp_df)/np.log(2)\n",
    "    \n",
    "    boot_df = pd.DataFrame(boot_array)\n",
    "    boot_df.columns = obj_df.index\n",
    "    boot_df.to_csv('bootstrap/h_boot_{}.csv'.format(b))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1551695d-ce7c-48f9-a34b-fbff302f9175",
   "metadata": {},
   "source": [
    "## read in bootstraps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "3b1b79e5-e146-4436-93b7-aba7172fffae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import glob\n",
    "import pandas as pd\n",
    "\n",
    "# Get CSV files list from a folder\n",
    "path = 'bootstrap'\n",
    "csv_files = glob.glob(path + \"/*.csv\")\n",
    "\n",
    "# Read each CSV file into DataFrame\n",
    "# This creates a list of dataframes\n",
    "df_list = (pd.read_csv(file) for file in csv_files)\n",
    "\n",
    "# Concatenate all DataFrames\n",
    "big_df   = pd.concat(df_list, ignore_index=True).iloc[:,1:]\n",
    "\n",
    "#big_df.columns = obj_df.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "f8841c94-9aa2-41d2-a538-74605c7fa5f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = table.index\n",
    "vars = big_df.var()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "c7b78bfd-054a-46e4-a34c-6aeeba23f13d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MSPE     0.553174\n",
       "U        0.374864\n",
       "MSE      0.022104\n",
       "NSE      2.827437\n",
       "MAE      0.000423\n",
       "MSLE     0.000049\n",
       "MARE     0.000050\n",
       "ZMSLE    0.000043\n",
       "MALE     0.000041\n",
       "ZMALE    0.000035\n",
       "dtype: float64"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vars[cols]"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9f45e02d-17eb-44bf-95a1-eda985b3b2eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>obs</th>\n",
       "      <th>NNDAR</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gage_id</th>\n",
       "      <th>date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">01013500</th>\n",
       "      <th>1980-10-01</th>\n",
       "      <td>509.0</td>\n",
       "      <td>649.55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-10-02</th>\n",
       "      <td>518.0</td>\n",
       "      <td>618.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-10-03</th>\n",
       "      <td>516.0</td>\n",
       "      <td>618.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-10-04</th>\n",
       "      <td>620.0</td>\n",
       "      <td>796.61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-10-05</th>\n",
       "      <td>759.0</td>\n",
       "      <td>1397.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">402114105350101</th>\n",
       "      <th>2017-09-26</th>\n",
       "      <td>23.4</td>\n",
       "      <td>23.32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-09-27</th>\n",
       "      <td>22.1</td>\n",
       "      <td>22.79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-09-28</th>\n",
       "      <td>23.8</td>\n",
       "      <td>26.72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-09-29</th>\n",
       "      <td>28.4</td>\n",
       "      <td>26.99</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-09-30</th>\n",
       "      <td>25.5</td>\n",
       "      <td>25.15</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>14067063 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                              obs    NNDAR\n",
       "gage_id         date                      \n",
       "01013500        1980-10-01  509.0   649.55\n",
       "                1980-10-02  518.0   618.91\n",
       "                1980-10-03  516.0   618.91\n",
       "                1980-10-04  620.0   796.61\n",
       "                1980-10-05  759.0  1397.14\n",
       "...                           ...      ...\n",
       "402114105350101 2017-09-26   23.4    23.32\n",
       "                2017-09-27   22.1    22.79\n",
       "                2017-09-28   23.8    26.72\n",
       "                2017-09-29   28.4    26.99\n",
       "                2017-09-30   25.5    25.15\n",
       "\n",
       "[14067063 rows x 2 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hyriver",
   "language": "python",
   "name": "hyriver"
  },
  "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.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}