I am attempting to approximate the posterior $\tilde{\pi_{G}}(z|\theta,Y)$ which is the Gaussian approximation to the full conditional of $z$, and in order to do this I need to find the mode $z^{*} \equiv z^{*}(\theta, Y).$

Firstly, giving the definition for parameters/variables, we have $$ Y_{t} | \log(\lambda_{t}) = x_{t} \sim^{i.i.d} Poi (\log(\lambda _{t})),$$ where $x_{t} = \log(\lambda_{t})$ can represent the pattern using the latent autoregressive process, $$ x_{t} = \beta_{0} + \beta_{1} x_{t-1} + \eta_{t},$$ where $\eta_{t} \sim^{i.i.d} N(0, \sigma_{\eta}^{2}).$ Then taking $z_t = x_t - \mu$, we can obtain a nicer form of, $z_{t} = \beta_{1}z_{t-1} + \eta_{t},$ which is the $z$ used above. This is $z \sim MVN (0, \Sigma_{z}),$ and we see that $z \sim AR(1)$. To simplify notation we have above, $\theta = (\beta_0,\beta_1,\sigma^{2}_{\eta})$ being the parameters, and $Y$ is the given observed data.

In order to calculate the mode, we can apply the newton method, $$z^{n} = z^{n-1} + (\nabla^{2}\log(\pi(z^{n-1}|\theta, Y)))^{-1} \nabla \log(\pi(z^{n-1}|\theta,Y)).$$

I have calculate these values as, $$ \nabla_{z}\log(\pi(z|\theta,Y)) = \nabla_{z}\log(\pi(z|\theta)) + \nabla_{z}\log(\pi(Y|z)), \quad \text{and}$$ $$ \nabla_{z}^{2}\log(\pi(z|\theta,Y)) = \nabla_{z}^{2}\log(\pi(z|\theta)) + \nabla_{z}^{2}\log(\pi(Y|z)),$$ which are simply the gradient and hessian of prior and likelihood function (as we can eliminate any constant terms not involving z, also taking log as it simplifies calculations without changing a stationary point). We can find these through simple calculations as, $$z^{n} = z^{n-1} + (-\Sigma^{-1} - \exp{(\mu)}diag(\exp(z))^{-1} (-\Sigma^{-1}z - \exp{(z)}\exp{(\mu)}+ Y).$$

I have created a code with some synthetic data in order to compute this, however, it seems that with more iterations the magnitude of mode increases, clearly incorrect. Can anyone see what is going wrong? Or any advice?

My code for this is:

# Synthetic data generation 
beta0 <- 0
beta1 <- 0.5
varwn <- 1.5
theta<- c(beta0,beta1,varwn)
n = 100
x <- NULL
x[1] <- rnorm(n=1, mean = (beta0/(1-beta1)), sd = sqrt(varwn/(1-beta1^2)))
for (i in 2:n) {
  e <- rnorm(n=1, mean = 0, sd = sqrt(varwn))
  x[i] <- beta0 +beta1*x[i-1] + e
lambda <- exp(x)
Y=rpois(n,lambda = lambda)
data <- Y

# Initialisation of parameters
n     <- length(data)
beta0 <- theta[1]
beta1 <- theta[2]
sigmasqe <- exp(theta[3])
mu       <- (beta0)/(1-beta1)
sigmasq  <- (sigmasqe)/(1-beta1^2)
z<- Matrix(rep(0,n))
znew <- z

# Precision matrix Q = Sigma^-1, (inverse of covariance matrix)
invsigMat            <- Matrix(diag(1 + beta1^2,n))
diag(invsigMat[-1,]) <- -beta1
diag(invsigMat[,-1]) <- -beta1
invsigMat[1,1]       <- 1; invsigMat[n,n] <- 1;
inverseSigMat         = Matrix((1/(sigmasqe))*invsigMat)

#Newton-Raphson, with M iterations
M <- 20
for (i in 1:M){
  z = znew
  gradLikelihood  <- -exp(mu)*exp(z) + data
  gradPrior       <- -Matrix(inverseSigMat%*%z)
  totgrad         <- gradPrior + gradLikelihood

  D <- Matrix(diag(as.vector(exp(z))))
  hessianLikelihood  <- -exp(mu)*D
  hessianPrior       <- -inverseSigMat
  tothessian         <- hessianPrior + hessianLikelihood

  NRstep   <- solve(tothessian, totgrad) # x <- A^{-1}b
  znew = z + NRstep

Your Answer

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

Browse other questions tagged or ask your own question.