Python scipy.stats 模块，multivariate_normal() 实例源码

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 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 probabilities(self, Y):
        """Probabilities of parameter vectors Y

        Parameters
        ----------
        Y : array-like, shape (num_samples, n_weights)
            2d array of parameter vectors

        Returns
        ----------
        resp : array-like, shape (num_samples)
            the probabilities of the samples under this policy

        """
        if not compatible_version("scipy", ">= 0.14"):
            raise ImportError(
                "SciPy >= 0.14 is required for "
                "'scipy.stats.multivariate_normal'.")
        from scipy.stats import multivariate_normal
        return multivariate_normal(mean=self.mean, cov=self.Sigma).pdf(Y)

def __call__(self, context, explore=True):
        """Evaluates policy for given context.

        Samples weight vector from distribution if explore is true, otherwise
        return the distribution's mean (which depends on the context).

        Parameters
        ----------
        context: array-like, (n_context_dims,)
            context vector

        explore: bool
            if true, weight vector is sampled from distribution. otherwise the
            distribution's mean is returned
        """
        if explore:
            return self.random_state.multivariate_normal(
                self.W.dot(context), self.Sigma, size=[1])[0]
        else:
            return self.W.dot(context)

def gen_blurred_diag_pxy(s):
    X = 1024
    Y = X

    # generate pdf
    from scipy.stats import multivariate_normal
    pxy = np.zeros((X,Y))
    rv = multivariate_normal(cov=s)
    for x in range(X):        
        pxy[x,:] = np.roll(rv.pdf(np.linspace(-X/2,X/2,X+1)[:-1]),int(X/2+x))
    pxy = pxy/np.sum(pxy)

    # plot p(x,y)
    import matplotlib.pyplot as plt
    plt.figure()
    plt.contourf(pxy)
    plt.ion()
    plt.title("p(x,y)")
    plt.show()

    return pxy

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 = 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 loss(self, nodes, model, beta, chunksize=150):
        """
        Forward function used to compute o3 loss
        :param input_labels: of the node present in the batch
        :param model: model containing all the shared data
        :param beta: trade off param
        """
        ret_loss = 0
        for node_index in chunkize_serial(map(lambda x: model.vocab(x).index, nodes), chunksize):
            input = model.node_embedding[node_index]

            batch_loss = np.zeros(len(node_index), dtype=np.float32)
            for com in range(model.k):
                rd = multivariate_normal(model.centroid[com], model.covariance_mat[com])
                # check if can be done as matrix operation
                batch_loss += rd.logpdf(input).astype(np.float32) * model.pi[node_index, com]

            ret_loss = abs(batch_loss.sum())

        return ret_loss * (beta/model.k)

def init_link_labels(self, x0, x1, p):

        pdfs = [multivariate_normal(mean=m, cov=c, allow_singular=True) for m, c in zip(self.means, self.covs)]

        neglog_costs = np.zeros((2, self.means.shape[0]))

        for i in xrange(len(self.means)):
            neglog_costs[0, i] = -np.log(pdfs[i].pdf(x0))
            neglog_costs[1, i] = -np.log(pdfs[i].pdf(x1))

        neglog_link_costs = np.clip(neglog_costs, 0, 9999999999)

        def _cost(k, l):
            return (-np.log(self.weights[k]) * neglog_link_costs[0, k]) + (-np.log(self.weights[l]) * neglog_link_costs[1, l])

        must = must_consistent(self.partition)
        cannot = cannot_consistent(self.partition)

        if p == 1:  # must link
            return must[np.argmin([_cost(k, l) for k, l in must])]
        elif p == -1:  # cannot link
            return cannot[np.argmin([_cost(k, l) for k, l in cannot])]
        else:
            raise ValueError("At this point link is expected to be -1 or 1 only!")

