Python scipy.stats.multivariate_normal() Examples

def norm_post_sim(modes,cov_matrix):
    post = multivariate_normal(modes,cov_matrix)
    nsims = 30000
    phi = np.zeros([nsims,len(modes)])

    for i in range(0,nsims):
         phi[i] = post.rvs()

    chain = np.array([phi[i][0] for i in range(len(phi))])
    for m in range(1,len(modes)):
        chain = np.vstack((chain,[phi[i][m] for i in range(len(phi))]))

    mean_est = [np.mean(np.array([phi[i][j] for i in range(len(phi))])) for j in range(len(modes))]
    median_est = [np.median(np.array([phi[i][j] for i in range(len(phi))])) for j in range(len(modes))]
    upper_95_est = [np.percentile(np.array([phi[i][j] for i in range(len(phi))]),95) for j in range(len(modes))]
    lower_95_est = [np.percentile(np.array([phi[i][j] for i in range(len(phi))]),5) for j in range(len(modes))]

    return chain, mean_est, median_est, upper_95_est, lower_95_est 

def setUp_configure(self):
        from scipy import stats
        self.dist = distributions.MultivariateNormal
        self.scipy_dist = stats.multivariate_normal
        self.scipy_onebyone = True
        d, = self.event_shape

        self.test_targets = set([
            'batch_shape', 'entropy', 'event_shape', 'log_prob',
            'support'])

        loc = numpy.random.uniform(
            -1, 1, self.shape + self.event_shape).astype(numpy.float32)
        cov = numpy.random.normal(
            size=(numpy.prod(self.shape),) + (d, d))
        cov = [cov_.dot(cov_.T) for cov_ in cov]
        cov = numpy.vstack(cov).reshape(self.shape + (d, d))
        scale_tril = numpy.linalg.cholesky(cov).astype(numpy.float32)
        self.params = {'loc': loc, 'scale_tril': scale_tril}
        self.scipy_params = {'mean': loc, 'cov': cov} 

def sample(transition_matrix,
           means, covs,
           start_state, n_samples,
           random_state):
    n_states, n_features, _ = covs.shape
    states = np.zeros(n_samples, dtype='int')
    emissions = np.zeros((n_samples, n_features))
    for i in range(n_samples):
        if i == 0:
            prev_state = start_state
        else:
            prev_state = states[i - 1]
        state = random_state.choice(n_states,
                                    p=transition_matrix[:, prev_state])
        emissions[i] = random_state.multivariate_normal(means[state],
                                                        covs[state])
        states[i] = state
    return emissions, states 

def precompute_marginals(self):
        sys.stderr.write('Precomputing marginals...\n')
        self._pdfs = [None] * self._num_instances
        # precomputing all possible marginals
        for i in xrange(self._num_instances):
            mean = self._corrected_means[i]
            cov = self._corrected_covs[i]
            self._pdfs[i] = [None] * (2 ** mean.shape[0])
            for marginal_pattern in itertools.product([False, True], repeat=mean.shape[0]):
                marginal_length = marginal_pattern.count(True)
                if marginal_length == 0:
                    continue
                m = np.array(marginal_pattern)
                marginal_mean = mean[m]
                mm = m[:, np.newaxis]
                marginal_cov = cov[np.dot(mm, mm.transpose())].reshape((marginal_length, marginal_length))
                self._pdfs[i][hash_bool_array(m)] = multivariate_normal(mean=marginal_mean, cov=marginal_cov) 

