# Gibbs Sampler Initialization

I am trying to implement a Gibbs sampler with three blocks, that's not actually a full conditional Gibbs sampler.

Say I partitition my parameter $$X$$ into three: {$$X_1$$, $$X_2$$, $$X_3$$}.

I inted to sample, in random order, $$X_1|X_2$$ and $$X_2|X_1$$ and then sample $$X_3|X_2,X_1$$.

I am running into trouble initializing the algorithm when the initial state $$X$$ has probability almost zero given the invariant distribution $$\pi$$ (i.e. $$\pi(X(0)) \approx0$$).

If I sample $$X_1|X_2$$ first, the values are so off that I get an almost degenerate $$p(X_2|X_1)$$, so that I cannot sample from it using numerical methods.

$$X_2$$ is actually a matrix and a vector conditional on the inverse of that matrix. The problem I'm having is that I can only sample the zero matrix---which doesn't have an inverse---because of floating point precision.

Am I simply running into a numerical problem, or is there something more fundamentally wrong with my Markov chain? If it were irreducible I should be able to reach a point within my target distribution in finite time.

Is it cheating if I pick $$X(0)$$ such that $$\pi(X(0))\gg0$$? Would I be risking pseudo-convergence?

• I'm not familiar with the idea of "pseudo-convergence" of Gibbs samplers, could you provide a reference? Nevertheless, it is common practice to initialize a chain with as good a guess as we can get for the model parameters, sometimes even via MLE. If you have a good guess, use it to initialize. But frankly, just on the basis of my experience and what you've described (particularly the line "If I sample $X1|X2$ first..."), it sounds like there is a bug in the code (but I haven't seen it!). Commented May 25, 2023 at 14:45
• @JohnMadden, the term pseudo-convergence I took from the Handbook of Markov Chain Monte Carlo by Chapter 1 Section 11.2 by Geyer. He actually advocates against the practice of starting with diffuse points to assess convergence, and maybe I should listen. But I wanted to compare the convergence of two different parametrizations of the same model. There is a bug, it's that sampling $X_2$ is not numerically feasible for some values $X_1$. I'm supposed to get a positive definite matrix but instead get the zero matrix. I'm using a constrained Wishart distribution. Commented May 25, 2023 at 15:00
• Whenever I suspect an issue is due to numerical error rather than a mistake I made in my code, I like to run the code on the smallest possible problem I can, where I don't expect numerical noise to be relevant, to double check that assumption. Have you verified that it is indeed numerical noise in this fashion? Commented May 25, 2023 at 15:03
• It does work if I initialize it somewhere near the OLS estimates. I suspect the numerical problem is because I'm using an inverse CDF method to sample a constrained chi square variable ( chi square distribution truncated at some value before the CDF actually has a strictly positive value). Commented May 25, 2023 at 17:04
• to confirm this you could try plugging in a less efficient but more reliable method for this step, such as using rejection sampling on top of a standard chi square sampler. Commented May 25, 2023 at 17:06