issue75.py

import numpy as np

import autoregressive.distributions as d
import autoregressive.models as m

def get_empirical_ar_params(train_datas, params):
    """
    Estimate the parameters of an AR observation model
    by fitting a single AR model to the entire dataset.
    """
    assert isinstance(train_datas, list) and len(train_datas) > 0
    datadimension = train_datas[0].shape[1]
    assert params["nu_0"] > datadimension + 1

    # Initialize the observation parameters
    obs_params = dict(nu_0=params["nu_0"],
                      S_0=params['S_0'],
                      M_0=params['M_0'],
                      K_0=params['K_0'],
                      affine=params['affine'])


    # Fit an AR model to the entire dataset
    obs_distn = d.AutoRegression(**obs_params)
    obs_distn.max_likelihood(train_datas)

    # Use the inferred noise covariance as the prior mean
    # E_{IW}[S] = S_0 / (nu_0 - datadimension - 1)
    obs_params["S_0"] = obs_distn.sigma * (params["nu_0"] - datadimension - 1)
    obs_params["M_0"] = obs_distn.A.copy()

    return obs_params

def genmodel_empirical_estimator(kappa):
    Nmax = 10
    affine = False
    nlags = 3
    D_obs = 3

    prior_ar_params = \
        dict(nu_0=D_obs+2,
             S_0=np.eye(D_obs),
             M_0=np.hstack((np.eye(D_obs), np.zeros((D_obs, D_obs*(nlags-1)+affine)))),
             K_0=np.eye(D_obs*nlags+affine),
             affine=affine)

    prior_ar_params = get_empirical_ar_params([reduced.T], prior_ar_params)

    model = m.ARWeakLimitStickyHDPHMM(
            alpha=4.,
            kappa=kappa,
            gamma=4,
            init_state_distn='uniform',
            obs_distns=[
                d.AutoRegression(**prior_ar_params)
                for state in range(Nmax)],
            )

    model.add_data(reduced.T)
    return model

reduced = np.random.randn(10, 3).T
model = genmodel_empirical_estimator(10**6)