def fit(self):
		data = np.squeeze(np.array(self._demonstrations[0])) #np.loadtxt(os.path.join(os.path.dirname(__file__),'example-data.txt'))[:T]
		Nmax = 25

		# and some hyperparameters
		obs_dim = data.shape[1]
		print data.shape
		obs_hypparams = {'mu_0':np.zeros(obs_dim),
                'sigma_0':np.eye(obs_dim),
                'kappa_0':0.25,
                'nu_0':obs_dim+2}
		dur_hypparams = {'alpha_0':2*30,
                 'beta_0':2}

		obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(Nmax)]
		dur_distns = [pyhsmm.distributions.PoissonDuration(**dur_hypparams) for state in range(Nmax)]

		posteriormodel = pyhsmm.models.WeakLimitHDPHSMM(
        	alpha=6.,gamma=6., # these can matter; see concentration-resampling.py
        	init_state_concentration=6., # pretty inconsequential
        	obs_distns=obs_distns,
        	dur_distns=dur_distns)
		
		for d in self._demonstrations:
			posteriormodel.add_data(np.squeeze(np.array(d)),trunc=60) # duration truncation speeds things up when it's possible

		for idx in progprint_xrange(50):
			posteriormodel.resample_model()

		new_segments = []
		for i in range(0, len(self._demonstrations)):
			new_segments.append(self.findTransitions(posteriormodel.states_list[i].stateseq))

		self.segmentation = new_segments
		self.model = posteriormodel

def fit_model(data, N_iters, args={}):
    # Now fit the model with a model using all the data

    default_args = dict(N=N,
                        K_max=K,
                        alpha=alpha,
                        gamma=gamma,
                        alpha_obs=alpha_obs,
                        beta_obs=beta_obs,
                        init_state_concentration=1.0)
    default_args.update(args)
    model = PoissonHDPHMM(**default_args)
    model.add_data(data)

    def _evaluate(model):
        ll = model.log_likelihood()
        return ll, K_used(model)

    def _step(model):
        model.resample_model()
        return _evaluate(model)

    results = [_step(model) for _ in progprint_xrange(N_iters)]
    lls = np.array([r[0] for r in results])
    Ks = np.array([r[1] for r in results])
    return lls, Ks

def train_sticky_hdp_hmm(data_sets, nmax, niter, nsamples):
    """

    :param data_sets: a list of ndarrays.
    :param nmax:
    :param niter:
    :param nsamples:
    :return:
    """
    interval = niter / nsamples
    obs_dim = data_sets[0].shape[-1]
    obs_hypparams = {'mu_0': np.zeros(obs_dim),
                     'sigma_0': np.eye(obs_dim),
                     'kappa_0': 0.25,
                     'nu_0': obs_dim + 2}

    ### HDP-HMM without the sticky bias

    obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(nmax)]
    posteriormodel = pyhsmm.models.WeakLimitStickyHDPHMM(kappa=50., alpha=6., gamma=6., init_state_concentration=1.,
                                                         obs_distns=obs_distns)
    for d in data_sets:
        posteriormodel.add_data(d)

    models = []
    last = None
    for idx in progprint_xrange(niter):
        posteriormodel.resample_model()
        if idx % interval == 0:
            models.append(copy.deepcopy(posteriormodel))
        last = posteriormodel
    models[nsamples - 1] = posteriormodel
    return models

def hsmm_segmenter(factors, width=MEDIAN_WIDTH):
    obs = pre.normalize(factors, axis=1)
    obs_dim = obs.shape[1]
    obs_len = obs.shape[0]
    obs_hypparams = {
        'mu_0': np.mean(obs, axis=0),
        'sigma_0': np.cov(obs, rowvar=0),
        'kappa_0': 0.25,
        'nu_0': obs_dim + 2}

    dur_hypparams = {'alpha_0': 48,
                     'beta_0': 2}

    obs_distns = [distributions.Gaussian(**obs_hypparams) for state in xrange(Nmax)]
    dur_distns = [distributions.PoissonDuration(**dur_hypparams) for state in range(Nmax)]

    posteriormodel = pyhsmm.models.WeakLimitHDPHSMM(
        alpha_a_0=1., alpha_b_0=1./4,
        gamma_a_0=1., gamma_b_0=1./4,
        init_state_concentration=1.,  # pretty inconsequential
        obs_distns=obs_distns,
        dur_distns=dur_distns)

    posteriormodel.add_data(factors, trunc=int(obs_len / 3))

    for idx in textutil.progprint_xrange(150):
        posteriormodel.resample_model()

    labels = posteriormodel.stateseqs
    boundaries, labels = find_boundaries(labels, width)
    best_n_types = len(np.unique(labels))

    if len(boundaries) < best_n_types + 1:
        best_n_types = len(boundaries) - 1

    best_labels = segment_labeling(factors, boundaries, c_method='kmeans', k=best_n_types)
    best_boundaries = np.array(boundaries)

    return best_boundaries, best_labels