def _testPDFShapes(self, mvn_dist, mu, sigma):
    with self.test_session() as sess:
      mvn = mvn_dist(mu, sigma)
      x = 2 * tf.ones_like(mu)

      log_pdf = mvn.log_pdf(x)
      pdf = mvn.pdf(x)

      mu_value = np.ones([3, 3, 2])
      sigma_value = np.zeros([3, 3, 2, 2])
      sigma_value[:] = np.identity(2)
      x_value = 2. * np.ones([3, 3, 2])
      feed_dict = {mu: mu_value, sigma: sigma_value}

      scipy_mvn = stats.multivariate_normal(mean=mu_value[(0, 0)],
                                            cov=sigma_value[(0, 0)])
      expected_log_pdf = scipy_mvn.logpdf(x_value[(0, 0)])
      expected_pdf = scipy_mvn.pdf(x_value[(0, 0)])

      log_pdf_evaled, pdf_evaled = sess.run([log_pdf, pdf], feed_dict=feed_dict)
      self.assertAllEqual([3, 3], log_pdf_evaled.shape)
      self.assertAllEqual([3, 3], pdf_evaled.shape)
      self.assertAllClose(expected_log_pdf, log_pdf_evaled[0, 0])
      self.assertAllClose(expected_pdf, pdf_evaled[0, 0]) 

def testLogPDFScalarBatch(self):
    with self.test_session():
      mu = self._rng.rand(2)
      chol, sigma = self._random_chol(2, 2)
      mvn = distributions.MultivariateNormalCholesky(mu, chol)
      x = self._rng.rand(2)

      log_pdf = mvn.log_pdf(x)
      pdf = mvn.pdf(x)

      scipy_mvn = stats.multivariate_normal(mean=mu, cov=sigma)

      expected_log_pdf = scipy_mvn.logpdf(x)
      expected_pdf = scipy_mvn.pdf(x)
      self.assertEqual((), log_pdf.get_shape())
      self.assertEqual((), pdf.get_shape())
      self.assertAllClose(expected_log_pdf, log_pdf.eval())
      self.assertAllClose(expected_pdf, pdf.eval()) 

def testLogPDFXIsHigherRank(self):
    with self.test_session():
      mu = self._rng.rand(2)
      chol, sigma = self._random_chol(2, 2)
      mvn = distributions.MultivariateNormalCholesky(mu, chol)
      x = self._rng.rand(3, 2)

      log_pdf = mvn.log_pdf(x)
      pdf = mvn.pdf(x)

      scipy_mvn = stats.multivariate_normal(mean=mu, cov=sigma)

      expected_log_pdf = scipy_mvn.logpdf(x)
      expected_pdf = scipy_mvn.pdf(x)
      self.assertEqual((3,), log_pdf.get_shape())
      self.assertEqual((3,), pdf.get_shape())
      self.assertAllClose(expected_log_pdf, log_pdf.eval())
      self.assertAllClose(expected_pdf, pdf.eval()) 

def testLogPDFXLowerDimension(self):
    with self.test_session():
      mu = self._rng.rand(3, 2)
      chol, sigma = self._random_chol(3, 2, 2)
      mvn = distributions.MultivariateNormalCholesky(mu, chol)
      x = self._rng.rand(2)

      log_pdf = mvn.log_pdf(x)
      pdf = mvn.pdf(x)

      self.assertEqual((3,), log_pdf.get_shape())
      self.assertEqual((3,), pdf.get_shape())

      # scipy can't do batches, so just test one of them.
      scipy_mvn = stats.multivariate_normal(mean=mu[1, :], cov=sigma[1, :, :])
      expected_log_pdf = scipy_mvn.logpdf(x)
      expected_pdf = scipy_mvn.pdf(x)

      self.assertAllClose(expected_log_pdf, log_pdf.eval()[1])
      self.assertAllClose(expected_pdf, pdf.eval()[1]) 

def testSampleWithSampleShape(self):
    with self.test_session():
      mu = self._rng.rand(3, 5, 2)
      chol, sigma = self._random_chol(3, 5, 2, 2)

      mvn = distributions.MultivariateNormalCholesky(mu, chol)
      samples_val = mvn.sample((10, 11, 12), seed=137).eval()

      # Check sample shape
      self.assertEqual((10, 11, 12, 3, 5, 2), samples_val.shape)

      # Check sample means
      x = samples_val[:, :, :, 1, 1, :]
      self.assertAllClose(
          x.reshape(10 * 11 * 12, 2).mean(axis=0),
          mu[1, 1], atol=1e-2)

      # Check that log_prob(samples) works
      log_prob_val = mvn.log_prob(samples_val).eval()
      x_log_pdf = log_prob_val[:, :, :, 1, 1]
      expected_log_pdf = stats.multivariate_normal(
          mean=mu[1, 1, :], cov=sigma[1, 1, :, :]).logpdf(x)
      self.assertAllClose(expected_log_pdf, x_log_pdf) 

