4
\$\begingroup\$

It is usually inefficient to find bugs in a longer code. Unfortunately, after spending too much time in debugging my code, I realize that a good programming habit is important. Please give me some advice about codestyle, design... anything important to write high-quality code.

import numpy as np
from numpy.random import standard_normal, chisquare, multivariate_normal, dirichlet, multinomial
from numpy.linalg import cholesky, inv
from math import sqrt
import math

class GibbsSampler(object):
    """Gibbs sampler for finite Gaussian mixture model

    Given a set of hyperparameters and observations, run Gibbs sampler to estimate the parameters of the model
    """
    def __init__(self, hyp_pi, mu0, kappa0, T0, nu0, y, prior_z):

        """Initialize the Gibbs sampler

        @para hyp_pi: hyperparameter of pi
        @para mu0, kappa0: parameter of Normal-Wishart distribution
        @para T0, nu0: parameter of Normal-Wishart distribution
        @para y: samples draw from Normal distributions
        """
        self.hyp_pi = hyp_pi
        self.mu0 = mu0
        self.kappa0 = kappa0
        self.T0 = T0
        self.nu0 = nu0
        self.y = y
        self.prior_z = prior_z

    def _clusters(self):
        return len(self.hyp_pi)

    def _dim(self):
        """
        Dimension of the data
        """
        return len(mu0)

    def _numbers(self):
        return len(self.y)

    def estimate_clusters(self, iterations, burn_in, lag):
        """
        """
        estimated_clusters = np.zeros(self._numbers(), float)
        for iteration, z, pi, mu, sigma in self.run(iterations, burn_in, lag):
            #  print "Precision = %f" % self._estimate_precision(z)
            count = [z.count(k) for k in range(self._clusters())]
            print count
            if iteration % 10 == 0:
                print "iteration = %d" %iteration

    def _estimate_precision(self, z):
        numbers = self._numbers()
        count = 0
        for i in range(numbers):
            if self.prior_z[i] == z[i]:
                count += 1
        return float(count)/float(numbers)

    def run(self, iterations, burn_in, lag):
        """
        Run the Gibbs sampler
        """
        self._initialize_gibbs_sampler()
        lag_counter = lag
        iteration = 1
        while iteration <= iterations:
            self._iterate_gibbs_sampler()
            if burn_in > 0:
                burn_in -= 1
            else:
                if lag_counter > 0:
                    lag_counter -= 1
                else:
                    lag_counter = lag
                    yield iteration, self.z, self.pi, self.mu, self.sigma
                    iteration += 1

    def _initialize_gibbs_sampler(self):
        """
        This sets the initial values of the parameters.
        """
        clusters = self._clusters()
        numbers = self._numbers()
        self.mu = np.array([self._sampling_Normal_Wishart()[0] for _ in range(clusters)])
        self.pi = dirichlet(hyp_pi, 1)[0]
        self.sigma = np.array([self._sampling_Normal_Wishart()[1] for _ in range(clusters)])
        self.z = np.array([self._multinomial_samples(pi) for _ in range(numbers)])

    def _sampling_Normal_Wishart(self):
        """
        Sampling mu and sigma from Normal-Wishart distribution.

        """
        # Create the matrix A of the Bartlett decomposition from a p-variate Wishart distribution
        d = self._dim()
        chol = np.linalg.cholesky(self.T0)

        if (self.nu0 <= 81+d) and (self.nu0 == round(self.nu0)):
            X = np.dot(chol, np.random.normal(size = (d, self.nu0)))
        else:
            A = np.diag(np.sqrt(np.random.chisquare(self.nu0 - np.arange(0, d), size = d)))
            A[np.tri(d, k=-1, dtype = bool)] = np.random.normal(size = (d*(d-1)/2.))
            X = np.dot(chol, A)
        inv_sigma = np.dot(X, X.T)
        mu = np.random.multivariate_normal(self.mu0, np.linalg.inv(self.kappa0*inv_sigma))

        return mu, np.linalg.inv(inv_sigma)

    def _norm_pdf_multivariate(self, index, cluster):
        """
        Calculate the probability density of multivariable normal distribution
        """
        d = self._dim()
        m = self.y[index] - self.mu[cluster]
        part1 = np.dot(m, np.linalg.inv(self.sigma[cluster]))
        part = np.dot(part1, m.T)
        value = 1.0 / (math.pow(2.0*math.pi, d*0.5) * math.sqrt(np.linalg.det \
            (self.sigma[cluster]))) * math.exp(-(0.5)*part)
        return value

    def _iterate_gibbs_sampler(self):
        """
        Updates the values of the z, pi, mu, sigma.
        """
        clusters = self._clusters()
        # sampling the indicator variables
        pos_z = []
        for i in range(len(self.y)):
            f_xi = np.array([self._norm_pdf_multivariate(i, k) for k in range(clusters)])
            prob_zi = (self.pi * f_xi) / np.dot(self.pi, f_xi)
            pos_zi = self._multinomial_samples(prob_zi)
            pos_z.append(pos_zi)

        # sampling new mixture weights
        count_k = np.array([pos_z.count(k) for k in range(clusters)])
        pos_pi = np.random.dirichlet(count_k + self.pi, 1)[0]

        # sampling parameters for each cluster
        pos_x = []
        for k in range(clusters):
            pos_xk = np.array([self.y[i] for i in range(len(pos_z)) if pos_z[i] == k ])
            pos_x.append(pos_xk)
        # calculate the posterior of multi-normal distribution
        pos_mu = []
        pos_sigma = []
        for k in range(clusters):
            if len(pos_x[k]) == 0: # No observations, no update.
                pos_T0 = self.T0
                pos_mu0 = self.mu0
                pos_kappa0 = self.kappa0
                pos_nu0 = self.nu0
            else:
                # Update the parameters of Normal-Wishart distribution.
                # mean_k is the sample mean in k-th cluster.
                # C is the sample covariance matrix.
                # D is the true covariance matrix.
                mean_k = np.mean(pos_x[k], axis=0)
                C = np.zeros((len(mean_k), len(mean_k)))
                for x_i in pos_x[k]:
                    C += (x_i - mean_k).reshape(len(mean_k), 1) * (x_i - mean_k)
                pos_nu0 = self.nu0 + len(pos_x[k])
                pos_kappa0 = self.kappa0 + len(pos_x[k])
                D = float(self.kappa0  * len(pos_x[k])) / (self.kappa0 + len(pos_x[k])) * \
                    (mean_k - self.mu0).reshape(len(mean_k), 1) * (mean_k - self.mu0)
                pos_T0 = np.linalg.inv(np.linalg.inv(self.T0) + C + D)
                pos_mu0 = (self.kappa0*self.mu0 + len(pos_x[k])*mean_k) / (self.kappa0 + len(pos_x[k]))
                # Update posterior parameters of Normal-Wishart distribution.
                # Then draw the new parameters pos_mu and pos_sigma for each cluster.
            self.mu0 = pos_mu0
            self.kappa0 = pos_kappa0
            self.T0 = pos_T0
            self.nu0 = pos_nu0
            pos_mu_k, pos_sigma_k = self._sampling_Normal_Wishart()
            print pos_mu_k
            pos_mu.append(pos_mu_k)
            pos_sigma.append(pos_sigma_k)

        # After all parameters updated, pass them to the initial values.
        self.z = pos_z
        self.pi = pos_pi
        self.mu = pos_mu
        self.sigma = pos_sigma

    def _multinomial_samples(self, distributions):
        return np.nonzero(multinomial(1, distributions))[0][0]