##################

Nmax = 10
affine = True
nlags = 2
model = m.ARWeakLimitStickyHDPHMM(
    alpha=4.,gamma=4.,kappa=50.,
    init_state_distn='uniform',
    obs_distns=[
        d.AutoRegression(
            nu_0=2.5,
            S_0=2.5*np.eye(2),
            M_0=np.zeros((2,2*nlags+affine)),
            K_0=10*np.eye(2*nlags+affine),
            affine=affine)
        for state in range(Nmax)],
)

model.add_data(data)

###############
#  inference  #
###############

fig = model.make_figure()
model.plot(fig=fig,draw=False)
for _ in progprint_xrange(300):
    model.resample_model()
    model.plot(fig=fig,update=True)

obs_hypparams = {'mu_0':np.zeros(obs_dim),
                'sigma_0':np.eye(obs_dim),
                'kappa_0':0.25,
                'nu_0':obs_dim+2}

obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(N)]

# Build the HMM model that will represent the fitmodel
fitmodel = pyhsmm.models.HMM(
        alpha=50.,init_state_concentration=50., # these are only used for initialization
        obs_distns=obs_distns)
fitmodel.add_data(data)

print('Gibbs sampling for initialization')

for idx in progprint_xrange(25):
    fitmodel.resample_model()

plt.figure()
fitmodel.plot()
plt.gcf().suptitle('Gibbs-sampled initialization')

print('EM')

likes = fitmodel.EM_fit()

plt.figure()
fitmodel.plot()
plt.gcf().suptitle('EM fit')

plt.figure()

        A=np.eye(D),sigma=0.1*np.eye(D), # TODO remove special case
        nu_0=5,S_0=np.eye(D),M_0=np.zeros((D,P)),K_0=np.eye(P))
    for _ in xrange(Nmax)]

init_dynamics_distns = [
    Gaussian(nu_0=5,sigma_0=np.eye(P),mu_0=np.zeros(P),kappa_0=1.)
    for _ in xrange(Nmax)]

model = WeakLimitStickyHDPHMMSLDS(
    dynamics_distns=dynamics_distns,
    emission_distns=emission_distns,
    init_dynamics_distns=init_dynamics_distns,
    kappa=50.,alpha=5.,gamma=5.,init_state_concentration=1.)


##################
#  run sampling  #
##################

def resample():
    model.resample_model()
    return model.stateseqs[0].copy()


model.add_data(data)
samples = [resample() for _ in progprint_xrange(1000)]

plt.matshow(np.vstack(samples+[np.tile(labels,(10,1))]))

plt.show()

temp = np.concatenate(((0,),truemodel.states_list[0].durations.cumsum()))
changepoints = zip(temp[:-1],temp[1:])
changepoints[-1] = (changepoints[-1][0],T) # because last duration might be censored



#########################
#  posterior inference  #
#########################

Nmax = 25

obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in range(Nmax)]
dur_distns = [pyhsmm.distributions.GeometricDuration(**dur_hypparams) for state in range(Nmax)]

posteriormodel = pyhsmm.models.GeoHSMMPossibleChangepoints(
        alpha=6.,
        init_state_concentration=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)
posteriormodel.add_data(data,changepoints=changepoints)

for idx in progprint_xrange(50):
    posteriormodel.resample_model()

plt.figure()
posteriormodel.plot()

plt.show()