def _int_var_rbf(self, X, hyp, jitter=1e-8):
        """
        Posterior integral variance of the Gaussian Process quadrature.
        X - vector (1, 2*xdim**2+xdim)
        hyp - kernel hyperparameters [s2, el_1, ... el_d]
        """
        # reshape X to SP matrix
        X = np.reshape(X, (self.n, self.d))
        # set kernel hyper-parameters
        s2, el = hyp[0], hyp[1:]
        self.kern.param_array[0] = s2  # variance
        self.kern.param_array[1:] = el  # lengthscale
        K = self.kern.K(X)
        L = np.diag(el ** 2)
        # posterior variance of the integral
        ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X)
        postvar = -ks.dot(solve(K + jitter * np.eye(self.n), ks.T))
        return postvar

def _int_var_rbf_hyp(self, hyp, X, jitter=1e-8):
        """
        Posterior integral variance as a function of hyper-parameters
        :param hyp: RBF kernel hyper-parameters [s2, el_1, ..., el_d]
        :param X: sigma-points
        :param jitter: numerical jitter (for stabilizing computations)
        :return: posterior integral variance
        """
        # reshape X to SP matrix
        X = np.reshape(X, (self.n, self.d))
        # set kernel hyper-parameters
        s2, el = 1, hyp  # sig_var hyper always set to 1
        self.kern.param_array[0] = s2  # variance
        self.kern.param_array[1:] = el  # lengthscale
        K = self.kern.K(X)
        L = np.diag(el ** 2)
        # posterior variance of the integral
        ks = s2 * np.sqrt(det(L + np.eye(self.d))) * multivariate_normal(mean=np.zeros(self.d), cov=L).pdf(X)
        postvar = s2 * np.sqrt(det(2 * inv(L) + np.eye(self.d))) ** -1 - ks.dot(
            solve(K + jitter * np.eye(self.n), ks.T))
        return postvar

