Gibbs sampling and Bayesian inference I wanted to get a more in-depth understanding of sampling algorithms, and so I thought I could start with the very simple example of a binomial or bernoulli likelihood with a beta prior (since it is solvable in closed form, I thought it might be suitable).
Assume we have trials where have draw $N$ times, and count $k$ successes. $N$ is a deterministic parameter predefined in the experimental setting. We model $k \mid \theta, N \sim \text{Binomial}(\theta, N)$ and model $\theta \sim \text{Beta}(\alpha_0, \beta_0)$. We would like to obtain $\Pr(\theta \mid k, N)$ by sampling, given several observations of $k$, without relying on the closed-form solution of a Beta-posterior.
Looking at MacKay's Information theory, inference and learning algorithms, we would start with an initial $x=(k_0, \theta_0)$, and use the $\text{Binomial}(k\mid \theta_0, N)$ to draw a new $k_1$. This part I can understand, because it is drawing a random sample from a binomial distribution. Then we would use $\Pr(\theta\mid k_1, N)$ to draw $\theta_1$. This part I do not understand anymore, because at this point, we don't have a posterior yet, don't we? So without a posterior distribution $\Pr(\theta\mid k, N)$, how can we sample a new $\theta$?. And furthermore, **shouldn't the algorithm also sample from the observed $k$ values? Or am I missing a step?
 A: If the distribution to simulate is the posterior
$$\pi(\theta|k_0) \propto {N \choose k_0} \theta^{k_0+\alpha_0 - 1}(1-\theta)^{N-k_0+\beta_0 - 1}\propto \theta^{k_0+\alpha_0 - 1}(1-\theta)^{N-k_0+\beta_0 - 1},$$a Gibbs sampler can be constructed by introducing auxiliary variables and making the above density function an integral over these auxiliary variables:
$$\pi(\theta|k_0) = \int \pi(\theta,u|k_0)\text{d}u$$
A standard completion algorithm is the slice sampler, where terms in the density are replaced by integrals of indicator functions:
\begin{align*}\theta^{k_0+\alpha_0-1}&=\int_0^\infty \mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\text{d}u_1\\
(1-\theta)^{N-k_0+\beta_0 - 1}&=\int_0^\infty \mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\text{d}u_2\\\end{align*}
which makes the joint distribution equal to
$$\pi(\theta,u|k_0)\propto\mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\qquad u=(u_1,u_2)$$
The associated Gibbs sampler is then
\begin{align*}U_1|\theta,k_0,U_2 &\sim \mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}} \\
U_2|\theta,k_0,U_2 &\sim  \mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\\
\theta|k_0,U_1,U_2 &\sim\mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}\\\end{align*}
Meaning
\begin{align*}U_1|\theta,k_0,U_2 &\sim \mathcal{U}(0,\theta^{N-k_0+\beta_0 - 1}) \\
U_2|\theta,k_0,U_2 &\sim \mathcal{U}(0,(1-\theta)^{N-k_0+\beta_0 - 1}) \\
\theta|k_0,U_1,U_2 &\sim \mathcal{U}(U_1^{1/{k_0+\alpha_0-1}},1-U_2^{1/{N-k_0+\beta_0 - 1}})\\\end{align*}
since $$\mathbb{I}_{u_1\le\theta^{k_0+\alpha_0-1}}\mathbb{I}_{u_2\le(1-\theta)^{N-k_0+\beta_0 - 1}}=\mathbb{I}_{u_1^{1/{k_0+\alpha_0-1}}\le\theta}\mathbb{I}_{u_2^{1/{N-k_0+\beta_0 - 1}}\le(1-\theta)}$$
An R code implementing this slice sampler is
  sliz=function(T=1e3){
   the0=k0/N #starting value
   thez=rep(0,T)
   for (t in 1:T){
     u1=runif(1)^{1/{k0+a0-1}}*the0
     u2=runif(1)^{1/{N-k0+b0-1}}*(1-the0)
     thez[t]=the0=runif(1,u1,1-u2)}
   return(thez)}

resulting in a neat fit to the intended target

A: The joint distribution 
$$
f(\theta, k) \propto {N \choose k} \theta^{k+\alpha_0 - 1}(1-\theta)^{N-k+\beta_0 - 1} \tag{1}.
$$
That means 
$$
f(\theta \mid k) = \text{Beta}(k + \alpha_0, N-k+\beta_0) \tag{2}.
$$
Gibbs sampling in this case would involve iterations that alternate by drawing from $f(\theta \mid k)$ and $f(k \mid \theta)$ (which is given). This would give you draws from a Markov chain with a stationary distribution equal to the joint distribution above. They would look something like:
$$
(\theta_0, k_0),(\theta_0, k_1)
,(\theta_1, k_1)
,(\theta_1, k_2)
,\ldots
$$ 
You are correct that it's strange to do Gibbs sampling in this situation because you are targeting a joint distribution, and usually in practice you target the posterior. The posterior is completely known (2), and the joint is also known ((1) is a Beta-Binomial). However, this is a common example used for for illustrative purposes.
