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:

enter image description here

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.

  • 2
    $\begingroup$ What would things look like, if you worked on the log-scale (assuming these values can never be 0 or negative)? The distribution looks a bit skewed like it could be log-normal (or something else, of course) and log-transforming might even reduce/get rid of the mu to SD relationship. $\endgroup$
    – Björn
    Commented Dec 23, 2022 at 14:32
  • $\begingroup$ Thanks! Indeed, transforming the data to try to get rid of correlation sounds like the simplest approach. I'll try logarithm and maybe a few other transformations. $\endgroup$ Commented Dec 24, 2022 at 16:06
  • 1
    $\begingroup$ @Fiodor1234 I'm sorry, I am not too familiar with terminology (i.e. I don't know what 'vector representation' is, sorry). There are $m$ individuals, and each individual produces a different number $n_i$ of datapoints $X_{i,j} \sim \mathcal{N}(\mu_i, \sigma_i)$. I want to infer the individual parameters $\mu_i$ and $\sigma_i$. I know that they themselves are sampled from the prior distribution plotted above. I actually don't care about $\sigma_i$ and only want the marginal posterior $P[\mu_i| Data]$ $\endgroup$ Commented Dec 24, 2022 at 17:19
  • 1
    $\begingroup$ Unfortunately, you're pretty limited if you want an analytic prior. If the lognormal recommendation of @Bjorn above doesn't work, you might want to consider non-analytic approaches, e.g., MCMC, which will really open up the possibilities. $\endgroup$
    – jbowman
    Commented Dec 28, 2022 at 20:43
  • 1
    $\begingroup$ @AleksejsFomins, I asked about multilevel models because this is the canonical example of when to use them. I thought there might be some reason why they wouldn't work that I wasn't seeing. I'll write something up or at least provide a link. $\endgroup$
    – Eli
    Commented Dec 28, 2022 at 22:24

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.