def testEntropy(self):
    with self.test_session():
      shift = np.array([[-1, 0, 1], [-1, -2, -3]], dtype=np.float32)
      diag = np.array([[1, 2, 3], [2, 3, 2]], dtype=np.float32)
      actual_mvn_entropy = np.concatenate([
          [stats.multivariate_normal(shift[i], np.diag(diag[i]**2)).entropy()]
          for i in range(len(diag))])
      fake_mvn = self._cls()(
          ds.MultivariateNormalDiag(
              array_ops.zeros_like(shift),
              array_ops.ones_like(diag),
              validate_args=True),
          bs.AffineLinearOperator(
              shift,
              scale=la.LinearOperatorDiag(diag, is_non_singular=True),
              validate_args=True),
          validate_args=True)
      self.assertAllClose(actual_mvn_entropy,
                          fake_mvn.entropy().eval())

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

      log_pdf = mvn.log_prob(x)
      pdf = mvn.prob(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_prob(x)
      pdf = mvn.prob(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_prob(x)
      pdf = mvn.prob(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 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 __init__(self, mean, covariance, randomstate=None):
        self.__mean = npu.immutablecopyof(npu.tondim1(mean))
        self.__covariance = npu.immutablecopyof(npp.checkshapeissquare(npu.tondim2(covariance)))
        self.__randomstate = npu.getrandomstate() if randomstate is None else randomstate
        self.__impl = stats.multivariate_normal(self.__mean, self.__covariance)

def sample(self, size=1):
        return self.__randomstate.multivariate_normal(self.__mean, self.__covariance, size)

def fit_gaussians(x_train_boxcox, y_train):
    """ Fit class-dependent multivariate gaussians on the training set.

    Parameters
    ----------
    x_train_boxcox : np.array [n_samples, n_features_trans]
        Transformed training features.
    y_train : np.array [n_samples]
        Training labels.

    Returns
    -------
    rv_pos : multivariate normal
        multivariate normal for melody class
    rv_neg : multivariate normal
        multivariate normal for non-melody class
    """
    pos_idx = np.where(y_train == 1)[0]
    mu_pos = np.mean(x_train_boxcox[pos_idx, :], axis=0)
    cov_pos = np.cov(x_train_boxcox[pos_idx, :], rowvar=0)

    neg_idx = np.where(y_train == 0)[0]
    mu_neg = np.mean(x_train_boxcox[neg_idx, :], axis=0)
    cov_neg = np.cov(x_train_boxcox[neg_idx, :], rowvar=0)
    rv_pos = multivariate_normal(mean=mu_pos, cov=cov_pos, allow_singular=True)
    rv_neg = multivariate_normal(mean=mu_neg, cov=cov_neg, allow_singular=True)
    return rv_pos, rv_neg

def define_parameters(self):
        params=Parameters()

        params.append(Parameter("trend", multivariate_normal, (self.size,1)))
        params.append(Parameter("sigma2", invgamma, (1,1)))
        params.append(Parameter("lambda2", gamma, (1,1)))
        params.append(Parameter("omega", invgauss, (self.size-self.total_variation_order,1)))

        self.parameters = params

def get_gauss_pdf(sigma):
    n = sigma * 8

    x, y = np.mgrid[0:n, 0:n]
    pos = np.empty(x.shape + (2,))
    pos[:, :, 0] = x
    pos[:, :, 1] = y

    rv = multivariate_normal([n / 2, n / 2], [[sigma ** 2, 0], [0, sigma ** 2]])
    pdf = rv.pdf(pos)

    heatmap = pdf / np.max(pdf)

    return heatmap

def simulation_smoother(self,beta):
        """ Koopman's simulation smoother - simulates from states given
        model parameters and observations

        Parameters
        ----------

        beta : np.array
            Contains untransformed starting values for latent variables

        Returns
        ----------
        - A simulated state evolution
        """         

        T, Z, R, Q, H = self._ss_matrices(beta)

        # Generate e_t+ and n_t+
        rnd_h = np.random.normal(0,np.sqrt(H),self.data.shape[0]+1)
        q_dist = ss.multivariate_normal([0.0, 0.0], Q)
        rnd_q = q_dist.rvs(self.data.shape[0]+1)

        # Generate a_t+ and y_t+
        a_plus = np.zeros((T.shape[0],self.data.shape[0]+1)) 
        a_plus[0,0] = np.mean(self.data[0:5])
        y_plus = np.zeros(self.data.shape[0])

        for t in range(0,self.data.shape[0]+1):
            if t == 0:
                a_plus[:,t] = np.dot(T,a_plus[:,t]) + rnd_q[t,:]
                y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]
            else:
                if t != self.data.shape[0]:
                    a_plus[:,t] = np.dot(T,a_plus[:,t-1]) + rnd_q[t,:]
                    y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]

        alpha_hat, _ = self.smoothed_state(self.data,beta)
        alpha_hat_plus, _ = self.smoothed_state(y_plus,beta)
        alpha_tilde = alpha_hat - alpha_hat_plus + a_plus

        return alpha_tilde

def simulation_smoother(self,beta):
        """ Koopman's simulation smoother - simulates from states given
        model parameters and observations

        Parameters
        ----------

        beta : np.array
            Contains untransformed starting values for latent variables

        Returns
        ----------
        - A simulated state evolution
        """         

        T, Z, R, Q, H = self._ss_matrices(beta)

        # Generate e_t+ and n_t+
        rnd_h = np.random.normal(0,np.sqrt(H),self.data.shape[0]+1)
        q_dist = ss.multivariate_normal([0.0, 0.0], Q)
        rnd_q = q_dist.rvs(self.data.shape[0]+1)

        # Generate a_t+ and y_t+
        a_plus = np.zeros((T.shape[0],self.data.shape[0]+1)) 
        a_plus[0,0] = np.mean(self.data[0:5])
        y_plus = np.zeros(self.data.shape[0])

        for t in range(0,self.data.shape[0]+1):
            if t == 0:
                a_plus[:,t] = np.dot(T,a_plus[:,t]) + rnd_q[t,:]
                y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]
            else:
                if t != self.data.shape[0]:
                    a_plus[:,t] = np.dot(T,a_plus[:,t-1]) + rnd_q[t,:]
                    y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]

        alpha_hat, _ = self.smoothed_state(self.data,beta)
        alpha_hat_plus, _ = self.smoothed_state(y_plus,beta)
        alpha_tilde = alpha_hat - alpha_hat_plus + a_plus

        return alpha_tilde

