v0.5.1

* fixes to seeding
* IndependentJoint distribution (replaces MixedDistribution)
* Mixture of transformed Gaussians
* updates to generators
* utils/math

delfi/distribution/BaseDistribution.py

@@ -71,6 +71,11 @@ def gen(self, n_samples=1):
         """
         pass
 
+    def reseed(self, seed):
+        """Reseeds the distribution's RNG"""
+        self.rng.seed(seed=seed)
+        self.seed = seed
+
     def gen_newseed(self):
         """Generates a new random seed"""
         if self.seed is None:

delfi/distribution/Gamma.py

@@ -4,7 +4,7 @@
 
 
 class Gamma(BaseDistribution):
-    def __init__(self, alpha=1., beta=1., seed=None):
+    def __init__(self, alpha=1., beta=1., offset=0., seed=None):
         """Univariate (!) Gamma distribution
 
         Parameters
@@ -25,6 +25,7 @@ def __init__(self, alpha=1., beta=1., seed=None):
         assert np.all(beta > 0.), 'Should be greater than zero.'
         self.alpha = alpha
         self.beta = beta
+        self.offset = offset
         self._gamma = gamma(a=alpha, scale=1./beta)
 
     @property
@@ -40,13 +41,13 @@ def std(self):
     @copy_ancestor_docstring
     def eval(self, x, ii=None, log=True):
         # univariate distribution only, i.e. ii=[0] in any case
-        return self._gamma.logpdf(x) if log else self._gamma.pdf(x)
+        return self._gamma.logpdf(x-self.offset) if log else self._gamma.pdf(x-self.offset)
 
     @copy_ancestor_docstring
     def gen(self, n_samples=1, seed=None):
         # See BaseDistribution.py for docstring
 
         x = self.rng.gamma(shape=self.alpha, 
                            scale=1./self.beta, 
-                           size=(n_samples, self.ndim))
+                           size=(n_samples, self.ndim)) + self.offset
         return x

delfi/distribution/IndependentJoint.py

@@ -0,0 +1,91 @@
+import numpy as np
+
+from .BaseDistribution import BaseDistribution
+
+
+class IndependentJoint(BaseDistribution):
+    """Joint distribution composed of statistically independent sub-
+    distributions.
+
+    This class defines a concatenation of a list of distributions.
+    It supports `eval()` and `gen()`.
+
+    Parameters
+    ----------
+    dists : array of distributions
+        Array of distributions
+    seed : int or None
+        If provided, random number generator will be seeded
+    """
+    def __init__(self, dists, seed=None):
+
+        for d in dists:
+            assert not isinstance(d, IndependentJoint), \
+                "IndependentJoint objects cannot be nested"
+        self.dists = [ d for d in dists if d.ndim > 0 ]
+        self.dimlist = [ d.ndim for d in self.dists]
+        ndim = np.sum(self.dimlist)
+
+        # self.dist_index_eachdim stores the index of the child distribution
+        # for each dimension of the full distribution.
+        self.dist_index_eachdim = np.zeros(ndim, dtype=int)
+        # self.ii_full2child stores the index into the child distribution for
+        # each index into the full distribution
+        self.ii_full2child = np.zeros(ndim, dtype=int)
+
+        # list of indices for the full distribution for each child distribution:
+        self.ii_list_eachdist = []
+        csdims = np.append(0, np.cumsum(self.dimlist))
+        for j in range(len(self.dists)):
+            ii_child = np.arange(csdims[j], csdims[j + 1])
+            self.dist_index_eachdim[ii_child] = j
+            self.ii_full2child[ii_child] = np.arange(self.dimlist[j])
+        super().__init__(ndim=ndim, seed=seed)
+
+    def mean(self):
+        return np.concatenate([d.mean for d in self.dists])
+
+    def std(self):
+        return np.concatenate([d.std for d in self.dists])
+
+    def eval(self, x, ii=None, log=True):
+        if ii is not None:
+            ii = np.atleast_1d(ii)
+            if ii.dtype == bool:  # convert to array of indices
+                assert ii.size == self.ndim, 'incorrectly sized binary mask'
+                ii = np.flatnozero(ii)
+            assert (np.diff(ii) > 0).all(), 'ii must be increasing'
+            assert (ii >= 0).all() and (ii < self.ndim).all(), 'invalid index'
+
+            xsplit, ds, ii_children = [], [], []
+            for j, d in enumerate(self.dists):
+                # is_childj stores whether each element of ii is in child j
+                is_childj = self.dist_index_eachdim[ii] == j
+                if not is_childj.any():
+                    continue
+                ds.append(d)
+                # ii_child stores indices into child distribution's dimensions
+                ii_child = self.ii_full2child[ii[is_childj]]
+                xsplit.append(x[:, is_childj])
+                if ii_child.size == self.dimlist[j]:
+                    ii_children.append(None)  # full child distribution
+                else:
+                    ii_children.append(ii_child)
+        else:
+            ds = self.dists
+            xsplit = np.split(x, np.cumsum(self.dimlist), axis=-1)
+            ii_children = [None for d in ds]  # full child distribution
+
+        logps = [d.eval(x, ii=ii, log=True)
+                 for d, x, ii in zip(ds, xsplit, ii_children)]
+        # each element of logps is a vector with one log prob per data point
+        logp = np.sum(np.vstack(logps), axis=0)
+        return logp if log else np.exp(logp)
+
+    def gen(self, n_samples=1):
+        return np.concatenate([ d.gen(n_samples) for d in self.dists ], axis=-1)
+
+    def reseed(self, seed):
+        super().reseed(seed)
+        for d in self.dists:
+            d.reseed(self.gen_newseed())

delfi/distribution/TransformedNormal.py

@@ -31,7 +31,7 @@ def __init__(self, m=None, P=None, U=None, S=None, Pm=None,
             List of flags for each variable whether it is untransformed (=0),
             log-transformed (=1) or logit-transformed (=2).
         lower: list or np.array, 1d
-            Lower bounds for logit-box. Defaults to 1 for each parameter.
+            Lower bounds for logit-box. Defaults to 0 for each parameter.
         upper: list or np.array, 1d
             Upper bounds for logit-box. Defaults to 1 for each parameter.
         seed : int or None

delfi/distribution/__init__.py

@@ -2,9 +2,9 @@
 from delfi.distribution.Gaussian import Gaussian
 from delfi.distribution.TransformedNormal import TransformedNormal
 from delfi.distribution.StudentsT import StudentsT
-from delfi.distribution.Uniform import Uniform, LogUniform
-
-from delfi.distribution.MixedDistribution import MixedDistribution
+from delfi.distribution.Uniform import Uniform
+from delfi.distribution.Gamma import Gamma
+from delfi.distribution.IndependentJoint import IndependentJoint
 
 from delfi.distribution.mixture.EllipsoidalMixture import MoE
 from delfi.distribution.mixture.GaussianMixture import MoG

delfi/distribution/mixture/BaseMixture.py

@@ -80,6 +80,20 @@ def gen_comp(self, n_samples):
         """Generate component index according to self.a"""
         return self.discrete_sample.gen(n_samples).reshape(-1)  # n_samples,
 
+    def reseed(self, seed):
+        """Reseeds the following RNGs in the following order:
+        1) Master RNG for the mixture object, using the input seed
+        2) RNG for the discrete distribution used to sample components. The seed
+        is generated using the master RNG.
+        3) RNG for each mixture component, in order. Each seed is generated by
+        the master RNG.
+        """
+        self.rng.seed(seed=seed)
+        self.seed = seed
+        self.discrete_sample.reseed(seed=self.gen_newseed())
+        for x in self.xs:
+            x.reseed(seed=self.gen_newseed())
+
     def gen_newseed(self):
         """Generates a new random seed"""
         if self.seed is None:
@@ -115,4 +129,3 @@ def prune_negligible_components(self, threshold):
         self.a = np.delete(self.a, ii)
         self.a += total_del_a / self.n_components
         self.xs = [x for i, x in enumerate(self.xs) if i not in ii]
-

delfi/distribution/mixture/GaussianMixture.py

@@ -47,7 +47,8 @@ def __init__(
             super().__init__(
                 a=np.asarray(a),
                 ncomp=len(ms),
-                ndim=len(ms[0]),
+                ndim=np.asarray(
+                    ms[0]).ndim,
                 seed=seed)
 
             if Ps is not None: