Updating multivariate priors

Questions
5

First of all, many thanks for the prompt replies. I really appreciate. It makes sense to me to use an approximation. Ultimately, I’ll be interested in non-negativity constrained parameters, and thus in something like a Dirichlet distribution approximation (as I read that the Wishart distribution is not recommended within PyMC).

As a first step, I have tried to implement a receding horizon estimator based on a Multivariate Normal-distribution, as suggested by @junpenglao. Based on other examples that I’ve found across the web, I think I got the overall structure right. But I’m not sure exactly how to update the parameters that characterize the distribution of the covariance matrix.

Below is a MWE, in which two lines are tagged with a FIXME note (the same lines appear a bit later again). Any help with correctly specifiying these update statements would be highly appreciated.

import numpy as np
import pymc3 as pm
import theano
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

μ_actual = np.array([1, -2])
Σ_actual = np.array([[0.5, -0.3],
                     [-0.3, 1.]])

def getNewObservations(N=100):
    return np.random.multivariate_normal(μ_actual, Σ_actual, size=N)

def printResults(trace):
    #mean
    print("deviation in μ: ", np.mean(trace["μ"],0) / μ_actual - 1)

    #covariance
    L = np.mean(trace["L"],0)
    print('deviation in Σ: ', L.dot(L.T) / Σ_actual - 1)

initialObservations = getNewObservations(1000)

with pm.Model() as model:
    packed_L = pm.LKJCholeskyCov('packed_L', n=2, eta=2., sd_dist=pm.HalfCauchy.dist(2.5))
    L = pm.Deterministic("L", pm.expand_packed_triangular(2, packed_L))

    μ = pm.Normal('μ', mu=0., sd=10., shape=2, testval=initialObservations.mean(axis=0))
    obs = pm.MvNormal('obs', μ, chol=L, observed=initialObservations)

    stepNUTS = pm.NUTS(vars=[packed_L, μ])
    trace = pm.sample(cores=1, step=[stepNUTS])

printResults(trace) #DEBUG

#### rolling-horizon model
traces = [trace]

#FIXME
eta_cov = theano.shared(???, name="eta_cov")
#FIXME
sd_dist_cov_beta = theano.shared(???, name="sd_dist_cov_beta")

mu_mean = theano.shared(trace['μ'].mean(axis=0), name="mu_mean")
sd_mean = theano.shared(trace['μ'].std(axis=0), name="sd_mean")

new_observations = getNewObservations()
observations_shared = theano.shared(new_observations)

#iteration 0
with pm.Model() as update_model:
    packed_L = pm.LKJCholeskyCov('packed_L', n=2, eta=eta_cov, sd_dist=pm.HalfCauchy.dist(sd_dist_cov_beta))
    L = pm.Deterministic("L", pm.expand_packed_triangular(2, packed_L))
    μ = pm.Normal('μ', mu=mu_mean, sd=sd_mean, shape=2, testval=observations_shared.mean(axis=0))

    obs = pm.MvNormal('obs', μ, chol=L, observed=observations_shared)

    stepNUTS = pm.NUTS(vars=[packed_L, μ])
    trace = pm.sample(cores=1, step=[stepNUTS])

    traces.append(trace)

    printResults(trace) #DEBUG

#iterations 1,2,...
for _ in range(1,3):
    #retrieve new observations
    observations_shared.set_value(getNewObservations())

    #FIXME
    eta_cov.set_value(?)
    sd_dist_cov_beta.set_value(?)

    mu_mean.set_value(trace['μ'].mean(axis=0))
    sd_mean.set_value(trace['μ'].std(axis=0))

    with update_model:
        stepNUTS = pm.NUTS(vars=[packed_L, μ])
        trace = pm.sample(cores=1, step=[stepNUTS])
        traces.append(trace)

    printResults(trace) #DEBUG

If there’s anything else wrong, I’d equally appreciate any comments. Many thanks in advance.