def multinomial_sample(distributions):
    return np.nonzero(multinomial(1, distributions))[0][0]

def generate_observations(clusters, numbers, hyp_pi = None):
    if hyp_pi == None:
        hyp_pi = [1]*clusters
    pi = dirichlet(hyp_pi, 1)[0]
    mu = []
    sigma = []
    observations = []
    prior_z = []
    for i in range(clusters):
        m, s = sampling_Normal_Wishart(mu0, kappa0, T0, nu0)
        mu.append(m)
        sigma.append(s)
    for i in range(clusters):
        cluster = multinomial_sample(pi)
        obs = multivariate_normal(mu[cluster], sigma[cluster], k_num)
        observations.extend(list(obs))
        prior_z.extend([cluster]*k_num)
    return observations, prior_z

def sampling_Normal_Wishart(mu0, kappa0, T0, nu0):
    """
    Sampling cluster parameters from normal inverse Wishart distribution.
    """
    d = len(mu0)
    chol = np.linalg.cholesky(T0)

    if (nu0 <= 81+d) and (nu0 == round(nu0)):
        X = np.dot(chol, np.random.normal(size = (d, nu0)))
    else:
        A = np.diag(np.sqrt(np.random.chisquare(nu0 - np.arange(0, d), size = d)))
        A[np.tri(d, k=-1, dtype = bool)] = np.random.normal(size = (d*(d-1)/2.))
        X = np.dot(chol, A)
    inv_sigma = np.dot(X, X.T)
    mu = np.random.multivariate_normal(mu0, np.linalg.inv(kappa0*inv_sigma))
    return mu, np.linalg.inv(inv_sigma)