def draw(X, pred, means, covariances, output):
    xp = cupy.get_array_module(X)
    for i in range(2):
        labels = X[pred == i]
        if xp is cupy:
            labels = labels.get()
        plt.scatter(labels[:, 0], labels[:, 1], c=np.random.rand(1, 3))
    if xp is cupy:
        means = means.get()
        covariances = covariances.get()
    plt.scatter(means[:, 0], means[:, 1], s=120, marker='s', facecolors='y',
                edgecolors='k')
    x = np.linspace(-5, 5, 1000)
    y = np.linspace(-5, 5, 1000)
    X, Y = np.meshgrid(x, y)
    for i in range(2):
        dist = stats.multivariate_normal(means[i], covariances[i])
        Z = dist.pdf(np.stack([X, Y], axis=-1))
        plt.contour(X, Y, Z)
    plt.savefig(output) 

def generate_dependent_normal_samples(n, mu, mcorr, dvar):
    d = mcorr.shape[1]
    mvar = np.copy(mcorr)
    np.fill_diagonal(mvar, 1.)
    # logger.debug("mvar:\n%s" % str(mvar))
    if d > 1:
        for i in range(0, d-1):
            for j in range(i+1, d):
                mvar[j, i] = mvar[i, j]
        p = np.diag(np.sqrt(dvar))
        mvar = p.dot(mvar).dot(p)
        rv = mvn(mu, mvar)
        s = rv.rvs(size=n)
    else:
        s = np.random.normal(loc=mu[0], scale=np.sqrt(dvar[0]), size=n)
        # logger.debug(str(list(s)))
        s = np.reshape(s, (-1, 1))

    if n == 1:
        # handle the case where numpy automatically makes column vector if #rows is 1
        s = np.reshape(s, (n, len(s)))
    return s 

def __init__(self, mean=2, cov=None, ndim_x=1, ndim_y=1, has_cdf=True, has_pdf=True, can_sample=True):
    self.ndim_x = ndim_x
    self.ndim_y = ndim_y
    self.ndim = self.ndim_x + self.ndim_y

    self.mean = mean
    # check if mean is scalar
    if isinstance(self.mean, list):
      self.mean = np.array(self.ndim_y * [self.mean])

    self.cov = cov
    if self.cov is None:
      self.cov = np.identity(self.ndim_y)
    assert self.cov.shape[0] == self.cov.shape[1] == self.ndim_y

    # safe data stats
    self.y_mean = self.mean
    self.y_std = np.sqrt(np.diagonal(self.cov))

    self.gaussian = stats.multivariate_normal(mean=self.mean, cov=self.cov)
    self.fitted = False

    self.can_sample = can_sample
    self.has_pdf = has_pdf
    self.has_cdf = has_cdf 

def __init__(self, mean=2, cov=None, ndim_x=1, ndim_y=1, has_cdf=True, has_pdf=True, can_sample=True):
    self.ndim_x = ndim_x
    self.ndim_y = ndim_y
    self.ndim = self.ndim_x + self.ndim_y

    self.mean = mean
    self.cov = cov
    # check if mean is scalar
    if isinstance(self.mean, list):
      self.mean = np.array(self.ndim_y*[self.mean])

    if self.cov is None:
      self.cov = np.identity(self.ndim_y)

    self.gaussian = stats.multivariate_normal(mean=self.mean, cov=self.cov)
    self.fitted = False

    self.can_sample = can_sample
    self.has_pdf = has_pdf
    self.has_cdf = has_cdf 

