Merge pull request #146 from QInfer/cgranade/unknown-T2

Add unknown T₂, Gaussian hyperparameterized models.

Author: cgranade

doc/source/apiref/derived_models.rst

@@ -4,9 +4,9 @@
     license, visit http://creativecommons.org/licenses/by-nc-sa/3.0/ or send a
     letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
     California, 94041, USA.
-    
+
 .. _derived_models:
-    
+
 .. currentmodule:: qinfer
 
 Derived Models
@@ -30,6 +30,12 @@ additional functionality or changing the behaviors of underlying models.
 .. autoclass:: BinomialModel
     :members:
 
+:class:`GaussianHyperparameterizedModel` - Model over Gaussian outcomes conditioned on two-outcome experiments
+--------------------------------------------------------------------------------------------------------------
+
+.. autoclass:: GaussianHyperparameterizedModel
+    :members:
+
 :class:`MultinomialModel` - Model over batches of D-outcome experiments
 ----------------------------------------------------------------------
 
@@ -41,13 +47,13 @@ additional functionality or changing the behaviors of underlying models.
 
 .. autoclass:: MLEModel
     :members:
-    
+
 :class:`RandomWalkModel` - Model for adding fixed random walk to parameters
 ---------------------------------------------------------------------------
 
 .. autoclass:: RandomWalkModel
     :members:
-    
+
 :class:`GaussianRandomWalkModel` - Model for adding gaussian random walk to parameters
 --------------------------------------------------------------------------------------
 

doc/source/apiref/test_models.rst

@@ -28,22 +28,28 @@ built on top of QInfer.
 .. autoclass:: SimplePrecessionModel
     :members:
 
+:class:`UnknownT2Model` - Model of a single qubit undergoing Larmor precession with finite decoherence
+------------------------------------------------------------------------------------------------------
+
+.. autoclass:: UnknownT2Model
+    :members:
+
 :class:`NoisyCoinModel` - Classical coin flip model corrupted by a noisy process
 --------------------------------------------------------------------------------
 
 .. autoclass:: NoisyCoinModel
     :members:
-    
-:class:`NDieModel` 
+
+:class:`NDieModel`
 ------------------
 
 .. autoclass:: NDieModel
     :members:
-    
+
 Custom Models
 -------------
 
-Writing custom models is standard practice for QInfer users. 
+Writing custom models is standard practice for QInfer users.
 See :ref:`CustomModels`.
 
 .. currentmodule:: qinfer.tests.base_test

src/qinfer/derived_models.py

@@ -46,6 +46,7 @@
     'DerivedModel',
     'PoisonedModel',
     'BinomialModel',
+    'GaussianHyperparameterizedModel',
     'MultinomialModel',
     'MLEModel',
     'RandomWalkModel',
@@ -59,14 +60,14 @@
 from past.builtins import basestring
 
 import numpy as np
-from scipy.stats import binom, multivariate_normal
+from scipy.stats import binom, multivariate_normal, norm
 from itertools import combinations_with_replacement as tri_comb
 
 from qinfer.utils import binomial_pdf, multinomial_pdf, sample_multinomial
 from qinfer.abstract_model import Model, DifferentiableModel
 from qinfer._lib import enum # <- TODO: replace with flufl.enum!
 from qinfer.utils import binom_est_error
-from qinfer.domains import IntegerDomain, MultinomialDomain
+from qinfer.domains import IntegerDomain, MultinomialDomain, RealDomain
 
 ## FUNCTIONS ###################################################################
 
@@ -245,9 +246,9 @@ def __init__(self, underlying_model):
         else:
             self._expparams_scalar = False
             self._expparams_dtype = underlying_model.expparams_dtype + [('n_meas', 'uint')]
-    
+
     ## PROPERTIES ##
-        
+
     @property
     def decorated_model(self):
         # Provided for backcompat only.
