# The proper way to compute the posterior distribution of a distribution

Suppose I am a Bayesian working with multi-level data, $$j$$ and $$t$$.

I run a model using $$t$$ that calculates the posterior distribution of a parameter $$\theta_j$$ for each $$j$$, which I then use to calculate the posterior distribution a function $$g(\theta_j)$$ for each $$j$$. Thus, afterwards, if I ran $$R$$ draws in my MCMC, I end up with:

$$\{g(\theta_j^r)\}_{r=1}^{R}$$ for each $$j$$.

Now, I am interested in the variation in $$E_{\theta_j}[g(\theta_j)]$$ across $$j$$ units, which can be characterized by a distribution, which I will call $$f(g)$$.

Of course, my best guess of $$f$$ would be to plot the empirical distribution of the $$E_{\theta_j}[g(\theta_j)]$$'s. However, I want to characterize the uncertainty in $$f$$. How can I do this empirically given all my draws? A bootstrap?

Here is what I tried before I realized this was vastly over-estimating uncertainty:

Basically, I plotted the distribution of $$g$$'s at each value of $$\theta$$. Then took the credible intervals vertically. However, I think this over-estimates the uncertainty and actually gives an incorrect distribution. For instance, I know this because in my problem, the $$E_{\theta_j}[g(\theta_j)]$$ must be strictly positive, and when I plot their empirical distribution, they are. However, when I use the single-draw-at-a-time approach, we get a lot of negative mass with too much precision. Sure, maybe some of the draws might turn out negative, but not the expectation itself. Maybe a bootstrap would solve this? But I was hoping for a more bayesian-theoretically guided approach here.

The problem with this question, upon further reflection, is that there is not a posterior distribution of $$E_{\theta_j}[g(\theta_j))]$$. If you are a true Bayesian, the expectation is just a number at this point. Therefore, it doesn't make sense fo try to characterize uncertainty for each expectation (and by extension, the distribution of the expectations).