code: def resample_zeroinflation_variables(self):
          
          # 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))
docs: There's no way around the fact that we have to look at every
          data point, even the zeros here.