I have been trying to implement a Bayesian inference procedure from scratch for a specific problem, but I have implemented the procedure, and it doesn't seem to work.
Since, I can't just post the code online and ask community to debug my code, I was wondering if someone could provide with a broader checklist when going about coding up a Bayesian inference procedure. (regardless of language)
EDIT: Specifics of the problem
I am trying to implement the procedure described in Section 5 of this paper on MATLAB . Briefly put, the procedure I've implemented is -
- I have 3 zero mean variables (i.e., $D = 3$ time series) for $500$ timepoints. I'm using initial $N = 350$ data points as training sample.
- The covariance function I'm using is a squared exponential kernel with 1 hyperparameter - characteristic length scale $l$. I'm assuming it to be the same for all 3 timeseries.
- I'm keeping degrees of freedom constant, $\nu = D + 1$.
- $L$, the lower Cholesky decomposition of the scale matrix $V$ is computed as the $D \times D$ covariance matrix of the $N \times D$ training dataset.
The sampling procedure essentially involves 2 steps (using Gibbs sampling)
5.1 Sample $u$ ($N \times D \times \nu$) dimensional vector, assuming Gaussian process prior (as defined in equation 19 of the paper). I've assumed a Gaussian likelihood function (as defined in equation 24). For this I'm using Elliptical Slice Sampling
5.2 Sample GP hyperparameter $l$, using a lognormal prior (assumption, $mean=1.5$, $var = 1$). I've used slice sampling for this with posterior as product of GP prior(eq. 19) and lognormal density.
I let this Gibbs sampler run for $10000$ iterations ($5000$ burn-in). But convergence plot of $u$ doesn't seem to converge.
I also tried this with smaller $N$ (~ $50$) and increased no. of iterations but didn't work.