Repository: astro-informatics/harmonic (language: python)

Issue:
Log-space variance (and update MAP website)
This just adds a very brief note to the documentation explaining how to consider the log(var(x)) in terms more directly comparable with var(log(x)).

I also updated my website link in the README.rst 


Discussion / hints:

Fix patch:
diff --git a/README.rst b/README.rst
index a3c2e7c6..841f7225 100644
--- a/README.rst
+++ b/README.rst
@@ -72,7 +72,15 @@ To check the install has worked correctly run the unit tests with
 .. code-block:: bash
 
     pytest 
-    
+
+To build the documentation from source run
+
+.. code-block:: bash
+
+    cd docs && make html
+
+Then open ``./docs/_build/html/index.html`` in a browser.
+
 Documentation
 =============
 
diff --git a/docs/background/Bayes_Factors/bayes_factors.rst b/docs/background/Bayes_Factors/bayes_factors.rst
index d06d78af..b20c2067 100644
--- a/docs/background/Bayes_Factors/bayes_factors.rst
+++ b/docs/background/Bayes_Factors/bayes_factors.rst
@@ -2,35 +2,40 @@
 Bayes factors
 **************
 
-We have so far considered the estimation of the reciprocal marginal likelihood and related variances only.  However it is the marginal likelihood itself (not its reciprocal), or the Bayes factors computed to compare two models, that is typically of direct interest.  We therefore consider how to compute these quantities of interest and a measure of their variance.
+We have so far considered the estimation of the reciprocal marginal likelihood and related variances only.  However, it is the marginal likelihood itself (not its reciprocal), or the Bayes factors computed to compare two models, that is typically of direct interest.  We therefore consider how to compute these quantities of interest and a measure of their variance (as discussed in `McEwen et al. 2021 <https://arxiv.org/abs/2111.12720>`_).
 
 First, consider the mean and variance of the function :math:`f(X,Y) = X/Y` of two uncorrelated random variables :math:`X` and :math:`Y`, which by Taylor expansion to second order are given by
 
 .. math::
 
   \mathbb{E}\biggl(\frac{X}{Y}\biggr) \simeq \frac{\mathbb{E}(X)}{\mathbb{E}(Y)} + \frac{\mathbb{E}(X)}{\mathbb{E}(Y)^3} \sigma_Y^2
-  \qquad \quad
-  \text{and}
-  \qquad \quad
+  
+and
+
+.. math::
+
   \text{var}\biggl(\frac{X}{Y}\biggr) \simeq \frac{1}{\mathbb{E}(Y)^2} \sigma_X^2 + \frac{\mathbb{E}(X)^2}{\mathbb{E}(Y)^4} \sigma_Y^2,
 
 respectively, where :math:`\sigma_X = \mathbb{E}\bigl[ (X - \mathbb{E}(X))^2 \bigr]` and :math:`\sigma_Y = \mathbb{E}\bigl[ (Y - \mathbb{E}(Y))^2 \bigr]`.
+
 Using this result the marginal likelihood and its variance can be estimated from the reciprocal estimators by making use of the relations
 
 .. math:: 
 
   \mathbb{E}( z )
   = \mathbb{E}\biggl(\frac{1}{{\rho}}\biggr)
-  \simeq \frac{1}{\mathbb{E}({\rho})} \biggl( 1 + \frac{\sigma_{\rho}^2}{\mathbb{E}({\rho})^2} \biggr)
-  = \hat{z}
-  \qquad \quad
-  \text{and}
+  \simeq \frac{1}{\mathbb{E}({\rho})} \biggl( 1 + \frac{\sigma_{\rho}^2}{\mathbb{E}({\rho})^2} \biggr)  
   \qquad \quad
+
+and 
+
+.. math:: 
+
+  \text{var}(z) =
   \text{var}\biggl(\frac{1}{{\rho}}\biggr)
-  \simeq \frac{\sigma_{\rho}^2}{\mathbb{E}({\rho})^4}
-  = \hat{\epsilon}^2,
+  \simeq \frac{\sigma_{\rho}^2}{\mathbb{E}({\rho})^4},
 
