Given that the $\tau_j$ are i.i.d. Gamma with parameters $\alpha, \lambda$ (where $f(x) \propto \lambda^{\alpha} x^{\alpha-1}\exp\{-\lambda x\}$, to define the parameterization of the Gamma), then $\sigma_j = 1/\sqrt{\tau_j}$ has a Stacy distribution, which is a special case of the Amoroso distribution. This has a closed form expression for the mean and standard deviation, assuming the parameters are such that the mean and standard deviation exist, in terms of $\alpha$ and $\lambda$, and hence in terms of your m
and d
.
Define $\theta = 1/\sqrt{\lambda}$. Then:
$\mathbb{E}\sigma_j = \alpha + \theta {\Gamma(\alpha-{1\over 2}) \over \Gamma(\alpha)} $
and
$\text{Sd}\space \sigma_j = \theta \left[{\Gamma(\alpha-1) \over \Gamma(\alpha)}-\left({\Gamma(\alpha-{1\over 2}) \over \Gamma(\alpha)}\right)^2 \right]^{1\over2}$
where $\alpha > 1/2$ for the mean to exist and $\alpha > 1$ for the variance to exist. (If you read the article, note the $\beta$ parameter of the distribution $= -2$, as $\tau = \sigma^{-2}$.)
One method for generating MCMC samples of the mean and std. deviation of the $\sigma_j$ would be to take your samples of m
and d
(or, more simply, rG
(= $\lambda$) and sG
(= $\alpha$)), and transform them as above. However, one suspects that you may indeed have samples where $\alpha \leq 1$ or even $\alpha \leq 1/2$. No doubt many solutions to this exist, but two that come to mind are to 1) throw out all the samples (in their entirety) where this happens, which corresponds to altering your prior on m
and d
so that they are confined to a region of $\mathbb{R}^{+2}$ where the mean and variance of the $\tau$ exist, or 2) put your prior on $\sigma$ instead of $\tau$, altering the hyperpriors to something appropriate. Without knowing how strongly you're attached to the particular shape of the priors on m
and d
, though, it's hard to give advice about what "something appropriate" might be. Maybe just put a Gamma distribution on the $\sigma_j$ instead of the $\tau_j$ and continue as you did... hopefully the choice of prior doesn't make much difference to your posterior.
If it doesn't happen all that often, and you like those priors as they are, I'd probably go with option 1. If it happens a lot, though, then those priors in combination with your data are perhaps telling you something about the upper tail of the posterior distribution of the $\sigma_j$.
To summarize why this works, this is equivalent to:
1) Put a Stacy prior on the $\sigma_i$, parameterized by rG
and sG
,
2) Put hyperpriors on rG
and sG
as they currently are,
3) Instead of bothering to generate Stacy-distributed $\sigma_i$ and transforming them to $\tau_i$, just generate the $\tau_i$ directly using the relationship between the Stacy prior on $\sigma$ and the equivalent Gamma prior on $\tau$,
4) Calculate samples of the first two moments of the Stacy prior from the sampled rG
and sG
.