Merge pull request #134 from ihincks/upgrade-vectorized-risk

ihincks · web-flow · commit 8170c84a0be1 · 2018-08-17T16:03:00.000-04:00
Rewrote bayes_risk to allow multiple expparams
diff --git a/src/qinfer/derived_models.py b/src/qinfer/derived_models.py
@@ -46,6 +46,7 @@
     'DerivedModel',
     'PoisonedModel',
     'BinomialModel',
+    'DifferentiableBinomialModel',
     'GaussianHyperparameterizedModel',
     'MultinomialModel',
     'MLEModel',
diff --git a/src/qinfer/smc.py b/src/qinfer/smc.py
@@ -552,88 +552,115 @@ def resample(self):
 
     def bayes_risk(self, expparams):
         r"""
-        Calculates the Bayes risk for a hypothetical experiment, assuming the
+        Calculates the Bayes risk for hypothetical experiments, assuming the
         quadratic loss function defined by the current model's scale matrix
         (see :attr:`qinfer.abstract_model.Simulatable.Q`).
 
-        :param expparams: The experiment at which to compute the Bayes risk.
+        :param expparams: The experiments at which to compute the risk.
         :type expparams: :class:`~numpy.ndarray` of dtype given by the current
             model's :attr:`~qinfer.abstract_model.Simulatable.expparams_dtype` property,
             and of shape ``(1,)``
 
-        :return float: The Bayes risk for the current posterior distribution
-            of the hypothetical experiment ``expparams``.
+        :return np.ndarray: The Bayes risk for the current posterior distribution
+            at each hypothetical experiment in ``expparams``, therefore 
+            has shape ``(expparams.size,)``
         """
-        # This subroutine computes the bayes risk for a hypothetical experiment
-        # defined by expparams.
-
-        # Assume expparams is a single experiment
-
-        # expparams =
-        # Q = np array(Nmodelparams), which contains the diagonal part of the
-        #     rescaling matrix.  Non-diagonal could also be considered, but
-        #     for the moment this is not implemented.
-        nout = self.model.n_outcomes(expparams) # This is a vector so this won't work
-        w, N = self.hypothetical_update(np.arange(nout), expparams, return_normalization=True)
-        w = w[:, 0, :] # Fix w.shape == (n_outcomes, n_particles).
-        N = N[:, :, 0] # Fix L.shape == (n_outcomes, n_particles).
-
-        xs = self.particle_locations.transpose([1, 0]) # shape (n_mp, n_particles).
-
-        # In the following, we will use the subscript convention that
-        # "o" refers to an outcome, "p" to a particle, and
-        # "i" to a model parameter.
-        # Thus, mu[o,i] is the sum over all particles of w[o,p] * x[i,p].
-
-        mu = np.transpose(np.tensordot(w,xs,axes=(1,1)))
-        var = (
-            # This sum is a reduction over the particle index and thus
-            # represents an expectation value over the diagonal of the
-            # outer product $x . x^T$.
-
-            np.transpose(np.tensordot(w,xs**2,axes=(1,1)))
-            # We finish by subracting from the above expectation value
-            # the diagonal of the outer product $mu . mu^T$.
-            - mu**2).T
 
