Implementing MultivariateNormalFullCovariance for parameter estimation

I have had success using the following model to parameterize a 2-coefficient function with linearly-combined input/output data from that system. I would like to create a 2-dimensional confidence region (ellipse) instead of having two independent normal distributions (requiring conversion to a covariance-inclusive model). I have tried using similar syntax, increasing the 'scalar' array to a k x k 'covariance' matrix, but I haven't had any luck. Does anybody have any ideas about how I can expand my model to parameterize a 2D normal distribution instead of two independent normals? Would be greatly appreciated. Thanks.

def mutable_model_nomulti(self):
"Probabilistic model for SISO data"
negloglik = lambda y, rv_y: -rv_y.log_prob(y[:])
model = tf.keras.Sequential([
tf.keras.layers.LSTM(100, activation='tanh', input_shape=(None, 1)),
tfp.layers.DenseVariational(4*1,Model.posterior_mean_field,Model.prior_trainable,activation='linear',kl_weight=1/100),
tfp.layers.DistributionLambda(lambda t: tfd.Normal(loc = t[..., :2],
scale = (1e-3 + tf.math.softplus(0.1 * t[...,2:])),)),])
model.compile(optimizer='adam', loss=negloglik,metrics=[Model.coeff_determination])
return model

Start from the example model in the issue body, which uses tf.keras.layers.LSTM, tfp.layers.DenseVariational, and tfd.Normal. No repository file, test, expected API, or acceptance criteria is identified; clarify the desired covariance-inclusive behavior and how a two-dimensional confidence region should be validated.