def simulation_smoother(self,beta):
        """ Koopman's simulation smoother - simulates from states given
        model latent variables and observations

        Parameters
        ----------

        beta : np.array
            Contains untransformed starting values for latent variables

        Returns
        ----------
        - A simulated state evolution
        """         

        T, Z, R, Q, H = self._ss_matrices(beta)

        # Generate e_t+ and n_t+
        rnd_h = np.random.normal(0,np.sqrt(H),self.data.shape[0]+1)
        q_dist = ss.multivariate_normal([0.0], Q)
        rnd_q = q_dist.rvs(self.data.shape[0]+1)

        # Generate a_t+ and y_t+
        a_plus = np.zeros((T.shape[0],self.data.shape[0]+1)) 
        a_plus[0,0] = np.mean(self.data[0:5])
        y_plus = np.zeros(self.data.shape[0])

        for t in range(0,self.data.shape[0]+1):
            if t == 0:
                a_plus[:,t] = np.dot(T,a_plus[:,t]) + rnd_q[t]
                y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]
            else:
                if t != self.data.shape[0]:
                    a_plus[:,t] = np.dot(T,a_plus[:,t-1]) + rnd_q[t]
                    y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]

        alpha_hat,_ = self.smoothed_state(self.data,beta)
        alpha_hat_plus,_ = self.smoothed_state(y_plus,beta)
        alpha_tilde = alpha_hat - alpha_hat_plus + a_plus

        return alpha_tilde

def simulation_smoother(self,beta):
        """ Koopman's simulation smoother - simulates from states given
        model latent variables and observations

        Parameters
        ----------

        beta : np.array
            Contains untransformed starting values for latent variables

        Returns
        ----------
        - A simulated state evolution
        """         

        T, Z, R, Q, H = self._ss_matrices(beta)

        # Generate e_t+ and n_t+
        rnd_h = np.random.normal(0,np.sqrt(H),self.data.shape[0]+1)
        q_dist = ss.multivariate_normal([0.0, 0.0], Q)
        rnd_q = q_dist.rvs(self.data.shape[0]+1)

        # Generate a_t+ and y_t+
        a_plus = np.zeros((T.shape[0],self.data.shape[0]+1)) 
        a_plus[0,0] = np.mean(self.data[0:5])
        y_plus = np.zeros(self.data.shape[0])

        for t in range(0,self.data.shape[0]+1):
            if t == 0:
                a_plus[:,t] = np.dot(T,a_plus[:,t]) + rnd_q[t,:]
                y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]
            else:
                if t != self.data.shape[0]:
                    a_plus[:,t] = np.dot(T,a_plus[:,t-1]) + rnd_q[t,:]
                    y_plus[t] = np.dot(Z,a_plus[:,t]) + rnd_h[t]

        alpha_hat,_ = self.smoothed_state(self.data,beta)
        alpha_hat_plus,_ = self.smoothed_state(y_plus,beta)
        alpha_tilde = alpha_hat - alpha_hat_plus + a_plus

        return alpha_tilde

def fit(self, X, Y):
        """ Fit class-dependent multivariate gaussians on the training set.

        Parameters
        ----------
        x_train_boxcox : np.array [n_samples, n_features_trans]
            Transformed training features.
        y_train : np.array [n_samples]
            Training labels.

        Returns
        -------
        rv_pos : multivariate normal
            multivariate normal for melody class
        rv_neg : multivariate normal
            multivariate normal for non-melody class
        """
        X_boxcox = self._fit_boxcox(X)
        pos_idx = np.where(Y == 1)[0]
        mu_pos = np.mean(X_boxcox[pos_idx, :], axis=0)
        cov_pos = np.cov(X_boxcox[pos_idx, :], rowvar=0)

        neg_idx = np.where(Y == 0)[0]
        mu_neg = np.mean(X_boxcox[neg_idx, :], axis=0)
        cov_neg = np.cov(X_boxcox[neg_idx, :], rowvar=0)
        rv_pos = multivariate_normal(
            mean=mu_pos, cov=cov_pos, allow_singular=True
        )
        rv_neg = multivariate_normal(
            mean=mu_neg, cov=cov_neg, allow_singular=True
        )
        self.rv_pos = rv_pos
        self.rv_neg = rv_neg