if __name__ == "__main__":
    # Generate the data set.
    # Initialize the parameters for the model.
    # d: dimension of the data.
    # mu0, kappa0, T0, nu0 are the parameters of the Normal-Wishart distribution.
    kappa0 = 4.0
    d = 2
    T0 = np.diag(np.ones(d))
    mu0 = np.zeros(d)
    nu0 = 14.0
    clusters = 6
    k_num = 50
    hyp_pi = [1]*clusters
    pi = dirichlet(hyp_pi, 1)[0]
    y, prior_z = generate_observations(clusters, k_num, hyp_pi = None)
    prior_count = [prior_z.count(k) for k in range(clusters)]
    sampler = GibbsSampler(hyp_pi, mu0, kappa0, T0, nu0, y, prior_z)
    sampler.estimate_clusters(200, 3, 0)
\$\endgroup\$
0

1 Answer 1

Reset to default
8
\$\begingroup\$ 
import numpy as np
from numpy.random import standard_normal, chisquare, multivariate_normal, dirichlet, multinomial
from numpy.linalg import cholesky, inv

You are importing a lot of stuff here, and you aren't even using much of it. I guess not importing these names unless you refer to them a lot.

from math import sqrt
import math

When using numpy, you generally don't want stuff from math. You should use the numpy versions.

class GibbsSampler(object):
    """Gibbs sampler for finite Gaussian mixture model

    Given a set of hyperparameters and observations, run Gibbs sampler to estimate the parameters of the model
    """
    def __init__(self, hyp_pi, mu0, kappa0, T0, nu0, y, prior_z):

        """Initialize the Gibbs sampler

        @para hyp_pi: hyperparameter of pi
        @para mu0, kappa0: parameter of Normal-Wishart distribution
        @para T0, nu0: parameter of Normal-Wishart distribution
        @para y: samples draw from Normal distributions
        """

You don't have a description of prior_z. You also have a lot of parameters here. Four of them describe the Normal-Wishart distribution suggesting perhaps that should be its own class.

        self.hyp_pi = hyp_pi
        self.mu0 = mu0
        self.kappa0 = kappa0
        self.T0 = T0
        self.nu0 = nu0
        self.y = y
        self.prior_z = prior_z

    def _clusters(self):
        return len(self.hyp_pi)

Rather than this, I suggest adding self._clusters = len(self.hyp_pi) to __init__. 

    def _dim(self):
        """
        Dimension of the data
        """
        return len(mu0)

    def _numbers(self):
        return len(self.y)

    def estimate_clusters(self, iterations, burn_in, lag):
        """
        """

Empty docstring, why?

        estimated_clusters = np.zeros(self._numbers(), float)

You create this, but don't seem to do anything with it

        for iteration, z, pi, mu, sigma in self.run(iterations, burn_in, lag):

That's a lot of elements in a tuple. your ignore most of them which raises questions

            #  print "Precision = %f" % self._estimate_precision(z)
            count = [z.count(k) for k in range(self._clusters())]
            print count

Could you use a collections.Counter here?

            if iteration % 10 == 0:
                print "iteration = %d" %iteration

You are mixing output and logic. ITs better to have this class do the calculations, and the calling code print stuff.

    def _estimate_precision(self, z):
        numbers = self._numbers()
        count = 0
        for i in range(numbers):
            if self.prior_z[i] == z[i]:
                count += 1

Use numpy's vector routines for this kind of thing. You should be doing something like:

count = (self.prior_z == z).sum()



        return float(count)/float(numbers)

Add from __future__ import division to the beginning of your script rather then explicitly forcing everything to float.

    def run(self, iterations, burn_in, lag):
        """
        Run the Gibbs sampler
        """
        self._initialize_gibbs_sampler()