-respectively, by considering the case :math:`X=1` and :math:`Y = \rho`.
+respectively, by considering the case :math:`X=1` and :math:`Y = \rho`.  These expressions are computed (in log space) in ``harmonic`` by the :code:`Evidence.compute_ln_evidence` method.
 
 Typically it is the Bayes factor given by the ratio of marginal likelihoods that is of most interest in order to compare models.  Again using the expressions above for the mean and variance of the function :math:`f(X,Y) = X/Y`, this time for the case :math:`X = \rho_2` and :math:`Y = \rho_1`, the Bayes factor and its variance can be estimated directly from the reciprocal marginal likelihood estimates and variances by making use of the relations
 
@@ -53,4 +58,6 @@ and
   \simeq
   \frac{\mathbb{E}({\rho_1})^2 \sigma_{\rho_2}^2+ \mathbb{E}({\rho_2})^2 \sigma_{\rho_1}^2}{\mathbb{E}({\rho_1})^4},
 
-respectively.
\ No newline at end of file
+respectively.  These expressions are computed (in log space) in ``harmonic`` by the :code:`evidence.compute_ln_bayes_factor` function.
+
+While we compute log quantities, as discussed above, the variance computed (i.e. error estimate) is related to the evidence itself (not the log evidence).  We consider how to compute an error estimate for the log evidence in the section :ref:`logvar`.
diff --git a/docs/background/Log_space_statistics/log_variance.rst b/docs/background/Log_space_statistics/log_variance.rst
new file mode 100644
index 00000000..3eb0eb9e
--- /dev/null
+++ b/docs/background/Log_space_statistics/log_variance.rst
@@ -0,0 +1,37 @@
+.. _logvar:
+
+******************
+Log-space variance
+******************
+
+``harmonic`` computes estimators of the reciprocal evidence and its variance, denoted by :math:`\hat{\rho}` and :math:`\hat{\sigma}^2` respectively. 
+
+We compute the natural logarithm of these quantities for numerical stability (using a number of numerical tricks, such as the `LogSumExp trick <https://en.wikipedia.org/wiki/LogSumExp>`_).  We therefore directly compute :math:`\log ( \hat{\rho} )` and :math:`\log ( \hat{\sigma} )`.  
+
+While quantities are computed in log space for numerical stability, the variance that is computed relates to the underlying reciprocal evidence (**not** the log reciprocal evidence), i.e. the terms computed can be considered as the estimate and error :math:`\hat{\rho} \pm \hat{\sigma}`.  This is of course not equivalent to a variance estimate of the log space estimate, which can be seen straightforwardly since :math:`\log ( \text{var}(x) ) \not = \text{var}(\log (x ))` as summation and the natural logarithm do not commute.
+
+In some settings one may be interested in an error estimate defined in log space, i.e. the log-space error :math:`\hat{\zeta}_\pm` defined by
+
+.. math::
+
+  \log ( \hat{\rho} \pm \hat{\sigma} ) = \log (\hat{\rho}) + \hat{\zeta}_\pm .
+
+The log-space error estimate can be computed by
+
+.. math::
+
+  \hat{\zeta}_\pm = \log (\hat{\rho} \pm \hat{\sigma} ) - \log (\hat{\rho}) .
+
+This may also be expressed as 
+
+.. math::
+
+  \hat{\zeta}_\pm = \log(1 \pm \hat{\sigma} / \hat{\rho}) ,
+
+where 
+
+.. math::
+
+  \hat{\sigma} / \hat{\rho} =  \exp \bigl( \log(\hat{\sigma}) - \log(\hat{\rho}) \bigr) .
+
+The advantage of this second approach is that it avoids explicitly computing :math:`\exp(\log(\hat{\rho})) \pm \exp(\log(\hat{\sigma}))`, which has the potental to over- or under-flow.  The second, more stable, approach is implemented in :code:`Evidence.compute_ln_inv_evidence_errors`.
diff --git a/docs/index.rst b/docs/index.rst
index 580bbb6a..82b9f207 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -102,6 +102,7 @@ Referencing
    background/Harmonic_Estimator/learnt_harmonic_mean
    background/Machine_Learning/index
    background/Bayes_Factors/bayes_factors
+   background/Log_space_statistics/log_variance
 
 
 .. toctree::