plt.figure()
model.plot()
plt.gcf().suptitle('subHSMM sampled model after {} iterations'.format(ITERATIONS))
plt.savefig('plots/' + testcase + '/subhmm.png')
plt.close()
s = model.states_list[0] 
"""

### HDP-HMM without the sticky bias

obs_distns = [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in xrange(Nmax)]
posteriormodel = pyhsmm.models.WeakLimitHDPHMM(alpha=6.,gamma=6.,init_state_concentration=1.,
                                   obs_distns=obs_distns)
posteriormodel.add_data(data)

for idx in progprint_xrange(ITERATIONS):
    posteriormodel.resample_model()

posteriormodel.plot()
plt.gcf().suptitle('HDP-HMM sampled model after {} iterations'.format(ITERATIONS))
plt.savefig('plots/' + testcase + '/hdp-hmm.png')
plt.close() 

# Some more hypparams

obs_hypparams = {'mu_0':np.zeros(obs_dim),
                'sigma_0':np.eye(obs_dim),
                'kappa_0':0.3,
                'nu_0':obs_dim+5}
dur_hypparams = {'alpha_0':2,
                 'beta_0':2}

    mu_0=np.zeros(obs_dim), sigma_0=np.eye(obs_dim),
    kappa_0=0.25, nu_0=obs_dim+2)

obs_distns = \
    [pyhsmm.distributions.Gaussian(**obs_hypparams) for state in xrange(Nmax)]
model = pyhsmm.models.WeakLimitStickyHDPHMM(
    kappa=50.,alpha=6.,gamma=6.,init_state_concentration=1.,
    obs_distns=obs_distns)
model.add_data(data)


# run inference!
fig = model.make_figure()
model.plot(fig=fig,draw=False)

for _ in progprint_xrange(250):
    model.resample_model()
    model.plot(fig=fig,update=True)



# from moviepy.video.io.bindings import mplfig_to_npimage
# from moviepy.editor import VideoClip

# fig = model.make_figure()
# model.plot(fig=fig,draw=False)

# def make_frame_mpl(t):
#     model.resample_model()
#     model.plot(fig=fig,update=True,draw=False)
#     return mplfig_to_npimage(fig)

    alpha_0=components_per_gmm,
    weights=row)
    for row in init_weights]

##############
#  build MM  #
##############

model = LibraryMM(
        alpha_0=6.,
        components=meta_components)

model.add_data(data)

##################
#  infer things  #
##################

for i in progprint_xrange(50):
    model.resample_model()

plt.figure()
truemodel.plot()
plt.gcf().suptitle('truth')

plt.figure()
model.plot()
plt.gcf().suptitle('inferred')

plt.show()

    weights=row)
    for row in init_weights]

################
#  build HSMM  #
################

dur_distns = [NegativeBinomialIntegerRVariantDuration(np.r_[0.,0,0,1,1,1,1,1],alpha_0=5.,beta_0=5.)
        for state in range(library_size)]

model = LibraryHSMMIntNegBinVariant(
        init_state_concentration=10.,
        alpha=6.,gamma=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)

model.add_data(data)

#####################
#  sample and save  #
#####################

models = [model.resample_and_copy() for itr in progprint_xrange(100)]

with open('models.pickle','w') as outfile:
    cPickle.dump(models,outfile,protocol=-1)

with open('model.pickle','w') as outfile:
    cPickle.dump(models[0],outfile,protocol=-1)

        for state in range(library_size)]

model = LibraryHSMMIntNegBinVariant(
        init_state_concentration=10.,
        alpha=6.,gamma=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)

for data in training_datas:
    model.add_data(data,left_censoring=True)
    # model.add_data_parallel(data,left_censoring=True)

##################
#  infer things  #
##################

train_likes = []
test_likes = []

for i in progprint_xrange(5):
    model.resample_model()
    # model.resample_model_parallel()
    train_likes.append(model.log_likelihood())
    # test_likes.append(model.log_likelihood(test_data,left_censoring=True))

newmodel = model.unfreeze()

for i in progprint_xrange(5):
    newmodel.resample_model()

            np.r_[0,0,0,0,0,1.,1.,1.], # discrete distribution uniform over {6,7,8}
            alpha_0=9,beta_0=1, # average geometric success probability 1/(9+1)
            ) for state in range(Nmax)]

model  = pyhsmm.models.HSMMIntNegBinVariant(
        init_state_concentration=10.,
        alpha=6.,gamma=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)
model.add_data(data,left_censoring=True)

##############
#  resample  #
##############

for itr in progprint_xrange(10):
    model.resample_model()

################
#  viterbi EM  #
################

model.Viterbi_EM_fit()

##########
#  plot  #
##########

plt.figure()
model.plot()

        for state in range(library_size)]

model = LibraryHSMMIntNegBinVariant(
        init_state_concentration=10.,
        alpha=6.,gamma=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)

for data in training_datas:
    model.add_data(data,left_censoring=True)

##################
#  infer things  #
##################

samples1 = [model.resample_and_copy() for i in progprint_xrange(1)]
samples2 = [model.resample_and_copy() for i in progprint_xrange(10)]
samples3 = [model.resample_and_copy() for i in progprint_xrange(100)]
# samples4 = [model.resample_and_copy() for i in progprint_xrange(1000)]

import cPickle
with open('samples1','w') as outfile:
    cPickle.dump(samples1,outfile,protocol=-1)

with open('samples2','w') as outfile:
    cPickle.dump(samples2,outfile,protocol=-1)

with open('samples3','w') as outfile:
    cPickle.dump(samples3,outfile,protocol=-1)

# with open('samples4','w') as outfile:

################

dur_distns = [NegativeBinomialIntegerRVariantDuration(np.r_[0.,0,0,1,1,1,1,1],alpha_0=5.,beta_0=5.)
        for state in range(Nmax)]

model = LibraryHSMMIntNegBinVariant(
        init_state_concentration=10.,
        alpha=6.,gamma=6.,
        obs_distns=obs_distns,
        dur_distns=dur_distns)

for data in datas:
    model.add_data_parallel(data,left_censoring=True)

##################
#  infer things  #
##################

for i in progprint_xrange(25):
    model.resample_model_parallel()

# plt.figure()
# truemodel.plot()
# plt.gcf().suptitle('truth')

# plt.figure()
# model.plot()
# plt.gcf().suptitle('inferred')

# plt.show()

obs_distns = [
        library_models.FrozenMixtureDistribution(
            components=component_library,
            a_0=1.0,b_0=1./5,weights=weights, # for initialization only
        ) for weights in init_weights]

hmm = library_models.LibraryHMM(
        init_state_concentration=10.,
        alpha_a_0=1.0,alpha_b_0=1./10,
        gamma_a_0=1,gamma_b_0=1,
        obs_distns=obs_distns)
hmm.add_data(training_data)

##########################
#  Gather model samples  #
##########################

for itr in progprint_xrange(num_iter):
    hmm.resample_model()

hmms = [hmm.resample_and_copy() for itr in progprint_xrange(1)]

##########
#  Save  #
##########

with open('/scratch/hmm_results.pickle','w') as outfile:
    cPickle.dump(hmms,outfile,protocol=-1)

    gamma_a_0=0.5,
    gamma_b_0=3.0,
    init_state_concentration=6.0,
    obs_distns=[
        d.MNIW(len(dmu) + 2, 5 * np.diag(np.diag(dcov)), np.zeros((len(dmu), 2 * len(dmu))), 10 * np.eye(2 * len(dmu)))
        for state in range(Nmax)
    ],
    # dur_distns=[pyhsmm.basic.distributions.GeometricDuration(100,2000) for state in range(Nmax)],
    # dur_distns=[pyhsmm.basic.distributions.PoissonDuration(3*20,3) for state in range(Nmax)],
    dur_distns=[
        pyhsmm.basic.distributions.NegativeBinomialDuration(10 * 10.0, 1.0 / 10.0, 3 * 10.0, 1 * 10.0)
        for state in range(Nmax)
    ],
)

for i, t in enumerate(tracks):
    model.add_data_parallel(i)

########################
#  try some inference  #
########################
plt.figure()
for itr in progprint_xrange(200):
    if itr % 5 == 0:
        plt.gcf().clf()
        model.plot()
        plt.ion()
        plt.draw()
        plt.ioff()
    model.resample_model_parallel()

Kmax = 10                           # number of latent discrete states
D_latent = 2                        # latent linear dynamics' dimension
D_obs = 2                           # data dimension

Cs = [np.eye(D_obs) for _ in range(Kmax)]                   # Shared emission matrices
sigma_obss = [0.05 * np.eye(D_obs) for _ in range(Kmax)]    # Emission noise covariances

model = DefaultSLDS(
    K=Kmax, D_obs=D_obs, D_latent=D_latent,
    Cs=Cs, sigma_obss=sigma_obss)

model.add_data(data)
model.resample_states()

for _ in progprint_xrange(10):
    model.resample_model()
model.states_list[0]._init_mf_from_gibbs()


####################
#  run mean field  #
####################

vlbs = []
for _ in progprint_xrange(50):
    vlbs.append(model.meanfield_coordinate_descent_step())

plt.figure()
plt.plot(vlbs)
plt.xlabel("Iteration")

    component.add_data(data=observations)
    components.append(component)



# Data structure to collect results
results = []#pd.DataFrame(columns=['Log-Likelihood', 'Accuracy', 'Precision', 'Recall', 'F1'])

# Collect stats into the dataframe
stats = collect_stats(components, real_states, observations, testing, weights, np.ones((1, testing.shape[1])), Model)
results.append(list(stats))


# This is Gibbs sampling
for itr in progprint_xrange(ITER):
    # In each resamle, we are going to fix all the chains except one, then resample
    # that chain given all the other fixed ones. After doing this, a single
    # resample of the factorial model is done

    # Resample the variances
    states = np.matrix(np.zeros((COMP+1, observations.shape[0])))
    states[-1, :] = 1

    for i, component in enumerate(components):
        states[i, :] = component.states_list[0].stateseq

    # Now compute the means
    means = np.matrix(weights)*states
    # Squared summed error
    sse = np.power(observations - means.T, 2).sum(axis=0)