@@ -380,6 +381,131 @@ def fisher_information(self, modelparams, expparams):
         )
         return two_outcome_fi * expparams['n_meas']
 
+class GaussianHyperparameterizedModel(DerivedModel):
+    """
+    Model representing a two-outcome model viewed through samples
+    from one of two distinct Gaussian distributions. This model adds four new
+    model parameters to its underlying model, respectively representing the
+    mean outcome conditioned on an underlying 0, the mean outcome conditioned
+    on an underlying 1, and the variance of outcomes conditioned in each case.
+
+    :param qinfer.abstract_model.Model underlying_model: An instance of a two-
+        outcome model to be viewed through Gaussian distributions.
+    """
+
+    def __init__(self, underlying_model):
+        super(GaussianHyperparameterizedModel, self).__init__(underlying_model)
+
+        if not (underlying_model.is_n_outcomes_constant and underlying_model.n_outcomes(None) == 2):
+            raise ValueError("Decorated model must be a two-outcome model.")
+
+        n_orig_mps = underlying_model.n_modelparams
+        self._orig_mps_slice = np.s_[:n_orig_mps]
+        self._mu_slice = np.s_[n_orig_mps:n_orig_mps + 2]
+        self._sigma2_slice = np.s_[n_orig_mps + 2:n_orig_mps + 4]
+
+    ## PROPERTIES ##
+
+    @property
+    def decorated_model(self):
+        # Provided for backcompat only.
+        return self.underlying_model
+
+    @property
+    def modelparam_names(self):
+        return self.underlying_model.modelparam_names + [
+            r'\mu_0', r'\mu_1',
+            r'\sigma_0^2', r'\sigma_1^2'
+        ]
+
+    @property
+    def n_modelparams(self):
+        return len(self.modelparam_names)
+
+    ## METHODS ##
+
+    def domain(self, expparams):
+        return [RealDomain()] * len(expparams) if expparams is not None else RealDomain()
+
+    def are_expparam_dtypes_consistent(self, expparams):
+        return True
+
+    def are_models_valid(self, modelparams):
+        orig_mps = modelparams[:, self._orig_mps_slice]
+        sigma2 = modelparams[:, self._sigma2_slice]
+
+        return np.all([
+            self.underlying_model.are_models_valid(orig_mps),
+            np.all(sigma2 > 0, axis=-1)
+        ], axis=0)
+
+    def underlying_likelihood(self, binary_outcomes, modelparams, expparams):
+        """
+        Given outcomes hypothesized for the underlying model, returns the likelihood
+        which which those outcomes occur.
+        """
+        original_mps = modelparams[..., self._orig_mps_slice]
+        return self.underlying_model.likelihood(binary_outcomes, original_mps, expparams)
+
+    def likelihood(self, outcomes, modelparams, expparams):
+        # By calling the superclass implementation, we can consolidate
+        # call counting there.
+        super(GaussianHyperparameterizedModel, self).likelihood(outcomes, modelparams, expparams)
+
+        # We want these to broadcast to the shape
+        #     (idx_underlying_outcome, idx_outcome, idx_modelparam, idx_experiment).
+        # Thus, we need shape
+        #     (idx_underlying_outcome,           1, idx_modelparam,              1).
+        mu = (modelparams[:, self._mu_slice].T)[:, np.newaxis, :, np.newaxis]
+        sigma = np.sqrt(
+            (modelparams[:, self._sigma2_slice].T)[:, np.newaxis, :, np.newaxis]
+        )
+
+        assert np.all(sigma > 0)
+
+        # Now we can rescale the outcomes to be random variates z drawn from N(0, 1).
+        scaled_outcomes = (outcomes[np.newaxis,:,np.newaxis,np.newaxis] - mu) / sigma
+
+        # We can then compute the conditional likelihood Pr(z | underlying_outcome, model).
+        conditional_L = norm(0, 1).pdf(scaled_outcomes)
+
+        # To find the marginalized likeihood, we now need the underlying likelihood
+        # Pr(underlying_outcome | model), so that we can sum over the idx_u_o axis.
+        # Note that we need to add a new axis to shift the underlying outcomes left
+        # of the real-valued outcomes z.
+        underlying_L = self.underlying_likelihood(
+            np.array([0, 1], dtype='uint'),
+            modelparams, expparams
+        )[:, None, :, :]
+
+        # Now we marginalize and return.
+        return (underlying_L * conditional_L).sum(axis=0)
+
+    def simulate_experiment(self, modelparams, expparams, repeat=1):
+        super(GaussianHyperparameterizedModel, self).simulate_experiment(modelparams, expparams)
+
+        # Start by generating a bunch of (0, 1) normalized random variates
+        # that we'll randomly rescale to the right location and shape.
+        zs = np.random.randn(modelparams.shape[0], expparams.shape[0])
+
+        # Next, we sample a bunch of underlying outcomes to figure out
+        # how to rescale everything.
+        underlying_outcomes = self.underlying_model.simulate_experiment(
+            modelparams[:, :-4], expparams, repeat=repeat
+        )
+
+        # We can now rescale zs to obtain the actual outcomes.
+        mu = (modelparams[:, self._mu_slice].T)[:, None, :, None]
+        sigma = np.sqrt(
+            (modelparams[:, self._sigma2_slice].T)[:, None, :, None]
+        )
+        outcomes = (
+            np.where(underlying_outcomes, mu[1], mu[0]) +
+            np.where(underlying_outcomes, sigma[1], sigma[0]) * zs
+        )
+
+        return outcomes[0,0,0] if outcomes.size == 1 else outcomes
+
 class MultinomialModel(DerivedModel):
     """
     Model representing finite numbers of iid samples from another model with 
@@ -680,7 +806,7 @@ def __init__(
             # therefore, we need to add modelparams
             self._has_fixed_covariance = False
             if self._diagonal:
-                self._srw_names = ["\sigma_{{{}}}".format(name) for name in self._rw_names]
+                self._srw_names = [r"\sigma_{{{}}}".format(name) for name in self._rw_names]
                 self._srw_idxs = (underlying_model.n_modelparams + \
                     np.arange(self._n_rw)).astype(np.int)
             else:
@@ -692,9 +818,9 @@ def __init__(
                 for idx1, name1 in enumerate(self._rw_names):
                     for name2 in self._rw_names[:idx1+1]:
                         if name1 == name2:
-                            self._srw_names.append("\sigma_{{{}}}".format(name1))
+                            self._srw_names.append(r"\sigma_{{{}}}".format(name1))
                         else:
-                            self._srw_names.append("\sigma_{{{},{}}}".format(name2,name1))
+                            self._srw_names.append(r"\sigma_{{{},{}}}".format(name2,name1))
         else:
             # In this case the covariance matrix is fixed and fully specified
             self._has_fixed_covariance = True

src/qinfer/test_models.py

@@ -43,6 +43,7 @@
 __all__ = [
     'SimpleInversionModel',
     'SimplePrecessionModel',
+    'UnknownT2Model',
     'CoinModel',
     'NoisyCoinModel',
     'NDieModel'
@@ -54,10 +55,10 @@
 
 import numpy as np
 
-from .utils import binomial_pdf
-
+from .utils import binomial_pdf, decorate_init
 from .abstract_model import FiniteOutcomeModel, DifferentiableModel
-    
+from ._due import due, Doi
+
 ## CLASSES ####################################################################
 
 class SimpleInversionModel(FiniteOutcomeModel, DifferentiableModel):
@@ -137,28 +138,28 @@ def likelihood(self, outcomes, modelparams, expparams):
         # will cause an error.
         pr0 = np.zeros((modelparams.shape[0], expparams.shape[0]))
         pr0[:, :] = np.cos(t * dw / 2) ** 2
-        
+
         # Now we concatenate over outcomes.
         return FiniteOutcomeModel.pr0_to_likelihood_array(outcomes, pr0)
 
     def score(self, outcomes, modelparams, expparams, return_L=False):
         if len(modelparams.shape) == 1:
             modelparams = modelparams[:, np.newaxis]
-            
+
         t = expparams['t']
         dw = modelparams - expparams['w_']
 
         outcomes = outcomes.reshape((outcomes.shape[0], 1, 1))
 
-        arg = dw * t / 2        
+        arg = dw * t / 2
         q = (
             np.power( t / np.tan(arg), outcomes) *
             np.power(-t * np.tan(arg), 1 - outcomes)
         )[np.newaxis, ...]
 
         assert q.ndim == 4
-        
-        
+
+
         if return_L:
             return q, self.likelihood(outcomes, modelparams, expparams)
         else:
@@ -188,15 +189,21 @@ def likelihood(self, outcomes, modelparams, expparams):
         # Pass the expparams to the superclass as a record array.
         new_eps = np.empty(expparams.shape, dtype=super(SimplePrecessionModel, self).expparams_dtype)
         new_eps['w_'] = 0
-        new_eps['t'] = expparams
+        try:
+            new_eps['t'] = expparams
+        except ValueError:
+            new_eps['t'] = expparams['t']
 
         return super(SimplePrecessionModel, self).likelihood(outcomes, modelparams, new_eps)
 
     def score(self, outcomes, modelparams, expparams, return_L=False):
         # Pass the expparams to the superclass as a record array.
         new_eps = np.empty(expparams.shape, dtype=super(SimplePrecessionModel, self).expparams_dtype)
         new_eps['w_'] = 0
-        new_eps['t'] = expparams
+        try:
+            new_eps['t'] = expparams
+        except ValueError:
+            new_eps['t'] = expparams['t']
 
         q = super(SimplePrecessionModel, self).score(outcomes, modelparams, new_eps, return_L=False)
 
@@ -205,6 +212,53 @@ def score(self, outcomes, modelparams, expparams, return_L=False):
         else:
             return q
 
+@decorate_init(
+    due.dcite(
+        Doi('10.1088/1367-2630/14/10/103013'),
+        description='Robust online Hamiltonian learning',
+        tags=['implementation']
+    )
+)
+class UnknownT2Model(FiniteOutcomeModel):
+    """
+    Describes the free evolution of a single qubit prepared in the
+    :math:`\left|+\right\rangle` state under a Hamiltonian
+    :math:`H = \omega \sigma_z / 2` with an unknown :math:`T_2` process,
+    as explored in [GFWC12]_.
+
+    :modelparam omega: The precession frequency :math:`\omega`.
+    :modelparam T2_inv: The decoherence strength :math:`T_2^{-1}`.
+    :scalar-expparam float: The evolution time :math:`t`.
+    """
+
+    @property
+    def n_modelparams(self): return 2
+
+    @property
+    def modelparam_names(self): return [r'\omega', r'T_2^{-1}']
+
+    @property
+    def expparams_dtype(self):
+        return [('t', 'float')]
+
+    def n_outcomes(self, modelparams):
+        return 2
+
+    def are_models_valid(self, modelparams):
+        return np.all(modelparams >= 0, axis=1)
+
+    def likelihood(self, outcomes, modelparams, expparams):
+        w, T2_inv = modelparams.T[:, :, None]
+        t = expparams['t']
+
+        visibility = np.exp(-t * T2_inv)
+
+        pr0 = np.empty((w.shape[0], t.shape[0]))
+        pr0[:, :] = visibility * np.cos(w * t / 2) ** 2 + (1 - visibility) / 2
+
+        return FiniteOutcomeModel.pr0_to_likelihood_array(outcomes, pr0)
+
+
 class CoinModel(FiniteOutcomeModel, DifferentiableModel):
     r"""
     Arguably the simplest possible model; the unknown model parameter