I would like to sample from: $$ p(\theta_2|x)=\int p(\theta_2|\theta_1,x) . p(\theta_1|x) . d\theta_1 $$ knowing that I can easily sample from $p(\theta_1|x)$ and (less easily) from $p(\theta_2|\theta_1,x)$. I have some “naive-brute-force-computationally-intractable” ideas but maybe there exists some more elegant ones?
EDIT:
1) Several answers suggest me to use the following procedure:
for i=1:N
draw $\theta_1^i$ from $p(\theta_1|x)$
draw $\theta_2^i$ from $p(\theta_2|\theta_1^i,x)$
Is this provide me i.i.d $\theta_2^i$ from $p(\theta_2|x)$ ?
2) In practice, I would like to compute an HPD interval over $\theta_2$ (e.g. using boa.hpd). Thus I do not really need i.i.d RVs. Is the following procedure valid :
for i=1:P
draw $\theta_1^{i}$ from $p(\theta_1|x)$
for j=1:N/P
draw $\theta_2^{(i-1)*P+j}$ from $p(\theta_2|\theta_1^i,x)$
i.e. drawing many $\theta_2^{(i-1)*P+j}$ from a single $\theta_1^{i}$, $P$ times ($P$ being a divisor of $N$ and being fixed and $N$ the number of samples that I want to generate) and use these $N$ samples to compute hpd. In practice such a solution would be largely more efficient than the first one in my application.