Autoregressive Poisson model: Unable to recover latent process parameters

Questions
1

Hi,

I am setting up an autoregressive Poisson model and I’m using synthetic data to ensure that I am recovering the parameters correctly. Unfortunately, that is not happening.

I am actually building a hierarchical model but let’s simplify for now with just one series. Here is the data generating process:

N = 52*6 # number of points per series

x = np.zeros((N))

# true k:
true_k = 0.6
# true process mean:
true_mu=5
# true variance of the innovation
true_tau = 8


x[0]=true_mu

for i in range(1, N):
    x[i] = true_k * (x[i-1] - true_mu) + stats.norm.rvs(0, np.sqrt(1/true_tau)) + true_mu

y = np.random.poisson(np.exp(x))
_, ax = plt.subplots(1, 1, figsize=(8,5))

ax.plot(y);

I use a non-zero mean autoregressive process:

import theano as T
import theano.tensor as tt
import pymc3 as pm
from pymc3.distributions.distribution import generate_samples, draw_values

class NonzeroMeanAR1(pm.distributions.continuous.PositiveContinuous):
    """
    Autoregressive process with 1 lag.
    Parameters
    ----------
    k : tensor
       effect of lagged value on current value
    tau_e : tensor
       precision for innovations
    """

    def __init__(self, k, tau_e, mu, *args, **kwargs):
        super(NonzeroMeanAR1, self).__init__(*args, **kwargs)
        self.k = k = tt.as_tensor_variable(k)
        self.tau_e = tau_e = tt.as_tensor_variable(tau_e)
        self.tau = tau_e * (1 - k ** 2)
        self.mu = tt.as_tensor_variable(mu)
        self.mean = self.mu

    def logp(self, x):
        k = self.k
        tau_e = self.tau_e
        mu = self.mu

        x_im1 = x[:-1]
        x_i = x[1:]
        boundary = pm.Normal.dist(mu=mu, tau=tau_e).logp(x[0])

        innov_like = pm.Normal.dist(mu=k * (x_im1 - mu) + mu, tau=tau_e).logp(x_i)
        return boundary + tt.sum(innov_like)

First, you can find the model that assumes that the AR process is observed:

with pm.Model() as par_model:
    k = pm.Normal('k', 0, 1)
    tau = pm.HalfNormal('tau', sd=1)
    mu = pm.Normal('mu', mu=5, sd=1)
    lat = NonzeroMeanAR1('lat', k=k, tau_e=tau, mu=mu, observed=x)
    #y_pred = pm.Poisson('y_pred', mu=tt.exp(lat), observed=y)
    trace = pm.sample(
        2000,
        tune=2000,
        target_accept=0.9,
        init="advi+adapt_diag",
    )

Next, you can find the model that assumes that the AR process is a latent process:

with pm.Model() as par_model:
    k = pm.Normal('k', 0, 1)
    tau = pm.HalfNormal('tau', sd=1)
    mu = pm.Normal('mu', mu=5, sd=1)
    lat = NonzeroMeanAR1('lat', k=k, tau_e=tau, mu=mu, shape=len(y))
    y_pred = pm.Poisson('y_pred', mu=tt.exp(lat), observed=y)
    trace = pm.sample(
        2000,
        tune=2000,
        target_accept=0.9,
        init="advi+adapt_diag",
    )

I can recover the parameters for a purely AR process but when performing the AR-Poisson (AR is a latent process) I am unable to get the AR process parameters back. Here you have the trace plots:

image
image1173×585 164 KB

It is mixing ok but the estimations are completely off. Already tried fewer data and more informative priors and the opposite. What am I missing?

Thank you!

2

What happens if you sample from your model using the recovered parameters? I suspect this model isn’t identified because you only have a single data point for each time point in the Poisson mean.

Another thing you can do is draw 1000 trajectories from your data generating process and fit the combined Poisson-AR1 model to all of that data. That gives you 1000 samples to fit the Poisson mean (you hopefully do not need that many). This (at least) will give you confidence that the inference is working to recover the parameters.

If the latter “works” then you need to rethink your priors if you really want to “recover” your parameters. That said, if your recovered parameters are producing data that “looks like” your observations, is that good enough?

Dan

1 Like
3

Thanks for the tips @mathDR. I made it work, i.e., I can recover the parameters. But your questions raised some concerns.

What happens if you sample from your model using the recovered parameters? I suspect this model isn’t identified because you only have a single data point for each time point in the Poisson mean.

I understand that I only have a data point for each time point to supply as the parameter to the Poisson likelihood, but what could be another approach? I will link this to your second point as I understood that it could be a suggestion.

Another thing you can do is draw 1000 trajectories from your data generating process and fit the combined Poisson-AR1 model to all of that data. That gives you 1000 samples to fit the Poisson mean (you hopefully do not need that many). This (at least) will give you confidence that the inference is working to recover the parameters.

I am just generating the data to make sure that the model is well defined. I will not be able to do this when applying it to real data. So, I think that this suggestion would only work when I control the generative process of the data, right?

If the latter “works” then you need to rethink your priors if you really want to “recover” your parameters. That said, if your recovered parameters are producing data that “looks like” your observations, is that good enough?

As I said previously, the goal is not to recover the parameters or produce data that looks like the observations. The goal is to make sure that the model is well defined and apply it to real data. And the presented model is just a small part of the overall model but somewhat an important one.

With all of this said, is there another possible approach to make it identified?

4

My comments were geared to running your model on “real” data. Is the purpose of the model to do forecasting? Latent variable inference? Etc.

This is important because there may be several sets of latent parameter values that lead to “accurate forecasts”, but if you want to say something like “the true process mean is 5” because that is what the inference told you, then that is something different.

May I ask: how did you change the model to recover the parameters for the Poisson variant?

5

Right! The purpose of the model is to do forecasting. Based on that, do you have any suggestion to change my approach?

May I ask: how did you change the model to recover the parameters for the Poisson variant?

Of course.

I simplified things by adding an intercept term (in this case log(alpha) approx. 100).

with pm.Model() as par_model:
    k = pm.Normal('k', 0, 1)
    tau = pm.HalfNormal('tau', sd=1)
    alpha = pm.Normal('alpha', mu=5, sd=1)
    lat = pm.AR1('lat', k=k, tau_e=tau, shape=len(y))
    y_pred = pm.Poisson('y_pred', mu=tt.exp(alpha + lat), observed=y)
    trace = pm.sample(
        2000,
        tune=2000,
        target_accept=0.9,
        init="advi+adapt_diag",
    )

The mixing for alpha is not that good though.

AR_intercept
AR_intercept1178×602 175 KB