mackelab/delfi

using uniform priors in a moG results in error

coschroeder opened this issue · 2 comments

I want to use SNPE in a more involved problem, and as a first step I constructed this toy problem to fit a multidimensional gaussian.
If I use a gaussian prior everything works quite well, but switching to an uniform prior, defined via the MixedDistribution class, I get an error message at the beginning of the training of round 2.


ValueError Traceback (most recent call last)
in ()
54 # define and run SNPE
55 inf_snpe = SNPE(generator=g, n_components=1, n_hiddens=[10], obs=xo, pilot_samples=100)
---> 56 logs, tds, posteriors = inf_snpe.run(n_train=[1000, 1000], n_rounds=10, stop_on_nan=True)
57 posterior = posteriors[-1]

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/inference/SNPE.py in run(self, n_train, n_rounds, epochs, minibatch, round_cl, stop_on_nan, proposal, monitor, **kwargs)
180 verbose = '(round {}) '.format(self.round) if self.verbose else False
181
--> 182 trn_data = self.gen(n_train_round, prior_mixin=self.prior_mixin, verbose=verbose)
183 n_train_round = trn_data[0].shape[0]
184

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/inference/BaseInference.py in gen(self, n_samples, n_reps, prior_mixin, verbose)
108 """
109 verbose = self.verbose if verbose is None else verbose
--> 110 params, stats = self.generator.gen(n_samples, prior_mixin=prior_mixin, verbose=verbose)
111
112 # z-transform params and stats

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/generator/BaseGenerator.py in gen(self, n_samples, n_reps, skip_feedback, prior_mixin, minibatch, keep_data, verbose)
106 skip_feedback=skip_feedback,
107 prior_mixin=prior_mixin,
--> 108 verbose = verbose)
109
110 # Run forward model for params (in batches)

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/generator/BaseGenerator.py in draw_params(self, n_samples, skip_feedback, prior_mixin, verbose)
52 proposed_param = self.prior.gen(n_samples=1) # dim params,
53 else:
---> 54 proposed_param = self.proposal.gen(n_samples=1)
55
56 # check if parameter vector is valid

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/distribution/mixture/StudentsTMixture.py in gen(self, n_samples)
64
65 ns = [np.sum((ii == i).astype(int)) for i in range(self.n_components)]
---> 66 samples = [x.gen(n) for x, n in zip(self.xs, ns)]
67 samples = np.concatenate(samples, axis=0)
68 self.rng.shuffle(samples)

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/distribution/mixture/StudentsTMixture.py in (.0)
64
65 ns = [np.sum((ii == i).astype(int)) for i in range(self.n_components)]
---> 66 samples = [x.gen(n) for x, n in zip(self.xs, ns)]
67 samples = np.concatenate(samples, axis=0)
68 self.rng.shuffle(samples)

~/Documents/Ribbon_PR_project/Fitting_New/fitting_git/DELFI/DELFI_package/delfi-master/delfi/distribution/StudentsT.py in gen(self, n_samples)
65 u = self.rng.chisquare(self.dof, n_samples) / self.dof
66 y = self.rng.multivariate_normal(np.zeros(self.ndim),
---> 67 self.S, (n_samples,))
68 return self.m + y / np.sqrt(u)[:, None]

mtrand.pyx in mtrand.RandomState.multivariate_normal()

/usr/local/lib/python3.5/dist-packages/scipy/linalg/decomp_svd.py in svd(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)
107
108 """
--> 109 a1 = _asarray_validated(a, check_finite=check_finite)
110 if len(a1.shape) != 2:
111 raise ValueError('expected matrix')

/usr/local/lib/python3.5/dist-packages/scipy/_lib/_util.py in _asarray_validated(a, check_finite, sparse_ok, objects_ok, mask_ok, as_inexact)
236 raise ValueError('masked arrays are not supported')
237 toarray = np.asarray_chkfinite if check_finite else np.asarray
--> 238 a = toarray(a)
239 if not objects_ok:
240 if a.dtype is np.dtype('O'):

/usr/local/lib/python3.5/dist-packages/numpy/lib/function_base.py in asarray_chkfinite(a, dtype, order)
459 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all():
460 raise ValueError(
--> 461 "array must not contain infs or NaNs")
462 return a
463

ValueError: array must not contain infs or NaNs

I guess it is because of NaNs in the Cov matrix of the MoG I want to fit.
But am I doing something wrong? Or is this due to another problem I am not aware?

Here is my code to reproduce the error:


import delfi
import delfi.distribution as dd
from delfi.inference import SNPE
from delfi.simulator.BaseSimulator import BaseSimulator
from delfi.summarystats import Identity
from delfi.generator import Default
from delfi.distribution import MixedDistribution


#  define simulator

class MultiDGauss(BaseSimulator):
    """
    draw a random number from a multidim gaussian with Id Cov matrix
    ----------
    params: 1dim array, mean of the distribution marginals
    seed :  int or None
            If set, randomness is seeded
    """
    
    def __init__(self, dim=1, seed=None):   
        super().__init__(dim_param=dim, seed=seed)
        self.dim = dim
        
        
    def gen_single(self, params):
        """
        draw one sample
        ------
        params : 1 dim array
        """        
        params = np.asarray(params).reshape(-1)
        assert params.ndim == 1
        assert params.shape[0] == self.dim_param

        sample = scp.stats.norm.rvs(loc=[params],size=self.dim)

        return {'data': sample.reshape(-1)}

# define model
dimensions = 2
m = MultiDGauss(dim =dimensions)

# define prior
pgauss = delfi.distribution.Gaussian(m=np.zeros(dimensions), S=np.eye(dimensions)*5)
puniform = MixedDistribution ([dd.Uniform(lower=[-5], upper=[5])]* dimensions)


# define summary statistics
s = Identity()

# initialize 
# choose prior
prior = puniform #pgauss
g = Default(model=m, prior=prior, summary=s)

# 'observed' value
xo = np.array([[-1,2]])

# check if data has right dimension
if np.shape(xo)[1] != dimensions:
    print('observed data does not match the dimension of parameters.')
else:
    # define and run SNPE
    inf_snpe = SNPE(generator=g, n_components=1, n_hiddens=[10], obs=xo, pilot_samples=100)
    logs, tds, posteriors = inf_snpe.run(n_train=[1000, 1000], n_rounds=10, stop_on_nan=True)
    posterior = posteriors[-1]`
```

Seems to be an issue with linear vs. log domain in the eval method of the mixeddistribution class. Will be fixed shortly.

In the meantime, the uniform distribution class can be defined using vector upper/lower bounds to do what you're going for here. When I change the line in your code defining the prior to

puniform = dd.Uniform(lower= -5.0 * np.ones(dimensions), upper= 5.0 * np.ones(dimensions))

The code runs without any errors.

Aside from this, I would recommend trying not to copy distributions or other objects using the python list multiply operator, as in

[dd.Uniform(lower=[-5], upper=[5])]* dimensions

This will create a list of 5 references to the same uniform distribution object. This can introduce spurious correlations due to a shared RNG, and break RNG seeding.

If you want to copy delfi objects without coupling their RNGs, copy.deepcopy is generally safe.

Thanks, that snippet fixed the problem.
Also thanks for your advise not copying lists.