@@ -0,0 +1,2022 @@
+%22%22%22 Going to try algorithm 1 from Neal (2000)%0A%22%22%22%0A%0Aimport numpy as np%0Aimport bisect%0A%0A%0Adef pick_discrete(p):%0A    %22%22%22Pick a discrete integer between 0 and len(p) - 1 with probability given by p array.%22%22%22%0A    c = np.cumsum(p)%0A    u = np.random.uniform()%0A    return bisect.bisect_left(c, u)%0A%0A%0Aclass DPMM(object):%0A    %22%22%22Dirichlet Process Mixture Model.%0A%0A    @param conjugate_prior  The conjugate_prior object for whatever model is being inferred.%0A    @param alpha            Concentration parameter.%0A    @param D                Data.%0A    @param theta            Optional initial state for sampler.  Will be drawn from conjugate_prior%0A                            if not specified.%0A    %22%22%22%0A    def __init__(self, conjugate_prior, alpha, D, theta=None):%0A        self.conjugate_prior = conjugate_prior%0A        self.alpha = alpha%0A        self.D = D  # data%0A        self.n = len(self.D)%0A%0A        if theta is None:%0A            # Draw from the prior%0A            theta = conjugate_prior.sample(size=len(D))%0A        self.theta = theta%0A%0A        # Initialize r_i array%0A        self.r_i = self.alpha * np.array(%5Bconjugate_prior.pred(x) for x in D%5D)%0A%0A    def q(self, i):%0A        # compute and return row of q_ij matrix (we only ever need one row at a time).%0A        qs = np.array(%5Bself.conjugate_prior.like1(th_j, x=self.D%5Bi%5D) for th_j in self.theta%5D)%0A        qs%5Bi%5D = self.r_i%5Bi%5D  # cheat by placing r_i at q_ii.%0A        return qs%0A%0A    def update_1_theta(self, i):%0A        x = self.D%5Bi%5D%0A        qs = self.q(i)%0A        p = qs/np.sum(qs)%0A        picked = pick_discrete(p)%0A        if picked == i:  # This corresponds to picking r_i in Neal (2000); i.e. get a new theta%0A            # Neal (2000) H_i is the posterior given a single observation x.%0A            self.theta%5Bi%5D = self.conjugate_prior.post(x).sample()%0A        else:  # reuse an existing theta%0A            self.theta%5Bi%5D = self.theta%5Bpicked%5D%0A%0A    def update_theta(self, n=1):%0A        for j in xrange(n):%0A            for i in xrange(len(self.D)):%0A                self.update_1_theta(i)%0A

Create Dirichlet Process Mixture Model class.