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.
. I’ve managed to get the MWE up and running (see below).