def __call__(self, context=None, explore=True):
        """Evaluates policy.

        Samples weight vector from distribution if explore is true, otherwise
        return the distribution's mean.

        Parameters
        ----------
        context : array-like, (n_context_dims,)
            context vector (ignored by this policy, defaults to None)

        explore : bool
            if true, weight vector is sampled from distribution. otherwise the
            distribution's mean is returned

        Returns
        -------
        parameter_vector: array-like, (n_weights,)
            the selected parameters
        """
        # Note: Context is ignored
        if not explore:
            return self.mean
        else:
            # Sample from normal distribution
            return self.random_state.multivariate_normal(
                mean=self.mean, cov=self.Sigma, size=1)[0]

def sample(self, k):
        #samples = self.distribution.rvs(k).reshape(k,p)
        samples = self.rng.multivariate_normal(self.mean, self.cov, k)
        return samples

def pdf(self, x):
        return multivariate_normal(self.mean, self.cov).pdf(x)

def sample(self, k):
        p = len(self.mean)
        if self.df == np.inf:
            chisq = 1.0
        else:
            chisq = self.rng.chisquare(self.df, k) / self.df
            chisq = chisq.reshape(-1,1).repeat(p, axis=1)
        mvn = self.rng.multivariate_normal(np.zeros(p), self.cov, k)
        return self.mean + np.divide(mvn, np.sqrt(chisq))

def _likelihood_statistics(self, img_descriptors):
        """
        :param img_descriptors: X
        :return: 0th order, 1st order, 2nd order statistics
                 as described by equation 20, 21, 22 in reference [1]
        """

        def likelihood_moment(x, posterior_probability, moment):
            x_moment = np.power(np.float32(x), moment) if moment > 0 else np.float32([1])
            return x_moment * posterior_probability

        def zeros(like):
            return np.zeros(like.shape).tolist()

        means, covariances, weights = self.gmm.means, self.gmm.covariances, self.gmm.weights
        normals = [multivariate_normal(mean=means[k], cov=covariances[k]) for k in range(0, len(weights))]
        """ Gaussian Normals """
        gaussian_pdfs = [np.array([g_k.pdf(sample) for g_k in normals]) for sample in img_descriptors]
        """ u(x) for equation 15, page 4 in reference 1 """
        statistics_0_order, statistics_1_order, statistics_2_order = zeros(weights), zeros(weights), zeros(weights)
        for k in range(0, len(weights)):
            for index, sample in enumerate(img_descriptors):
                posterior_probability = FisherVector.posterior_probability(gaussian_pdfs[index], weights)
                statistics_0_order[k] = statistics_0_order[k] + likelihood_moment(sample, posterior_probability[k], 0)
                statistics_1_order[k] = statistics_1_order[k] + likelihood_moment(sample, posterior_probability[k], 1)
                statistics_2_order[k] = statistics_2_order[k] + likelihood_moment(sample, posterior_probability[k], 2)

        return np.array(statistics_0_order), np.array(statistics_1_order), np.array(statistics_2_order)