def _build_model(self, X, Y):
    # save mean and std of data for normalization
    self.x_std = np.std(X, axis=0)
    self.x_mean = np.mean(X, axis=0)
    self.y_mean = np.std(Y, axis=0)
    self.y_std = np.std(Y, axis=0)

    self.n_train_points = X.shape[0]

    # lazy learner - just store training data
    self.X_train = self._normalize_x(X)
    self.Y_train = Y

    # prepare Gaussians centered in the Y points
    self.locs_array = np.vsplit(Y, self.n_train_points)
    self.log_kernel = multivariate_normal(mean=np.ones(self.ndim_y)).logpdf

    # select / properly initialize bandwidth and epsilon
    if isinstance(self.bandwidth, (int, float)):
      self.bandwidth = self.y_std * self.bandwidth

    if self.param_selection == 'normal_reference':
      self.bandwidth = self._normal_reference()
    elif self.param_selection == 'cv_ml':
      self.bandwidth, self.epsilon = self._cv_ml() 

def _build_model(self, X, Y):
    # save mean and variance of data for normalization
    self.x_mean, self.y_mean = np.mean(X, axis=0), np.mean(Y, axis=0)
    self.x_std, self.y_std = np.std(X, axis=0),  np.std(Y, axis=0)

    # get locations of the gaussian kernel centers
    if self.center_sampling_method == 'all':
      self.n_centers = X.shape[0]
    else:
      self.n_centers = min(self.n_centers, X.shape[0])

    n_locs = self.n_centers
    X_Y_normalized = np.concatenate(list(self._normalize(X, Y)), axis=1)
    centroids = sample_center_points(X_Y_normalized, method=self.center_sampling_method, k=n_locs,
                                     keep_edges=self.keep_edges, random_state=self.random_state)
    self.centr_x = centroids[:, 0:self.ndim_x]
    self.centr_y = centroids[:, self.ndim_x:]

    #prepare gaussians for sampling
    self.gaussians_y = [stats.multivariate_normal(mean=center, cov=self.bandwidth) for center in self.centr_y]

    assert self.centr_x.shape == (n_locs, self.ndim_x) and self.centr_y.shape == (n_locs, self.ndim_y) 

def forward_simulate(self, input_values, k, rng=np.random.RandomState(), mpi_comm=None):
        """
        Samples from a multivariate normal distribution using the current values for each probabilistic model from which the
        model derives.

        Parameters
        ----------
        input_values: list
            List of input parameters, in the same order as specified in the InputConnector passed to the init function
        k: integer
            The number of samples that should be drawn.
        rng: Random number generator
            Defines the random number generator to be used. The default value uses a random seed to initialize the generator.

        Returns
        -------
        list: [np.ndarray]
            A list containing the sampled values as np-array.
        """

        dim = self.get_output_dimension()
        mean = np.array(input_values[0:dim])
        cov = np.array(input_values[dim:dim+dim**2]).reshape((dim, dim))
        result = rng.multivariate_normal(mean, cov, k)
        return [np.array([result[i,:]]).reshape(-1,) for i in range(k)] 

def gaussian_to_pi_density(psi_mesh, mu, Sigma):
    pi_mesh = np.array(list(map(psi_to_pi, psi_mesh)))
    valid_pi = np.all(np.isfinite(pi_mesh), axis=1)
    pi_mesh = pi_mesh[valid_pi,:]

    # Compute the det of the Jacobian of the inverse mapping
    det_jacobian = np.array(list(map(det_jacobian_pi_to_psi, pi_mesh)))
    det_jacobian = det_jacobian[valid_pi]

    # Compute the multivariate Gaussian density
    # pi_pdf = np.exp(log_dirichlet_density(pi_mesh, alpha=alpha))
    from scipy.stats import multivariate_normal
    psi_dist = multivariate_normal(mu, Sigma)
    psi_pdf = psi_dist.pdf(psi_mesh)
    psi_pdf = psi_pdf[valid_pi]

    # The psi density is scaled by the det of the Jacobian
    pi_pdf = psi_pdf * det_jacobian

    return pi_mesh, pi_pdf 

