@@ -0,0 +1,2165 @@
+from distributions.dbg.models import bb%0Afrom microscopes.models.mixture.dp import DirichletProcess%0Afrom microscopes.common.util import almost_eq%0Afrom scipy.misc import logsumexp%0A%0Aimport numpy as np%0Aimport numpy.ma as ma%0Aimport itertools as it%0A%0Afrom common import KL_discrete%0Afrom nose.plugins.attrib import attr%0A%0Adef test_sample_post_pred_no_missing_data():%0A    D = 5%0A    N = 1000%0A    alpha = 2.0%0A%0A    mm = DirichletProcess(N, %7B'alpha':alpha%7D, %5Bbb%5D*D, %5B%7B'alpha':1.,'beta':1.%7D%5D*D)%0A%0A    Y_clustered, _ = mm.sample(N)%0A    Y = np.hstack(Y_clustered)%0A    assert Y.shape%5B0%5D == N%0A    mm.fill(Y_clustered)%0A%0A    Y_samples = %5Bmm.sample_post_pred(y_new=None) for _ in xrange(10000)%5D%0A    Y_samples = np.hstack(Y_samples)%0A%0A    empty_groups = list(mm.empty_groups())%0A    if len(empty_groups):%0A        for egid in empty_groups%5B1:%5D:%0A            mm.delete_group(egid)%0A    else:%0A        mm.create_group()%0A    assert len(mm.empty_groups()) == 1%0A%0A    def score_post_pred(y):%0A        %22%22%22compute log p(y %7C C, Y)%22%22%22%0A        def score_for_group(gid):%0A            ck = mm.nentities_in_group(gid)%0A            ctotal = mm.nentities()%0A            top = ck if ck else alpha%0A            score_assign = np.log(top/(ctotal + alpha))%0A            score_value = sum(g.score_value(mm.get_feature_hp_shared(fi), yi) for fi, (g, yi) in enumerate(zip(mm.get_suff_stats_for_group(gid), y)))%0A            return score_assign + score_value%0A        return logsumexp(np.array(%5Bscore_for_group(gid) for gid in mm.groups()%5D))%0A%0A    scores = np.array(list(map(score_post_pred, it.product(%5BFalse, True%5D, repeat=D))))%0A    scores = np.exp(scores)%0A    assert almost_eq(scores.sum(), 1.0)%0A%0A    # lazy man%0A    idmap = %7B y : i for i, y in enumerate(it.product(%5BFalse, True%5D, repeat=D)) %7D%0A%0A    smoothing = 1e-5%0A    sample_hist = np.zeros(len(idmap), dtype=np.int)%0A    for y in Y_samples:%0A        sample_hist%5Bidmap%5Btuple(y)%5D%5D += 1.%0A    #print 'hist', sample_hist%0A%0A    sample_hist = np.array(sample_hist, dtype=np.float) + smoothing%0A    sample_hist /= sample_hist.sum()%0A%0A    #print 'actual', scores%0A    #print 'emp', sample_hist%0A    kldiv = KL_discrete(scores, sample_hist)%0A    print 'KL:', kldiv%0A%0A    assert kldiv %3C= 0.05%0A%0A

add post pred sampling test for DP