def gen_easytest(plot=True):

    # set name
    name = "easytest"

    n = 10
    # set generative parameters  
    mu1 = np.array([0,0])
    sig1 = np.eye(2)
    n1 = n
    mu2 = np.array([math.sqrt(75),5])
    sig2 = np.eye(2)
    n2 = n
    mu3 = np.array([0,10])
    sig3 = np.eye(2)
    n3 = n
    param = {'mu1': mu1, 'sig1': sig1, 'n1': n1,
             'mu2': mu2, 'sig2': sig2, 'n2': n2,
             'mu3': mu3, 'sig3': sig3, 'n3': n3}

    # make labels
    labels = np.array([0]*n1+[1]*n2+[2]*n3)

    # make coordinates
    coord = np.concatenate((np.random.multivariate_normal(mu1,sig1,n1),
                            np.random.multivariate_normal(mu2,sig2,n2),
                            np.random.multivariate_normal(mu3,sig3,n3)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_blob(plot=True):

    # set name
    name = "blob"

    # set generative parameters  
    mu1 = np.array([0,0])
    sig1 = np.eye(2)
    n1 = 90
    param = {'mu1': mu1, 'sig1': sig1, 'n1': n1}

    # make labels
    labels = np.array([0]*n1)

    # make coordinates
    coord = np.random.multivariate_normal(mu1,sig1,n1)

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_3sph_evensamp_evenspacing(plot=True):

    # set name
    name = "3sph_evensamp_evenspacing"

    # set generative parameters  
    mu1 = np.array([0,0])
    sig1 = np.eye(2)
    n1 = 30
    mu2 = np.array([math.sqrt(75),5])
    sig2 = np.eye(2)
    n2 = 30
    mu3 = np.array([0,10])
    sig3 = np.eye(2)
    n3 = 30
    param = {'mu1': mu1, 'sig1': sig1, 'n1': n1,
             'mu2': mu2, 'sig2': sig2, 'n2': n2,
             'mu3': mu3, 'sig3': sig3, 'n3': n3}

    # make labels
    labels = np.array([0]*n1+[1]*n2+[2]*n3)

    # make coordinates
    coord = np.concatenate((np.random.multivariate_normal(mu1,sig1,n1),
                            np.random.multivariate_normal(mu2,sig2,n2),
                            np.random.multivariate_normal(mu3,sig3,n3)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_3sph_unevensamp_evenspacing(plot=True):

    # set name
    name = "3sph_unevensamp_evenspacing"

    # set generative parameters  
    mu1 = np.array([0,0])
    sig1 = np.eye(2)
    n1 = 10
    mu2 = np.array([math.sqrt(75),5])
    sig2 = np.eye(2)
    n2 = 30
    mu3 = np.array([0,10])
    sig3 = np.eye(2)
    n3 = 60
    param = {'mu1': mu1, 'sig1': sig1, 'n1': n1,
             'mu2': mu2, 'sig2': sig2, 'n2': n2,
             'mu3': mu3, 'sig3': sig3, 'n3': n3}

    # make labels
    labels = np.array([0]*n1+[1]*n2+[2]*n3)

    # make coordinates
    coord = np.concatenate((np.random.multivariate_normal(mu1,sig1,n1),
                            np.random.multivariate_normal(mu2,sig2,n2),
                            np.random.multivariate_normal(mu3,sig3,n3)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_3sph_evensamp_unevenspacing(plot=True):

    # set name
    name = "3sph_evensamp_unevenspacing"

    # set generative parameters  
    mu1 = np.array([0,2.5])
    sig1 = np.eye(2)
    n1 = 30
    mu2 = np.array([0,-2.5])
    sig2 = np.eye(2)
    n2 = 30
    mu3 = np.array([15,0])    
    sig3 = np.eye(2)
    n3 = 30
    param = {'mu1': mu1, 'sig1': sig1, 'n1': n1,
             'mu2': mu2, 'sig2': sig2, 'n2': n2,
             'mu3': mu3, 'sig3': sig3, 'n3': n3}

    # make labels
    labels = np.array([0]*n1+[1]*n2+[2]*n3)

    # make coordinates
    coord = np.concatenate((np.random.multivariate_normal(mu1,sig1,n1),
                            np.random.multivariate_normal(mu2,sig2,n2),
                            np.random.multivariate_normal(mu3,sig3,n3)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_2cigars(plot=True):

    # set name
    name = "2cigars"

    # set generative parameters   
    mu1 = np.array([0,-4])
    sig1 = np.array([[25,0],[0,1]])
    n1 = 50
    mu2 = np.array([0,4])
    sig2 = np.array([[25,0],[0,1]])
    n2 = 50
    param = {'mu1': mu1, 'sig1': sig1, 'n1': n1,
             'mu2': mu2, 'sig2': sig2, 'n2': n2}

    # make labels
    labels = np.array([0]*n1+[1]*n2)

    # make coordinates
    coord = np.concatenate((np.random.multivariate_normal(mu1,sig1,n1),
                            np.random.multivariate_normal(mu2,sig2,n2)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_halfconcentric(plot=True):

    # set name
    name = "halfconcentric"

    # set generative parameters
    nt = 80 # number of thetas
    nd = 1 # number of samples per theta
    no = nd*nt # number of samples for outer circle
    ni = 20 # number of samples for inner circle
    r = 5 # radius of outer loop
    so = .25 # gaussian noise variance of outer circle
    si = .25 # gaussian noise variance of inner circle
    thetas = -np.linspace(0,math.pi,nt)
    x = [r*math.cos(theta) for theta in thetas]
    y = [r*math.sin(theta) for theta in thetas]
    param = {'nt': nt, 'nd': nd, 'no': no, 'ni': ni, 'r': r, 'so': so, 'si': si}

    # make labels
    labels = np.array([0]*ni+[1]*no)

    # make coordinates
    coord = np.random.multivariate_normal(np.array([0,0]),si*np.eye(2),ni)
    for i in range(len(x)):
        coord = np.concatenate((coord,np.random.multivariate_normal(np.array([x[i],y[i]]),so*np.eye(2),nd)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

def gen_concentric(plot=True):

    # set name
    name = "concentric"

    # set generative parameters
    nt = 80 # number of thetas
    nd = 1 # number of samples per theta
    no = nd*nt # number of samples for outer circle
    ni = 20 # number of samples for inner circle
    r = 8 # radius of outer loop
    so = .25 # gaussian noise variance of outer circle
    si = .25 # gaussian noise variance of inner circle
    thetas = -np.linspace(0,2*math.pi,nt)
    x = [r*math.cos(theta) for theta in thetas]
    y = [r*math.sin(theta) for theta in thetas]
    param = {'nt': nt, 'nd': nd, 'no': no, 'ni': ni, 'r': r, 'so': so, 'si': si}

    # make labels
    labels = np.array([0]*ni+[1]*no)

    # make coordinates
    coord = np.random.multivariate_normal(np.array([0,0]),si*np.eye(2),ni)
    for i in range(len(x)):
        coord = np.concatenate((coord,np.random.multivariate_normal(np.array([x[i],y[i]]),so*np.eye(2),nd)))

    # make dataset
    ds = dataset(coord = coord, labels = labels, gen_param = param, name = name)

    # plot coordinates
    if plot: ds.plot_coord()

    # normalize
    ds.normalize_coord()
    if plot: ds.plot_coord()

    return ds

# CODE BELOW NOT YET ADAPTED TO USE NEW IB DATASET CLASS
# GENERATE ONLY P(X,Y)

def create_mnormal_distr(mu, cov):
    '''
    mu, var are parameters of GMM
    '''
    K = mu.shape[0]
    mnormal_distrs = []
    for k in range(K):
        mnormal_distrs.append(mnormal(mean=mu[k, :], cov=cov[k, ...]))
    return mnormal_distrs

def exp6(num_data=1000):
    if num_data < 2:
        raise ValueError('num_data should be larger than 2. (num_data = {})'.format(num_data))

    center = -5 
    sigma_x = 7 
    sigma_y = 7 

    n1 = num_data 

    # init data 
    d1x = torch.FloatTensor(n1, 1)
    d1y = torch.FloatTensor(n1, 1)
    d1x.normal_(center, sigma_x)
    d1y.normal_(center, sigma_y)

    d1 = torch.cat((d1x, d1y), 1)

    d = d1 

    # label
    label = torch.IntTensor(num_data).zero_()
    label[:] = 0 

    # shuffle
    #shuffle = torch.randperm(d.size()[0])
    #d = torch.index_select(d, 0, shuffle)
    #label = torch.index_select(label, 0, shuffle)

    # pdf
    rv1 = multivariate_normal([ center,  center], [[math.pow(sigma_x, 2), 0.0], [0.0, math.pow(sigma_y, 2)]])

    def pdf(x):
        prob = (float(n1) / float(num_data)) * rv1.pdf(x)
        return prob

    def sumloglikelihood(x):
        return np.sum(np.log((pdf(x) + 1e-10)))

    return d, label, sumloglikelihood

def simulate_a(self, N):
        """
        Assuming that the u's are normal, this method draws a random path 
        for x^N  
        """
        V = self.construct_V(N + 1)
        d = spst.multivariate_normal(np.zeros(N + 1), V)

        return d.rvs()

def _init_params(self, X):
        init = self.init
        n_samples, n_features = X.shape
        n_components = self.n_components

        if (init == 'kmeans'):
            km = Kmeans(n_components)
            clusters, mean, cov = km.cluster(X)
            coef = sp.array([c.shape[0] / n_samples for c in clusters])
            comps = [multivariate_normal(mean[i], cov[i], allow_singular=True)
                     for i in range(n_components)]
        elif (init == 'rand'):
            coef = sp.absolute(sprand.randn(n_components))
            coef = coef / coef.sum()
            means = X[sprand.permutation(n_samples)[0: n_components]]
            clusters = [[] for i in range(n_components)]
            for x in X:
                idx = sp.argmin([spla.norm(x - mean) for mean in means])
                clusters[idx].append(x)

            comps = []
            for k in range(n_components):
                mean = means[k]
                cov = sp.cov(clusters[k], rowvar=0, ddof=0)
                comps.append(multivariate_normal(mean, cov, allow_singular=True))

        self.coef = coef
        self.comps = comps

def from_gaussian(z, R, N):
        """Init new PF from gaussian prior."""
        var = multivariate_normal(z, R)
        s = var.rvs(N)
        w = np.full((N,), 1/N)

        return PF((w, s))

def sample_normal(Q, N):
    """Generate random samples and their weights."""
    var = multivariate_normal(np.zeros(Q.shape[0]), Q)
    s = var.rvs(N)
    return s

def sample(self, n, seed=3):
        with util.NumpySeedContext(seed=seed):
            mvn = stats.multivariate_normal(self.mean, self.cov)
            X = mvn.rvs(size=n)
            if len(X.shape) ==1:
                # This can happen if d=1
                X = X[:, np.newaxis]
            return Data(X)

def load_lip_data_t2(image_id, flip=False, is_test=False):
    fine_size=64
    image_id = image_id[:-1] 
    image_path = './datasets/human/masks/{}.png'.format(image_id)
    img_A = scipy.misc.imread(image_path).astype(np.float)
    rows = img_A.shape[0]
    cols = img_A.shape[1]
    img_A = scipy.misc.imresize(img_A, [fine_size, fine_size])
    img_B = np.zeros((fine_size, fine_size), dtype=np.float64)
    with open('./datasets/human/pose/{}.txt'.format(image_id), 'r') as f:
        lines = f.readlines()
    points = lines[0].split(',')
    for idx, point in enumerate(points):
        if idx % 2 == 0:
            c_ = int(point)
            c_ = min(c_, cols-1)
            c_ = max(c_, 0)
            c_ = int(fine_size * 1.0 * c_ / cols)
        else:
            r_ = int(point)
            r_ = min(r_, rows-1)
            r_ = max(r_, 0)
            r_ = int(fine_size * 1.0 * r_ / rows)
            if c_ + r_ == 0:
                continue
            var = multivariate_normal(mean=[r_, c_], cov=2)
            for i in xrange(fine_size):
                for j in xrange(fine_size):
                    img_B[i, j] += var.pdf([i, j]) * 1.0
    img_A = img_A/127.5 - 1.
    img_BA = np.concatenate((img_B[:,:,np.newaxis], img_A), axis=2)
    # print img_BA.shape
    # img_AB shape: (fine_size, fine_size, input_c_dim + output_c_dim)
    return img_BA

#------------------------------------------------------------

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