def test_calc_logp_mariginal_likelihood():
    np.random.seed(1234)

    n_k = 100
    K = 5
    dim = 10

    data_dictionary = generate_data(n_k, K, dim)
    X = data_dictionary['data']
    Y = data_dictionary['labels']
    model = plda.Model(X, Y)

    prior_mean = model.prior_params['mean']
    prior_cov_diag = model.prior_params['cov_diag']

    logpdf = gaussian(prior_mean, np.diag(prior_cov_diag + 1)).logpdf

    data = np.random.random((n_k, prior_mean.shape[-1]))
    expected_logps = logpdf(data)
    actual_logps = model.calc_logp_marginal_likelihood(data[:, None])

    assert_allclose(actual_logps, expected_logps) 

def illustrate_preprocessing():
    x = np.random.multivariate_normal(np.array([5.0,5.0]),
            np.array([[5.0,3.0],[3.0,4.0]]),size=1000)
    x_demean = x - np.mean(x, axis=0)
    x_unitsd = x_demean/(np.std(x_demean,axis=0))
    x_whiten = np.dot(x_demean, whitening_matrix(x_demean))

    fig = pl.figure(figsize=(10,10))
    
    def mk_subplot(n, data, label):
        ax = fig.add_subplot(2,2,n)
        ax.scatter(data[:,0], data[:,1])
        ax.set_xlim((-10,10))
        ax.set_ylim((-10,10))
        ax.set_xlabel(label)

    mk_subplot(1, x, "Original")
    mk_subplot(2, x_demean, "De-meaned")
    mk_subplot(3, x_unitsd, "Unit SD")
    mk_subplot(4, x_whiten, "Whitened")
    pl.show() 

def test_calc_logp_prior_():
    def calc_error(truth_dict):
        model = plda.Model(truth_dict['data'], truth_dict['labels'])

        Phi_b = truth_dict['Phi_b']
        prior_mean = truth_dict['prior_mean']
        dim = prior_mean.shape[0]

        random_vectors = np.random.randint(-100, 100, (10, dim))
        expected = gaussian(prior_mean, Phi_b).logpdf(random_vectors)

        latent_vectors = model.transform(random_vectors, 'D', 'U_model')
        predicted = model.calc_logp_prior(latent_vectors)

        error = calc_mean_squared_error(expected, predicted, as_log=True)
        print(error)

        return error

    ks = [10, 100, 1000]  # List of numbers of categories.

    np.random.seed(1234)
    assert_error_falls_as_K_increases(calc_error, n_k=1000, D=2, k_list=ks)
    assert_error_falls_as_K_increases(calc_error, n_k=1000, D=50, k_list=ks) 

def gmm_sample(self, mean=None, w=None, N=10000,n=10,d=2,seed=10):
        np.random.seed(seed)
        self.d = d
        if mean is None:
            mean = np.random.randn(n,d)*10
        if w is None:
            w = np.random.rand(n)
        w = old_div(w,sum(w))
        multi = np.random.multinomial(N,w)
        X = np.zeros((N,d))
        base = 0
        for i in range(n):
            X[base:base+multi[i],:] = np.random.multivariate_normal(mean[i,:], np.eye(self.d), multi[i])
            base += multi[i]
        
        llh = np.zeros(N)
        for i in range(n):
            llh += w[i] * stats.multivariate_normal.pdf(X, mean[i,:], np.eye(self.d))
        #llh = llh/sum(llh)
        return X, llh 

def differential_entropies(X, labels):
    n_samples, n_features = X.shape
    
    labels = np.array(labels)
    names = sorted(set(labels))

    entropies = []
    
    for name in names:
        name_idx = np.where(labels == name)[0]

        gm = GaussianMixture().fit(X[name_idx, :])

        mn = multivariate_normal(
            mean=gm.means_.flatten(),
            cov=gm.covariances_.reshape(n_features, n_features)
        )

        entropies.append(mn.entropy())

    probs = softmax(entropies)

    for name, entropy, prob in zip(names, entropies, probs):
        #print('{}\t{}\t{}'.format(name, entropy, prob))
        print('{}\t{}'.format(name, entropy)) 