My general rule is that all initialization happens in a constructor. The gibbs sample should be initalized in the constructor.

        lag_counter = lag
        iteration = 1
        while iteration <= iterations:
            self._iterate_gibbs_sampler()
            if burn_in > 0:
                burn_in -= 1
            else:
                if lag_counter > 0:
                    lag_counter -= 1
                else:
                    lag_counter = lag
                    yield iteration, self.z, self.pi, self.mu, self.sigma
                    iteration += 1

Your looping logic isn't very clear. I suggest something like:

 for _ in xrange(burn_in):
     self._iterate_gibbs_sampler()

 for _ in xrange(burn_in, iterations, lag):
     for _ in xrange(lag):
         self._iterate_gibbs_sampler()
         yield stuff

I think that more clearly conveys what the code is doing.

    def _initialize_gibbs_sampler(self):
        """
        This sets the initial values of the parameters.
        """
        clusters = self._clusters()
        numbers = self._numbers()
        self.mu = np.array([self._sampling_Normal_Wishart()[0] for _ in range(clusters)])
        self.pi = dirichlet(hyp_pi, 1)[0]
        self.sigma = np.array([self._sampling_Normal_Wishart()[1] for _ in range(clusters)])
        self.z = np.array([self._multinomial_samples(pi) for _ in range(numbers)])

    def _sampling_Normal_Wishart(self):
        """
        Sampling mu and sigma from Normal-Wishart distribution.

        """
        # Create the matrix A of the Bartlett decomposition from a p-variate Wishart distribution
        d = self._dim()
        chol = np.linalg.cholesky(self.T0)

        if (self.nu0 <= 81+d) and (self.nu0 == round(self.nu0)):
            X = np.dot(chol, np.random.normal(size = (d, self.nu0)))
        else:
            A = np.diag(np.sqrt(np.random.chisquare(self.nu0 - np.arange(0, d), size = d)))
            A[np.tri(d, k=-1, dtype = bool)] = np.random.normal(size = (d*(d-1)/2.))
            X = np.dot(chol, A)
        inv_sigma = np.dot(X, X.T)
        mu = np.random.multivariate_normal(self.mu0, np.linalg.inv(self.kappa0*inv_sigma))

        return mu, np.linalg.inv(inv_sigma)

    def _norm_pdf_multivariate(self, index, cluster):
        """
        Calculate the probability density of multivariable normal distribution
        """
        d = self._dim()
        m = self.y[index] - self.mu[cluster]
        part1 = np.dot(m, np.linalg.inv(self.sigma[cluster]))
        part = np.dot(part1, m.T)
        value = 1.0 / (math.pow(2.0*math.pi, d*0.5) * math.sqrt(np.linalg.det \
            (self.sigma[cluster]))) * math.exp(-(0.5)*part)
        return value

    def _iterate_gibbs_sampler(self):
        """
        Updates the values of the z, pi, mu, sigma.
        """
        clusters = self._clusters()
        # sampling the indicator variables
        pos_z = []
        for i in range(len(self.y)):
            f_xi = np.array([self._norm_pdf_multivariate(i, k) for k in range(clusters)])

Do you need to create a python list and then convert to numpy? Can you vectorize that function?

            prob_zi = (self.pi * f_xi) / np.dot(self.pi, f_xi)
            pos_zi = self._multinomial_samples(prob_zi)
            pos_z.append(pos_zi)

When using numpy, you probably shouldn't be appending into python lists. Rewrite to use numpy vectorization.

        # sampling new mixture weights
        count_k = np.array([pos_z.count(k) for k in range(clusters)])
        pos_pi = np.random.dirichlet(count_k + self.pi, 1)[0]

        # sampling parameters for each cluster
        pos_x = []
        for k in range(clusters):
            pos_xk = np.array([self.y[i] for i in range(len(pos_z)) if pos_z[i] == k ])
            pos_x.append(pos_xk)
        # calculate the posterior of multi-normal distribution
        pos_mu = []
        pos_sigma = []
        for k in range(clusters):
            if len(pos_x[k]) == 0: # No observations, no update.
                pos_T0 = self.T0
                pos_mu0 = self.mu0
                pos_kappa0 = self.kappa0
                pos_nu0 = self.nu0
            else:
                # Update the parameters of Normal-Wishart distribution.
                # mean_k is the sample mean in k-th cluster.
                # C is the sample covariance matrix.
                # D is the true covariance matrix.
                mean_k = np.mean(pos_x[k], axis=0)
                C = np.zeros((len(mean_k), len(mean_k)))
                for x_i in pos_x[k]:
                    C += (x_i - mean_k).reshape(len(mean_k), 1) * (x_i - mean_k)
                pos_nu0 = self.nu0 + len(pos_x[k])
                pos_kappa0 = self.kappa0 + len(pos_x[k])
                D = float(self.kappa0  * len(pos_x[k])) / (self.kappa0 + len(pos_x[k])) * \
                    (mean_k - self.mu0).reshape(len(mean_k), 1) * (mean_k - self.mu0)
                pos_T0 = np.linalg.inv(np.linalg.inv(self.T0) + C + D)
                pos_mu0 = (self.kappa0*self.mu0 + len(pos_x[k])*mean_k) / (self.kappa0 + len(pos_x[k]))

