I'm looking for an example of using Gibbs sampling with a bivariate normal, where the correlation parameter is not fixed or known. In other words, what is the conditional distribution of the correlation $\rho$?
There are many examples out there showing Gibbs sampling with the bivariate normal, but they all assume that $\rho$ is known.
Borrowing the LaTeX from Wikipedia, the joint density of a bivariate $(X,Y)$ is given by $$f(x,y) = \frac{1}{2 \pi \sigma_X \sigma_Y \sqrt{1-\rho^2}}\\ \times \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_X)^2}{\sigma_X^2} + \frac{(y-\mu_Y)^2}{\sigma_Y^2} - \frac{2\rho(x-\mu_X)(y-\mu_Y)}{\sigma_X \sigma_Y} \right] \right)$$ Therefore, assuming a flat prior on $\rho\in(-1,1)$, the full conditional posterior of $\rho$ has density proportional to $$\pi(\rho)\propto\frac{1}{\sqrt{1-\rho^2}}\times \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_X)^2}{\sigma_X^2} + \frac{(y-\mu_Y)^2}{\sigma_Y^2}\right]\right)\\\times\exp\left(-\frac{1}{2(1-\rho^2)}\left[-\frac{2\rho(x-\mu_X)(y-\mu_Y)}{\sigma_X \sigma_Y}\right]\right)$$ that is $$\pi(\rho)\propto\frac{1}{\sqrt{1-\rho^2}}\times \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_X)^2}{\sigma_X^2} + \frac{(y-\mu_Y)^2}{\sigma_Y^2}\right]\right)\\ \times \exp\left(\frac{1}{2(1-\rho^2)}\frac{2\rho(x-\mu_X)(y-\mu_Y)}{\sigma_X \sigma_Y}\right)$$ which is thus a density of the form $$g(\rho)\propto (1-\rho^2)^{-1/2}\exp\{-\beta/(1-\rho^2)+\alpha\rho/(1-\rho^2)\}\mathbb{I}_{(-1,1)}(\rho)$$ with $|\alpha|\le\beta$. Since this does not appear to be a standard distribution, one solution is to run Metropolis within Gibbs.
I however answered a very similar question a while ago, using accept reject to simulate exactly this full conditional.