def get_pdfs(self):
        ''' Return the pdf of all the emission probabilities '''
        return [multivariate_normal(self.mu[k], np.diag(self.Sigma[k])) for k in range(self.K)] 

def get_pdfs(self):
        ''' Return the pdf of all the emission probabilities '''
        return [multivariate_normal(self.mu[k], self.Sigma[k]) for k in range(self.K)] 

def pdf(self, accepted_parameters_manager, kernel_index, mean, x):
        """Calculates the pdf of the kernel.
        Commonly used to calculate weights.

        Parameters
        ----------
        accepted_parameters_manager: abcpy.AcceptedParametersManager object
            The AcceptedParametersManager to be used.
        kernel_index: integer
            The index of the kernel in the list of kernels in the joint kernel.
        index: integer
            The row to be considered in the accepted_parameters_bds matrix.
        x: The point at which the pdf should be evaluated.

        Returns
        -------
        float
            The pdf evaluated at point x.
        """

        if isinstance(mean[0], (np.float, np.float32, np.float64, np.int, np.int32, np.int64)):
            mean = np.array(mean).astype(float)
            cov = np.array(accepted_parameters_manager.accepted_cov_mats_bds.value()[kernel_index]).astype(float)
            return multivariate_normal(mean, cov, allow_singular=True).pdf(np.array(x).astype(float))
        else:
            mean = np.array(np.concatenate(mean)).astype(float)
            cov = np.array(accepted_parameters_manager.accepted_cov_mats_bds.value()[kernel_index]).astype(float)
            return multivariate_normal(mean, cov, allow_singular=True).pdf(np.concatenate(x)) 

def generate_deformation_field_gauss(shape, points, max_deform=DEFORMATION_MAX,
                                     deform_smooth=DEFORMATION_SMOOTH):
    """ generate deformation field as combination of positive and
    negative Galatians densities scaled in range +/- max_deform

    :param tuple(int,int) shape: tuple of size 2
    :param points: <nb_points, 2> list of landmarks
    :param float max_deform: maximal deformation distance in any direction
    :param float deform_smooth: smoothing the deformation by Gaussian filter
    :return: np.array<shape>
    """
    ndim = len(shape)
    x, y = np.mgrid[0:shape[0], 0:shape[1]]
    pos_grid = np.rollaxis(np.array([x, y]), 0, 3)
    # initialise the deformation
    deform = np.zeros(shape)
    for point in points:
        sign = np.random.choice([-1, 1])
        cov = np.random.random((ndim, ndim))
        cov[np.eye(ndim, dtype=bool)] = 100 * np.random.random(ndim)
        # obtain a positive semi-definite matrix
        cov = np.dot(cov, cov.T) * (0.1 * np.mean(shape))
        gauss = stats.multivariate_normal(point, cov)
        deform += sign * gauss.pdf(pos_grid)
    # normalise the deformation and multiply by the amplitude
    deform *= max_deform / np.abs(deform).max()
    # set boundary region to zeros
    fix_deform_bounds = DEFORMATION_BOUNDARY_COEF * deform_smooth
    deform[:fix_deform_bounds, :] = 0
    deform[-fix_deform_bounds:, :] = 0
    deform[:, :fix_deform_bounds] = 0
    deform[:, -fix_deform_bounds:] = 0
    # smooth the deformation field
    deform = ndimage.gaussian_filter(deform, sigma=deform_smooth, order=0)
    return deform 

def forward_simulate(self, input_values, k, rng=np.random.RandomState(), mpi_comm=None):
        """
        Samples from a multivariate Student's T-distribution using the current values for each probabilistic model from
        which the model derives.

        Parameters
        ----------
        input_values: list
            List of input parameters, in the same order as specified in the InputConnector passed to the init function
        k: integer
            The number of samples that should be drawn.
        rng: Random number generator
            Defines the random number generator to be used. The default value uses a random seed to initialize the
            generator.

        Returns
        -------
        list: [np.ndarray]
            A list containing the sampled values as np-array.
        """

        # Extract input_parameters
        dim = self.get_output_dimension()
        mean = np.array(input_values[0:dim])
        cov = np.array(input_values[dim:dim+dim**2]).reshape((dim, dim))
        df = input_values[-1]

        if (df == np.inf):
            chisq = 1.0
        else:
            chisq = rng.chisquare(df, k) / df
            chisq = chisq.reshape(-1, 1).repeat(dim, axis=1)
        mvn = rng.multivariate_normal(np.zeros(dim), cov, k)
        result = (mean + np.divide(mvn, np.sqrt(chisq)))
        return [np.array([result[i, :]]).reshape(-1, ) for i in range(k)] 