Lots and lots of formulas. Suggests that you can perhaps break some of it into functions to greater explicate the high level logic.

                # Update posterior parameters of Normal-Wishart distribution.
                # Then draw the new parameters pos_mu and pos_sigma for each cluster.
            self.mu0 = pos_mu0
            self.kappa0 = pos_kappa0
            self.T0 = pos_T0
            self.nu0 = pos_nu0
            pos_mu_k, pos_sigma_k = self._sampling_Normal_Wishart()
            print pos_mu_k
            pos_mu.append(pos_mu_k)
            pos_sigma.append(pos_sigma_k)

        # After all parameters updated, pass them to the initial values.
        self.z = pos_z
        self.pi = pos_pi
        self.mu = pos_mu
        self.sigma = pos_sigma

    def _multinomial_samples(self, distributions):
        return np.nonzero(multinomial(1, distributions))[0][0]

def multinomial_sample(distributions):
    return np.nonzero(multinomial(1, distributions))[0][0]

def generate_observations(clusters, numbers, hyp_pi = None):
    if hyp_pi == None:
        hyp_pi = [1]*clusters
    pi = dirichlet(hyp_pi, 1)[0]
    mu = []
    sigma = []
    observations = []
    prior_z = []
    for i in range(clusters):
        m, s = sampling_Normal_Wishart(mu0, kappa0, T0, nu0)
        mu.append(m)
        sigma.append(s)
    for i in range(clusters):
        cluster = multinomial_sample(pi)
        obs = multivariate_normal(mu[cluster], sigma[cluster], k_num)
        observations.extend(list(obs))
        prior_z.extend([cluster]*k_num)
    return observations, prior_z

def sampling_Normal_Wishart(mu0, kappa0, T0, nu0):
    """
    Sampling cluster parameters from normal inverse Wishart distribution.
    """
    d = len(mu0)
    chol = np.linalg.cholesky(T0)

    if (nu0 <= 81+d) and (nu0 == round(nu0)):
        X = np.dot(chol, np.random.normal(size = (d, nu0)))
    else:
        A = np.diag(np.sqrt(np.random.chisquare(nu0 - np.arange(0, d), size = d)))
        A[np.tri(d, k=-1, dtype = bool)] = np.random.normal(size = (d*(d-1)/2.))
        X = np.dot(chol, A)
    inv_sigma = np.dot(X, X.T)
    mu = np.random.multivariate_normal(mu0, np.linalg.inv(kappa0*inv_sigma))
    return mu, np.linalg.inv(inv_sigma)

Is this duplicating the method in the class?

if __name__ == "__main__":
    # Generate the data set.
    # Initialize the parameters for the model.
    # d: dimension of the data.
    # mu0, kappa0, T0, nu0 are the parameters of the Normal-Wishart distribution.
    kappa0 = 4.0
    d = 2
    T0 = np.diag(np.ones(d))
    mu0 = np.zeros(d)
    nu0 = 14.0
    clusters = 6
    k_num = 50

Either these should be constants at the beginning of the file and in ALL_CAPS, or you should make them inside a function and thus not accessible as globals.

    hyp_pi = [1]*clusters
    pi = dirichlet(hyp_pi, 1)[0]
    y, prior_z = generate_observations(clusters, k_num, hyp_pi = None)
    prior_count = [prior_z.count(k) for k in range(clusters)]
    sampler = GibbsSampler(hyp_pi, mu0, kappa0, T0, nu0, y, prior_z)
    sampler.estimate_clusters(200, 3, 0)
\$\endgroup\$
1
  • \$\begingroup\$ Thank you very much for helpful advice. You are my teacher! \$\endgroup\$
    – fishiwhj
    Commented Jul 11, 2013 at 6:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.