Help select prior function

Question

I have a certain population of individuals, each of which produces a certain 1D output described by a random variable $X$. To a reasonable approximation, for an individual, $X$ is normal-distributed with some unknown mean and standard deviation (STD). The goal is to estimate those mean and STD on an individual level. The problem is that individuals sometimes produce only a dozen or so datapoints, and there is a risk of severely underestimating the STD because the sample is too small to be representative. I would prefer to slightly overestimate the STD, rather than underestimate it. My idea is to introduce somewhat conservative priors for the population mean and STD, and find a posterior for mean and STD which is a compromise between the prior and the likelihood. Since performance is a concern, I would prefer an analytical solution for the posterior. Since I do not need iterative updating, I do not have to have a conjugate prior, but merely a prior that would allow me to find an analytical expression for the posterior. Further, I do not need a perfect solution. The ultimate goal is to get the marginal posterior distribution for the individual mean, and thus make a qualitative estimate of its uncertainty via, e.g. HPDI.

In order to find a reasonable guess for the prior, I have computed the distribution of the sample means and STDs over the individuals, for which a large number of samples was available. It looks like this:

Originally I was hoping to fit it with normal-gamma distribution, but, as far as I understand, this model does not account for correlation between mean and STD which is clearly present.

Could you recommend a 2D prior function for (mu, std) that would be an ok fit for the above distribution, and that would result in an analytic posterior given gaussian likelihood? I am not 100% married to the gaussian likelihood either. I would at least consider another unimodal likelihood if it allows me to get to an analytical posterior.

EDIT 1: The data is guaranteed to be positive.

EDIT 2: I have now tried the advice of taking the logarithm and it works like a charm - the mean and the variance are no longer correlated. Log-transformed data seem to more or less fit to the normal-gamma 2D distribution, although the mean has slightly long tail. Perhaps gamma-gamma 2D would be better, but since it is not that easy to invert I'll live with normal-gamma, since it is not too bad really.

Answer 1

This is a very common use-case for hierarchical models. These models have a lot of different names, so you'll sometimes see them called multi-level models or mixed-effects models. Hierarchical models aren't strictly Bayesian, but, in my experience, they are more common fit as Bayesian models than Frequentist.

The basic idea is that you define population parameters then have member of the population distributed around the population parameters with their own variability.

In your example, let $i$ denote the person and $j$ denote the measurement. $X_{ij}$ refers to measurement $j$ on person $i$. Each person does not need to have the same number of measurements, $j$. Each person has their own mean, $\theta_i$, and variance $\sigma^2_i$. The person-level means are distributed around an overall mean, $\mu$.

• Individual Level: $X_{i, 1}, X_{i ,2}, \dots, X_{i, j} \sim N(\theta_i, \sigma^2_i)$ for $i = 1, \dots, n$
• Population Level: $\theta_1, \dots \theta_n \sim N(\mu, \tau) \\ \mu \sim N(0, 10), \; \tau \sim U[0, 10] \; \text{(or some other noninformative priors)}$.

You'll be to get a posterior distribution of each person-level mean, as well as the population level mean and variance (or precision).

Hierarchical models are too big a topic to explain in a Stack Exchange post. You can look into some resources to help you fit these models.