+        # for models whose outcome number changes with experiment, we 
+        # take the easy way out and for-loop over experiments
+        n_eps = expparams.size
+        if n_eps > 1 and not self.model.is_n_outcomes_constant:
+            risk = np.empty(n_eps)
+            for idx in range(n_eps):
+                risk[idx] = self.bayes_risk(expparams[idx, np.newaxis])
+            return risk
+        
+        # outcomes for the first experiment
+        os = self.model.domain(expparams[0,np.newaxis])[0].values
+
+        # compute the hypothetical weights, likelihoods and normalizations for
+        # every possible outcome and expparam
+        # the likelihood over outcomes should sum to 1, so don't compute for last outcome
+        w_hyp, L, N = self.hypothetical_update(
+                os[:-1], 
+                expparams, 
+                return_normalization=True, 
+                return_likelihood=True
+            )
+        w_hyp_last_outcome = (1 - L.sum(axis=0)) * self.particle_weights[np.newaxis, :]
+        N = np.concatenate([N[:,:,0], np.sum(w_hyp_last_outcome[np.newaxis,:,:], axis=2)], axis=0)
+        w_hyp_last_outcome = w_hyp_last_outcome / N[-1,:,np.newaxis]
+        w_hyp = np.concatenate([w_hyp, w_hyp_last_outcome[np.newaxis,:,:]], axis=0)
+        # w_hyp.shape == (n_out, n_eps, n_particles)
+        # N.shape == (n_out, n_eps)
+
+        # compute the hypothetical means and variances given outcomes and exparams
+        # mu_hyp.shape == (n_out, n_eps, n_models)
+        # var_hyp.shape == (n_out, n_eps)
+        mu_hyp = np.dot(w_hyp, self.particle_locations)
+        var_hyp = np.sum(
+            w_hyp * 
+            np.sum(self.model.Q * (
+                self.particle_locations[np.newaxis,np.newaxis,:,:] - 
+                mu_hyp[:,:,np.newaxis,:]
+            ) ** 2,  axis=3),
+            axis=2
+        )
 
-        rescale_var = np.sum(self.model.Q * var, axis=1)
-        # Q has shape (n_mp,), therefore rescale_var has shape (n_outcomes,).
-        tot_norm = np.sum(N, axis=1)
-        return np.dot(tot_norm.T, rescale_var)
+        # the risk of a given expparam can be calculated as the mean posterior
+        # variance weighted over all possible outcomes
+        return np.sum(N * var_hyp, axis=0)
 
     def expected_information_gain(self, expparams):
         r"""
-        Calculates the expected information gain for a hypothetical experiment.
+        Calculates the expected information gain for each hypothetical experiment.
 
-        :param expparams: The experiment at which to compute expected
+        :param expparams: The experiments at which to compute expected
             information gain.
         :type expparams: :class:`~numpy.ndarray` of dtype given by the current
             model's :attr:`~qinfer.abstract_model.Simulatable.expparams_dtype` property,
-            and of shape ``(1,)``
+            and of shape ``(n,)``
 
-        :return float: The Bayes risk for the current posterior distribution
-            of the hypothetical experiment ``expparams``.
+        :return float: The expected information gain for each 
+            hypothetical experiment in ``expparams``.
         """
-
-        nout = self.model.n_outcomes(expparams)
-        w, N = self.hypothetical_update(np.arange(nout), expparams, return_normalization=True)
-        w = w[:, 0, :] # Fix w.shape == (n_outcomes, n_particles).
-        N = N[:, :, 0] # Fix N.shape == (n_outcomes, n_particles).
-
         # This is a special case of the KL divergence estimator (see below),
         # in which the other distribution is guaranteed to share support.
-        #
-        # KLD[idx_outcome] = Sum over particles(self * log(self / other[idx_outcome])
-        # Est. KLD = E[KLD[idx_outcome] | outcomes].
+        
+        # for models whose outcome number changes with experiment, we 
+        # take the easy way out and for-loop over experiments
+        n_eps = expparams.size
+        if n_eps > 1 and not self.model.is_n_outcomes_constant:
+            risk = np.empty(n_eps)
+            for idx in range(n_eps):
+                risk[idx] = self.expected_information_gain(expparams[idx, np.newaxis])
+            return risk
+        
+        # number of outcomes for the first experiment
+        os = self.model.domain(expparams[0,np.newaxis])[0].values
+
+        # compute the hypothetical weights, likelihoods and normalizations for
+        # every possible outcome and expparam
+        # the likelihood over outcomes should sum to 1, so don't compute for last outcome
+        w_hyp, L, N = self.hypothetical_update(
+                os[:-1], 
+                expparams, 
+                return_normalization=True, 
+                return_likelihood=True
+            )
+        w_hyp_last_outcome = (1 - L.sum(axis=0)) * self.particle_weights[np.newaxis, :]
+        N = np.concatenate([N[:,:,0], np.sum(w_hyp_last_outcome[np.newaxis,:,:], axis=2)], axis=0)
+        w_hyp_last_outcome = w_hyp_last_outcome / N[-1,:,np.newaxis]
+        w_hyp = np.concatenate([w_hyp, w_hyp_last_outcome[np.newaxis,:,:]], axis=0)
+        # w_hyp.shape == (n_out, n_eps, n_particles)
+        # N.shape == (n_out, n_eps)
 
