def METHOD_NAME():
    # Generate correlated data with mean and variance
    np.random.seed(3434)
    x = np.random.normal(size=(1000, 2))
    r = 0.5
    x[:, 1] = r*x[:, 0] + np.sqrt(1-r**2)*x[:, 1]
    x[:, 0] *= 2
    x[:, 1] *= 3
    x[:, 0] += 1
    x[:, 1] -= 2
    # Introduce some missing values
    u = np.random.normal(size=x.shape[0])
    x[u > 1, 0] = np.nan
    u = np.random.normal(size=x.shape[0])
    x[u > 1, 1] = np.nan
    bm = BayesGaussMI(x)
    # Burn-in
    for k in range(500):
        bm.update()
    # Estimate the posterior mean
    mean = 0
    cov = 0
    dmean = 0
    dcov = 0
    for k in range(500):
        bm.update()
        mean += bm.mean
        cov += bm.cov
        dmean += bm.data.mean(0)
        dcov += np.cov(bm.data.T)
    mean /= 500
    cov /= 500
    dmean /= 500
    dcov /= 500
    assert_allclose(mean, np.r_[1, -2], 0.1)
    assert_allclose(dmean, np.r_[1, -2], 0.1)
    assert_allclose(cov, np.asarray([[4, 6*r], [6*r, 9]]), 0.1)
    assert_allclose(dcov, np.asarray([[4, 6*r], [6*r, 9]]), 0.1)