2-step sampling from a conditional density

The setting is as follows:

We are given two random variables $X : \Omega \to \mathbb{R}$ and $\Theta : \Omega \to T$ for some 'parameter space' $T \subset \mathbb{R}$, and

1) we know the density of $\Theta$ and $f_\Theta$, and

2) for every $\theta \in T$, we know the conditional density $f_{X|\Theta} = f_{X,\Theta}/f_{\Theta}$.

(NB: Let us assume that $f_\Theta(\theta) > 0$ for all $\theta$). The intuition behind it is that the density of $X$ has a parameter (such as the normal distribution with parameter $\theta = (\mu, \sigma^2)$) and this parameter is yet again a value of a random variable.

Now I want to sample from the two random variables $(X, \Theta)$ and $X$. I could theoretically compute the densities $f_{X|\Theta}$ [and $f_X$] from the given inputs and sample from these functions. This however, can be analytically 'difficult' up to intractable. So instead people do it by repeating the following steps:

1) Sample a single value $\theta$ from the known density $f_\Theta$.

2) Take the concrete value and insert it into $f_{X|\Theta}$, i.e. compute the density $f_{X|\Theta}(\cdot,\theta)$.

3) Sample a single value $x$ from this density.

4) Output $(x, \theta)$.

My questions are rather simple:

Why does the method above yield a sampling that 'follows' the density $f_{X,\Theta}$?

How to show that a sampling process 'follows' a distribution?

NOTES:

I tried to define a sampling process as a map $a : \mathbb{N} \to \mathbb{R}$ [i.e. a sequence]. The parameter essentially tells us how often we have called this function already. In reality, this function also depends on the random seed and so on, but let us ignore this for now. I suspect that 'a sampling process follows a density $f$' means that continuously drawing histograms with increasing amount of samples makes the histogram get 'closer and closer' to the density.

So let me try to define what I mean: Let us assume that we are given an infinite sequence of partitionings of $\mathbb{R}$, $\mathcal{A}_n$ with $\sup_{A \in \mathcal{A}_n} \text{diam} A \rightarrow 0$ as $n \rightarrow \infty$ (meaning that the histograms become finer and finer. For every $x \in \mathbb{R}$ we have an assigned 'box', $x$ belongs to in the partition $\mathcal{A}_n$, lets call this $\text{box}_n(x)$

So one possible definition would be

$a$ follows $f$ iff. for every $x \in \mathbb{R}$ we have $$\frac{|\{t \in \{1,...,n\} : a(t) \in \text{box}_n(x)\}}{n} \rightarrow f(x)$$

I am aware of the fact that $f_{X,\Theta} = f_{X|\Theta} f_\Theta$. So the question should actually be: Why does the sampling algorithm above 'follow' $f_{X|\Theta} f_\Theta$? Tried though I have, I was unable to prove that if a sampling process $a$ follows $f_\Theta$ and a function $b(n,\theta)$ is such that for every $\theta \in T$, $b(\cdot, \theta)$ follows $f_{X|\Theta}(\cdot,\theta)$, then $c(n) := (b(n, a(n)), a(n))$ follows $f_{X,\Theta}$. The difficulty here is that I do not see the relation between the sets $$\{t : c(t) \in \text{box}_n(x) \times \text{box}_n(y)\}$$ and $$\{(t_1, t_2) : b(t_1, a(t_1)) \in \text{box}_n(x) ~\text{and}~ a(t_2) \in \text{box}_n(y)\}$$ The latter one would indicate a convergence towards $f_{X|\Theta} f_\Theta$.

Furthermore, there are some more problems with this definition, namely that the histograms produce values between $0$ and $1$ while a density can have values $> 1$. Maybe one needs a constant $c_n$ in front of $\frac{|\{t \in \{1,...,n\} : a(t) \in \text{box}_n(x)\}}{n}$? Moreover, what happens if the histogram becomes finer faster than the sampling process produces samples, like when I divide $\mathbb{R}$ into pieces of size $2^{-n}$ or so then I get exponentially many boxed but only linearly many samples... maybe one should have two limits? LIke so

$$\lim_{m \rightarrow \infty} \lim_{n \rightarrow \infty} c_n \frac{|\{t \in \{1,...,n\} : a(t) \in \text{box}_m(x)\}}{n} = f(x)$$

(??)

Concrete example:

$X$ is a binary variable mapping to $\{0,1\}$ and is Bernoulli distributed with Parameter $\theta$. $\Theta$ is Beta-distributed with fixed Parameters $\alpha_0 = 3/2, \beta_0 = 2$. This means that we actually model the conditional density

$$f_{X|\Theta}(x, \theta) = \theta^x (1-\theta)^{1-x}$$ and $$f_\Theta(\theta) = \frac{\Gamma(\alpha_0 + \beta_0)}{\Gamma(\alpha_0)\Gamma(\beta_0)} \theta^{\alpha_0 - 1} (1-\theta)^{\beta_0 - 1}$$

We first sample a value $\theta$ from $f_\Theta$ and then we insert this into $f_{X|\Theta}(\cdot, \cdot)$ and then we sample from $f_{X|\Theta}(\cdot, \theta)$. The following R-code produces the picture: par( mfrow = c( 1, 1 ) )
alpha_0 = 1.5
beta_0 = 2

density = function(x) {
return(dbeta(x, alpha_0, beta_0))
}
curve(density, from = 0, to=1)

thetaTest = rbeta(1000000, alpha_0, beta_0)
hist(thetaTest, breaks = 100)

xData = c()
for (i in 1:100000) {
theta = rbeta(1, alpha_0, beta_0)
x = sample(c(TRUE, FALSE), 1, prob = c(theta, 1-theta))

xData = c(xData, x)
}

f_X_Theta = function(x, theta) {
return(theta^x * (1-theta)^(1-x) * dbeta(theta, alpha_0, beta_0))
}

f_X_Theta_x_zero = function(theta) {return(f_X_Theta(0, theta))}
f_X_Theta_x_one = function(theta) {return(f_X_Theta(1, theta))}

par( mfrow = c( 2, 2 ) )

curve(f_X_Theta_x_zero, from=0, to=1)
curve(f_X_Theta_x_one, from=0, to=1)