def resample_zeroinflation_variables(self):
        """
        There's no way around the fact that we have to look at every
        data point, even the zeros here.
        """
        # TODO: move this to cython?
        T, N, C, D, b = self.T, self.D_emission, self.C, self.D, self.emission_distn.b
        indptr = [0]
        indices = []
        vals = []
        offset = 0
        X = np.hstack((self.gaussian_states, self.inputs))
        for t in range(T):
            # Evaluate probability of data
            y_t = np.zeros(N)
            ns_t = self.data.indices[self.data.indptr[t]:self.data.indptr[t+1]]
            y_t[ns_t] = self.data.data[self.data.indptr[t]:self.data.indptr[t+1]]
            ll = self.emission_distn._elementwise_log_likelihood((X[t], y_t))
            ll = ll.ravel()

            # Evaluate the probability that each emission was "exposed",
            # i.e. p(z_tn = 1 | y_tn, x_tn)
            log_p_exposed = np.log(self.rho) + ll
            log_p_exposed -= np.log(np.exp(log_p_exposed) + (1-self.rho) * (y_t == 0))

            # Sample zero inflation mask
            z_t = np.random.rand(N) < np.exp(log_p_exposed)

            # Construct the sparse matrix
            t_inds = np.where(z_t)[0]
            indices.append(t_inds)
            vals.append(y_t[t_inds])
            offset += t_inds.size
            indptr.append(offset)

        # Construct a sparse matrix
        vals = np.concatenate(vals)
        indices = np.concatenate(indices)
        indptr = np.array(indptr)
        self.masked_data = csr_matrix((vals, indices, indptr), shape=(T, N))
There's no way around the fact that we have to look at every
        data point, even the zeros here.