def genGTmap(h, w, pos_x, pos_y, sigma_h=1, sigma_w=1, init=None):
    """
    Compute the heat-map of size (w x h) with a gaussian distribution fit in
    position (pos_x, pos_y) and a convariance matix defined by the related
    sigma values.
    The resulting heat-map can be summed to a given heat-map init.
    """
    if pos_x>0 and pos_y>0:
        init = init if init is not None else []

        cov_matrix = np.eye(2) * ([sigma_h**2, sigma_w**2])

        x, y = np.mgrid[0:h, 0:w]
        pos = np.dstack((x, y))
        rv = multivariate_normal([pos_x, pos_y], cov_matrix)

        tmp = rv.pdf(pos)
        hmap = np.multiply(
            tmp, np.sqrt(np.power(2 * np.pi, 2) * np.linalg.det(cov_matrix))
        )
        idx = np.where(hmap.flatten() <= np.exp(-4.6052))
        hmap.flatten()[idx] = 0

        if np.size(init) == 0:
            return hmap

        assert (np.shape(init) == hmap.shape)
        hmap += init
        idx = np.where(hmap.flatten() > 1)
        hmap.flatten()[idx] = 1
        return hmap
    else:
        return np.zeros((h, w)) 

def gaussian_heatmap(h, w, pos_x, pos_y, sigma_h=1, sigma_w=1, init=None):
    """
    Compute the heat-map of size (w x h) with a gaussian distribution fit in
    position (pos_x, pos_y) and a convariance matix defined by the related
    sigma values.
    The resulting heat-map can be summed to a given heat-map init.
    """
    init = init if init is not None else []

    cov_matrix = np.eye(2) * ([sigma_h**2, sigma_w**2])

    x, y = np.mgrid[0:h, 0:w]
    pos = np.dstack((x, y))
    rv = multivariate_normal([pos_x, pos_y], cov_matrix)

    tmp = rv.pdf(pos)
    hmap = np.multiply(
        tmp, np.sqrt(np.power(2 * np.pi, 2) * np.linalg.det(cov_matrix))
    )
    idx = np.where(hmap.flatten() <= np.exp(-4.6052))
    hmap.flatten()[idx] = 0

    if np.size(init) == 0:
        return hmap

    assert (np.shape(init) == hmap.shape)
    hmap += init
    idx = np.where(hmap.flatten() > 1)
    hmap.flatten()[idx] = 1
    return hmap 

def update(self, accepted_parameters_manager, kernel_index, row_index, rng=np.random.RandomState()):
        """Updates the parameter values contained in the accepted_paramters_manager using a multivariate normal distribution.

        Parameters
        ----------
        accepted_parameters_manager: abcpy.AcceptedParametersManager object
            Defines the AcceptedParametersManager to be used.
        kernel_index: integer
            The index of the kernel in the list of kernels in the joint kernel.
        row_index: integer
            The index of the row that should be considered from the accepted_parameters_bds matrix.
        rng: random number generator
            The random number generator to be used.

        Returns
        -------
        np.ndarray
            The perturbed parameter values.

        """
        # Get all current parameter values relevant for this model
        continuous_model_values = accepted_parameters_manager.kernel_parameters_bds.value()[kernel_index]

        # Perturb
        continuous_model_values = np.array(continuous_model_values)
        cov = accepted_parameters_manager.accepted_cov_mats_bds.value()[kernel_index]
        perturbed_continuous_values = rng.multivariate_normal(correctly_ordered_parameters[row_index], cov)

        return perturbed_continuous_values