diff --git a/harmonic/evidence.pyx b/harmonic/evidence.pyx
index e8d2a163..efc1db3b 100644
--- a/harmonic/evidence.pyx
+++ b/harmonic/evidence.pyx
@@ -404,6 +404,45 @@ class Evidence:
         return (ln_evidence, ln_evidence_std)
 
 
+    def compute_ln_inv_evidence_errors(self):
+        r"""Compute lower and uppper errors on the log_e of the inverse evidence. 
+
+        Compute the log-space error :math:`\hat{\zeta}_\pm` defined by
+
+        .. math::
+
+            \log ( \hat{\rho} \pm \hat{\sigma} ) = \log (\hat{\rho}) + \hat{\zeta}_\pm .
+
+        Computed in a numerically stable way by
+
+        .. math::
+
+            \hat{\zeta}_\pm = \log(1 \pm \hat{\sigma} / \hat{\rho}) .
+
+        Returns:
+
+            (double, double): Tuple containing the following.
+
+                - ln_evidence_err_neg (double): Lower error for log_e of inverse evidence.
+
+                - ln_evidence_err_pos (double): Upper error for log_e of inverse evidence.
+
+        """
+
+        ln_ratio = 0.5*self.ln_evidence_inv_var - self.ln_evidence_inv
+
+        ratio = np.exp(ln_ratio)
+
+        if np.abs(ratio - 1.0) > 1e-8:
+            ln_evidence_err_neg = np.log( 1.0 - ratio )
+        else:
+            ln_evidence_err_neg = np.NINF
+
+        ln_evidence_err_pos = np.log( 1.0 + ratio )
+
+        return (ln_evidence_err_neg, ln_evidence_err_pos)
+
+
     def serialize(self, filename):
         """Serialize evidence object.
 
diff --git a/notebooks/checkpointing.ipynb b/notebooks/checkpointing.ipynb
index 5c804db6..dab60dba 100644
--- a/notebooks/checkpointing.ipynb
+++ b/notebooks/checkpointing.ipynb
@@ -1,29 +1,10 @@
 {
-  "nbformat": 4,
-  "nbformat_minor": 0,
-  "metadata": {
-    "colab": {
-      "name": "checkpointing.ipynb",
-      "provenance": [],
-      "collapsed_sections": [],
-      "toc_visible": true,
-      "authorship_tag": "ABX9TyPeQz+AuMy/ROmUgYTB/C7J",
-      "include_colab_link": true
-    },
-    "kernelspec": {
-      "name": "python3",
-      "display_name": "Python 3"
-    },
-    "language_info": {
-      "name": "python"
-    }
-  },
   "cells": [
     {
       "cell_type": "markdown",
       "metadata": {
-        "id": "view-in-github",
-        "colab_type": "text"
+        "colab_type": "text",
+        "id": "view-in-github"
       },
       "source": [
         "<a href=\"https://colab.research.google.com/github/astro-informatics/harmonic/blob/main/notebooks/checkpointing.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
@@ -58,31 +39,31 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "XSCBJ_S4XvDe"
       },
+      "outputs": [],
       "source": [
         "%%capture\n",
         "# Install packages\n",
         "%pip install harmonic emcee corner getdist"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "1OtVkBf7OAIF"
       },
+      "outputs": [],
       "source": [
         "%%capture\n",
         "# Import modules\n",
         "import numpy as np\n",
         "import harmonic as hm\n",
         "import emcee"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -98,9 +79,11 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "8PgO6f4VQSpD"
       },
+      "outputs": [],
       "source": [
         "def ln_posterior(x, inv_cov):\n",
         "    \"\"\"Compute log_e of n dimensional multivariate gaussian.\n",
@@ -132,9 +115,7 @@
         "    diag_cov = np.ones(ndim)\n",
         "    np.fill_diagonal(cov, diag_cov)\n",
         "    return cov"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -149,9 +130,11 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "xjsxOkX0QSzL"
       },
+      "outputs": [],
       "source": [
         "# Define parameters for emcee sampling\n",
         "ndim = 10                   # number of dimensions\n",
@@ -165,9 +148,7 @@
         "# Create covariance matrix \n",
         "cov = init_cov(ndim)\n",
         "inv_cov = np.linalg.inv(cov) "
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -180,9 +161,11 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "Ju_YJlk2QS7x"
       },
+      "outputs": [],
       "source": [
         "# Instantiate hyper-sphere model\n",
         "domains_sphere = [np.array([1E0,1E0])] # Not used since not training model.\n",
@@ -191,9 +174,7 @@
         "# Fit model by hand by fixing radius\n",
         "model.set_R(np.sqrt(ndim)) # A good rule of thumb for a Gaussian with a diagonal covariance\n",
         "model.fitted = True"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -209,9 +190,11 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "us1kBuWlQZTy"
       },
+      "outputs": [],
       "source": [
         "# Set initial random position and state\n",
         "pos = np.random.rand(ndim * nchains).reshape((nchains, ndim)) * 0.1   \n",
@@ -246,9 +229,7 @@
         "    \n",
         "    # Clear memory \n",
         "    del chains, samples, lnprob, sampler, prob"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -272,14 +253,14 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "Eel3bSORQZW0"
       },
+      "outputs": [],
       "source": [
         "ln_evidence, ln_evidence_std = cal_ev.compute_ln_evidence()"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -295,9 +276,11 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "eSxxNW1KQZZc"
       },
+      "outputs": [],
       "source": [
         "def ln_analytic_evidence(ndim, cov):\n",
         "    \"\"\"Compute analytic ln_e evidence.\n",
@@ -315,9 +298,7 @@
         "    \"\"\"\n",
         "    ln_norm_lik = -0.5*ndim*np.log(2*np.pi)-0.5*np.log(np.linalg.det(cov))   \n",
         "    return -ln_norm_lik"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -330,14 +311,14 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "efPNgW8qQZcW"
       },
+      "outputs": [],
       "source": [
         "ln_evidence_analytic = ln_analytic_evidence(ndim, cov)"
-      ],
-      "execution_count": null,
-      "outputs": []
+      ]
     },
     {
       "cell_type": "markdown",
@@ -350,40 +331,34 @@
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
         "id": "WrB47hA3QZfM"
       },
+      "outputs": [],
       "source": [
-        "ln_evidence_err = np.log(np.exp(ln_evidence) + np.exp(ln_evidence_std)) - ln_evidence"
-      ],
-      "execution_count": null,
-      "outputs": []
+        "ln_evidence_err_neg, ln_evidence_err_pos = cal_ev.compute_ln_inv_evidence_errors()\n",
+        "ln_evidence_err_avg = (ln_evidence_err_pos - ln_evidence_err_neg) / 2.0\n",
+        "print(f\"ln_evidence_err_neg = {ln_evidence_err_neg}\")\n",
+        "print(f\"ln_evidence_err_pos = {ln_evidence_err_pos}\")\n",
+        "print(f\"ln_evidence_err_avg = {ln_evidence_err_avg}\")"
+      ]
     },
     {
       "cell_type": "code",
+      "execution_count": null,
       "metadata": {
-        "id": "yP3AjUrCQZhh",
         "colab": {
           "base_uri": "https://localhost:8080/"
         },
+        "id": "yP3AjUrCQZhh",
         "outputId": "cc01f114-07e0-4abe-ea2f-72fede5943f8"
       },
+      "outputs": [],
       "source": [
-        "print('ln_evidence (harmonic) = {} +/- {}'.format(ln_evidence, ln_evidence_err))\n",
+        "print('ln_evidence (harmonic) = {} (+ {} / {})'.format(ln_evidence, ln_evidence_err_pos, ln_evidence_err_neg))\n",
         "print('ln_evidence (analytic) = {}'.format(ln_evidence_analytic))\n",
-        "print('nsigma = {}'.format(np.abs(ln_evidence - ln_evidence_analytic) / ln_evidence_err))"
-      ],
-      "execution_count": null,
-      "outputs": [
-        {
-          "output_type": "stream",
-          "name": "stdout",
-          "text": [
-            "ln_evidence (harmonic) = 9.191328725973255 +/- 0.0032251202241564414\n",
-            "ln_evidence (analytic) = 9.189385332046726\n",
-            "nsigma = 0.6025803044405069\n"
-          ]
-        }
+        "print('nsigma = {}'.format(np.abs(ln_evidence - ln_evidence_analytic) / ln_evidence_err_avg))"
       ]
     },
     {
@@ -395,5 +370,33 @@
         "As expected, the evidence computed by `harmonic` is close to that computed analytically."
       ]
     }
-  ]
-}
\ No newline at end of file
+  ],
+  "metadata": {
+    "colab": {
+      "authorship_tag": "ABX9TyPeQz+AuMy/ROmUgYTB/C7J",
+      "collapsed_sections": [],
+      "include_colab_link": true,
+      "name": "checkpointing.ipynb",
+      "provenance": [],
+      "toc_visible": true
+    },
+    "kernelspec": {
+      "display_name": "Python 3",
+      "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.9.16"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}
diff --git a/requirements/requirements-examples.txt b/requirements/requirements-examples.txt
index 4f3a32d9..fb2adbee 100644
--- a/requirements/requirements-examples.txt
+++ b/requirements/requirements-examples.txt
@@ -1,4 +1,4 @@
-emcee==3.0.2
+emcee==3.1.4
 matplotlib==3.2.1
 corner==2.0.1
-getdist==1.1.0
\ No newline at end of file
+getdist==1.4.3


Test patch:
diff --git a/tests/test_evidence.py b/tests/test_evidence.py
index 90f77b8f..41cc974a 100644
--- a/tests/test_evidence.py
+++ b/tests/test_evidence.py
@@ -235,6 +235,63 @@ def test_compute_evidence():
     assert evidence_std == pytest.approx(np.exp(ln_evidence_std))
 
 
+def test_compute_ln_inv_evidence_errors():
+    ndim = 2
+    nchains = 100
+
+    domain = [np.array([1e-1, 1e1])]
+    sphere = md.HyperSphere(ndim, domain)
+    sphere.fitted = True
+
+    # Check boundary case where ratio 1.0
+    # (ln_ev_inv_var = 2 * ln_ev_inv)
+    ln_ev_inv = 10
+    ln_ev_inv_var = 2 * ln_ev_inv
+    ev = cbe.Evidence(nchains, sphere)
+    ev.evidence_inv = np.exp(ln_ev_inv)
+    ev.evidence_inv_var = np.exp(ln_ev_inv_var)
+    ev.ln_evidence_inv = ln_ev_inv
+    ev.ln_evidence_inv_var = ln_ev_inv_var
+
+    zeta_neg, zeta_pos = ev.compute_ln_inv_evidence_errors()
+    assert zeta_neg == np.NINF
+    assert zeta_pos == pytest.approx(np.log(2.0))
+
+    # Check case where ln_ev_inv_var = ln_ev_inv
+    ln_ev_inv = 10
+    ln_ev_inv_var = ln_ev_inv
+    ev = cbe.Evidence(nchains, sphere)
+    ev.evidence_inv = np.exp(ln_ev_inv)
+    ev.evidence_inv_var = np.exp(ln_ev_inv_var)
+    ev.ln_evidence_inv = ln_ev_inv
+    ev.ln_evidence_inv_var = ln_ev_inv_var
+
+    zeta_neg, zeta_pos = ev.compute_ln_inv_evidence_errors()
+    assert zeta_neg == pytest.approx(
+        np.log(1.0 - np.exp(0.5 * ln_ev_inv_var - ln_ev_inv))
+    )
+    assert zeta_pos == pytest.approx(
+        np.log(1.0 + np.exp(0.5 * ln_ev_inv_var - ln_ev_inv))
+    )
+
+    # Check case where ln_ev_inv_var = 0.5 * ln_ev_inv
+    ln_ev_inv = 10
+    ln_ev_inv_var = 0.5 * ln_ev_inv
+    ev = cbe.Evidence(nchains, sphere)
+    ev.evidence_inv = np.exp(ln_ev_inv)
+    ev.evidence_inv_var = np.exp(ln_ev_inv_var)
+    ev.ln_evidence_inv = ln_ev_inv
+    ev.ln_evidence_inv_var = ln_ev_inv_var
+
+    zeta_neg, zeta_pos = ev.compute_ln_inv_evidence_errors()
+    assert zeta_neg == pytest.approx(
+        np.log(1.0 - np.exp(0.5 * ln_ev_inv_var - ln_ev_inv))
+    )
+    assert zeta_pos == pytest.approx(
+        np.log(1.0 + np.exp(0.5 * ln_ev_inv_var - ln_ev_inv))
+    )
+
+
 def test_compute_bayes_factors():
     ndim = 2
     nchains = 100