-        KLD = np.sum(
-            w * np.log(w / self.particle_weights ),
-            axis=1 # Sum over particles.
-        )
+        # compute the Kullback-Liebler divergence for every experiment and possible outcome
+        # KLD.shape == (n_out, n_eps)
+        KLD = np.sum(w_hyp * np.log(w_hyp / self.particle_weights), axis=2)
 
-        tot_norm = np.sum(N, axis=1)
-        return np.dot(tot_norm, KLD)
+        # return the expected KLD (ie expected info gain) for every experiment
+        return np.sum(N * KLD, axis=0)
 
     ## MISC METHODS ###########################################################
 
@@ -1212,4 +1239,3 @@ def update(self, outcome, expparams,check_for_resample=True):
 
         # We now can update as normal.
         SMCUpdater.update(self, outcome, expparams,check_for_resample=check_for_resample)
-
diff --git a/src/qinfer/tests/test_metrics.py b/src/qinfer/tests/test_metrics.py
@@ -0,0 +1,165 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+##
+# test_metrics.py: Tests various metrics like risk and information gain.
+##
+# © 2014 Chris Ferrie (csferrie@gmail.com) and
+#        Christopher E. Granade (cgranade@gmail.com)
+#     
+# This file is a part of the Qinfer project.
+# Licensed under the AGPL version 3.
+##
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Affero General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+##
+
+## FEATURES ###################################################################
+
+from __future__ import division # Ensures that a/b is always a float.
+from __future__ import absolute_import
+## IMPORTS ####################################################################
+
+import numpy as np
+from numpy.testing import assert_equal, assert_almost_equal, assert_array_less,assert_approx_equal
+
+from qinfer.tests.base_test import DerandomizedTestCase
+from qinfer import (BinomialModel,CoinModel,BetaDistribution,DifferentiableBinomialModel)
+
+from qinfer.smc import SMCUpdater,SMCUpdaterBCRB
+
+class TestBayesRisk(DerandomizedTestCase):
+    # Test the implementation of numerical Bayes Risk by comparing to 
+    # numbers which were derived by doing analytic/numeric
+    # integrals of simple models in Mathematica. This test trusts that 
+    # these calculations were done correctly.
+
+    ALPHA = 1.
+    BETA = 3.
+    PRIOR_BETA = BetaDistribution(alpha=ALPHA, beta=BETA)
+    N_PARTICLES = 10000
+    NMEAS_EXPPARAMS = np.arange(1, 11, dtype=int)
+    
+    def setUp(self):
+
+        super(TestBayesRisk,self).setUp()
+        
+        # Set up relevant models.
+        self.coin_model = CoinModel()
+        self.binomial_model = BinomialModel(self.coin_model)
+
+        # Set up updaters for these models using particle approximations 
+        # of conjugate priors
+        self.updater_binomial = SMCUpdater(self.binomial_model,
+                TestBayesRisk.N_PARTICLES,TestBayesRisk.PRIOR_BETA)
+
+    def test_finite_outcomes_risk(self):
+        # The binomial model has a finite number of outcomes. Test the 
+        # risk calculation in this case.
+
+        expparams = self.NMEAS_EXPPARAMS.astype(self.binomial_model.expparams_dtype)
+
+        # estimate the risk
+        est_risk = self.updater_binomial.bayes_risk(expparams)
+
+        # compute exact risk
+        a, b = TestBayesRisk.ALPHA, TestBayesRisk.BETA
+        exact_risk = a * b / ((a + b) * (a + b + 1) * (a + b + expparams['n_meas']))
+
+        # see if they roughly match
+        assert_almost_equal(est_risk, exact_risk, decimal=3)
+
+class TestInformationGain(DerandomizedTestCase):
+    # Test the implementation of numerical information gain by comparing to 
+    # numbers which were derived by doing analytic/numeric
+    # integrals of simple models (binomialm, poisson, and gaussian) in 
+    # Mathematica. This test trusts that these calculations
+    # were done correctly.
+
+    ALPHA = 1
+    BETA = 3
+    PRIOR_BETA = BetaDistribution(alpha=ALPHA, beta=BETA)
+    N_PARTICLES = 10000
+    # Calculated in Mathematica, IG for the binomial model and the given expparams
+    NMEAS_EXPPARAMS = np.arange(1, 11, dtype=int)
+    BINOM_IG = np.array([0.104002,0.189223,0.261496,0.324283,0.379815,0.429613,0.474764,0.516069,0.554138,0.589446])
+    
+    def setUp(self):
+
+        super(TestInformationGain,self).setUp()
+        
+        # Set up relevant models.
+        self.coin_model = CoinModel()
+        self.binomial_model = BinomialModel(self.coin_model)
+        
+        # Set up updaters for these models using particle approximations 
+        # of conjugate priors
+        self.updater_binomial = SMCUpdater(self.binomial_model,
+                TestInformationGain.N_PARTICLES,TestInformationGain.PRIOR_BETA)
+
+
+    def test_finite_outcomes_ig(self):
+        # The binomial model has a finite number of outcomes. Test the 
+        # ig calculation in this case.
+
+        expparams = self.NMEAS_EXPPARAMS.astype(self.binomial_model.expparams_dtype)
+
+        # estimate the information gain
+        est_ig = self.updater_binomial.expected_information_gain(expparams)
+
+        # see if they roughly match
+        assert_almost_equal(est_ig, TestInformationGain.BINOM_IG, decimal=2)
+
+class TestFisherInformation(DerandomizedTestCase):
+    # Test the implementation of numerical Fisher Information by comparing to 
+    # numbers which were derived by doing analytic/numeric
+    # integrals of simple models (binomialm, poisson, and gaussian) in 
+    # Mathematica. This test trusts that these calculations
+    # were done correctly.
+
+    ALPHA = 1
+    BETA = 3
+    PRIOR_BETA = BetaDistribution(alpha=ALPHA, beta=BETA)
+    N_PARTICLES = 10000
+
+    BIN_FI_MODELPARAMS = np.linspace(0.01,0.99,5)
+    NMEAS_EXPPARAMS = np.arange(1, 11, dtype=int)
+    
+    def setUp(self):
+
+        super(TestFisherInformation,self).setUp()
+        
+        # Set up relevant models.
+        self.coin_model = CoinModel()
+        self.binomial_model = DifferentiableBinomialModel(self.coin_model)
+
+        # Set up updaters for these models using particle approximations 
+        # of conjugate priors
+        self.updater_binomial = SMCUpdater(self.binomial_model,
+                TestFisherInformation.N_PARTICLES,TestFisherInformation.PRIOR_BETA)
+
+
+    def test_finite_outcomes_fi(self):
+        # The binomial model has a finite number of outcomes. Test the 
+        # ig calculation in this case.
+
+        expparams = self.NMEAS_EXPPARAMS.astype(self.binomial_model.expparams_dtype)
+        p = TestFisherInformation.BIN_FI_MODELPARAMS
+        # estimate the information gain
+        est_fi = self.binomial_model.fisher_information(TestFisherInformation.BIN_FI_MODELPARAMS,expparams)[0,0]
+        p = p[:,np.newaxis]
+        n = expparams.astype(np.float32)[np.newaxis,:]
+
+        exact_fi = n/((1-p)*p) 
+        # see if they roughly match
+        
+        assert_almost_equal(est_fi,exact_fi, decimal=3)
