0
$\begingroup$

Is it possible to tell if the parameters can be uniquely estimated in a Bayesian state-space models from the system equations (beyond redundant parameterisations). If so, how?

For example, should it be possible to estimate the parameters in the following model?

\begin{align} x_t &= N(r\times x_{t-1}, \sigma _x)\\ y_t &= Po(\exp(\mu + p\times x_t)) \end{align}

where $x$ is the state variable and $y$ the observations.

This parameterisation is taken from this paper -- and my previous question.


If interested, I am motivated by trying to fit the univariate state-space model with a Gaussian state and Poisson distributed observations below but different chains converge to different values which I think indicates an issue with the parameterisation. Issues remain if I restrict the parameters $p$ and $\mu$ to be positive.

library(rjags)

mod <-
  "model{
    for (i in 2:N){
            x[i] ~ dnorm(r* x[i-1], tau)
            log(my[i]) <- mu + p*x[i]
            y[i] ~ dpois(my[i])
    }
    
    x[1] <- (log(y[1]) - mu)

    r ~ dunif(-1, 1)
    p ~ dunif(-1, 1)
    mu ~ dnorm(0, 0.1)
    sig ~ dunif(0, 10)
    tau <- 1/(sig*sig)
}"


m <- jags.model(textConnection(mod), data=list(y=y, N=NROW(y)), n.chains=2)
update(m, 5000)
s <- coda.samples(m, variable.names=c("r", "p", "mu", "sig"), 1e4)

The data

set.seed(1); 
r = 0.8; sx = 2 
mu = 5; p = 0.5; 
tm = 20
x = y = rep(0,tm); x[1] = rnorm(1, 0, sx) ; y[1] = rpois(1, exp(mu + p*x[1]))
for(i in 2:tm){
  x[i] = rnorm(1, r*x[i-1], sx)
  y[i] = rpois(1, exp(mu + p*x[i]))
}
$\endgroup$
5
  • $\begingroup$ The distribution of the data only depend on $\sigma_x$ and $p$ through the product $\sigma_x p$. So the likelihood has the shape of a ridge and you have to very careful about your choice of prior. $\endgroup$ – Jarle Tufto Jun 26 '20 at 13:23
  • $\begingroup$ Thank you for your comment @JarleTufto, I wonder if you would be willing to expand on them please (although I think they may move too far away fromthe question so I understand if not). Why does the data only depend on $\sigma_x$ and $p$, and how to know the shape of the likelihood? $\endgroup$ – user2957945 Jun 26 '20 at 13:37
  • $\begingroup$ ... and so do i understand correctly that a better choice of prior could see the above model estimated? $\endgroup$ – user2957945 Jun 26 '20 at 13:44
  • $\begingroup$ $x_t$ is a Gaussian AR(1) process with autocovariance function $\operatorname{Cov}(x_t,x_{t-h})=\sigma_x^2 r^{|h|}/(1-r^2)$. The term $p x_t$ is therefore Gaussian with autocovariance function $\operatorname{Cov}(p x_t,p x_{t-h})=p^2 \sigma_x^2 r^{|h|}/(1-r^2)$. Hence, the distribution of the data (and the likelihood) only depend on $\sigma_x$ and $p$ through their product. Along curves at which this product is constant the likelihood is flat, hence it has the shape of a ridge. $\endgroup$ – Jarle Tufto Jun 26 '20 at 13:50
  • $\begingroup$ Thanks @JarleTufto $\endgroup$ – user2957945 Jun 26 '20 at 13:57
2
$\begingroup$

Adding to @Jarle Tufto's comment, the likelihood can be written as

$$ p(y_{1:t} \mid r, \sigma_x, \mu, p) = \int p(y_{1:t} \mid x_{1:t}, \mu, p) \overbrace{p(x_{1:t} \mid r, \sigma_x)}^{{\text{AR(1)}}} dx_{1:t} \tag{1}. $$

Using conditional independence of the observed data, the first factor of the integrand can be written as $$ \prod_{i=1}^t \frac{\exp\left[-\exp(\mu + p\times x_i) \right] \exp(\mu + p\times x_t)^{y_i}}{y_i!}. $$

Assuming $-1 < r < 1$ and that the first $x_1$ is distributed according to the stationary distribution, the second factor of the integrand can be written as

$$ \left(2 \pi \frac{\sigma^2_x}{1 - r^2} \right)^{-1/2}\exp\left[-\frac{x_1^2}{\left(2 \frac{\sigma^2_x}{1 - r^2} \right)} \right] \prod_{i=2}^t \left(2 \pi \sigma^2_x \right)^{-1/2}\exp \left[- \frac{(x_i - rx_{i-1})^2}{2\sigma^2_x} \right]. $$

Now, for $c > 0$, replace $\sigma_x$ with $c \sigma_x$ and $p$ with $p / c$. If you get the same value for (1) (for a given set of observed data $y_{1:t}$), then this is nonidentifiability. If you get two values for (1) that are close, then this can also be a problem, and this might be called weak identifiability. Sorry, I haven't checked the algebra myself :/

Regarding the prior selection, start thinking about (1) as a function in the parameters. It will help convergence if you choose one point on the ridge and favor that with your prior. After you multiply the likelihood by the prior, you will no longer have the nonidentifiability problem. For example, you could reason that it would be nice to have $p = 1$ and have a super informative prior on that, then you wouldn't have this issue after you look at the posterior instead of the likelihood.

$\endgroup$
2
  • 1
    $\begingroup$ Taylor; thank you for working through this with me. I'll try to digest. $\endgroup$ – user2957945 Jun 26 '20 at 15:51
  • 1
    $\begingroup$ Hi: I dealt with a bayesian problem ( issue had to do with the prior ) a looooooooooong time ago and remember this paper being interesting and useful. Possibly it's related to what you are all discussing. I'm not sure to be honest. If not, apologies for noise. stat.columbia.edu/~gelman/research/published/GelmanMeng1991.pdf $\endgroup$ – mlofton Jun 26 '20 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.