Updating multivariate priors

Questions
7

Dear Chris,

Many thanks for your prompt reply. Your suggestion to directly approximate L makes a lot of sense to me. I have tried to implement it in my MWE as follows:

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]

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

packed_L_mean = theano.shared(trace['packed_L'].mean(axis=0), name="packed_L_mean")

packed_L_cov = theano.shared(np.cov(trace['packed_L'], rowvar=False) + 1e-3 * np.diag(np.ones(packed_L_mean.get_value().shape[0])), name="packed_L_cov") # a little regularization
packed_L_cholL = theano.shared(np.linalg.cholesky(packed_L_cov.get_value()), name="packed_L_cholL")

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

#iteration 0
with pm.Model() as update_model:
    packed_L = pm.MvNormal('packed_L', mu=packed_L_mean, chol=packed_L_cholL)
    L = pm.Deterministic("L", pm.expand_packed_triangular(2, packed_L))

    μ = pm.Normal('μ', mu=mu_mean, sd=mu_std, 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

#iteration 1+
for _ in range(1,3):
    #retrieve new observations

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

    packed_L_mean.set_value(trace['packed_L'].mean(axis=0))
    packed_L_cov.set_value(np.cov(trace['packed_L'], rowvar=False) + 1e-3 * np.diag(np.ones(packed_L_mean.get_value().shape[0])))
    packed_L_cholL.set_value(np.linalg.cholesky(packed_L_cov.get_value()))

    observations_shared.set_value(getNewObservations())

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

    printResults(trace) #DEBUG

The first line in the update_model (line 57)

packed_L = pm.MvNormal('packed_L', mu=packed_L_mean, chol=packed_L_cholL)

yields the following error:

TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Do you have any idea where this error comes from, and how I could solve it? I am still a bit confused about when to work with L and L_packed, but I don’t think this is the issue above.

Many thanks in advance!