A statistical conundrum. Combining two distributions I have a statistical problem in a domain that I can not talk about due to an NDA. However, I have worked out how to describe an exactly analogous problem as follows.
You are in charge of record keeping at a factory that makes girders. Your job is to measure the strength of every girder that comes out of a factory. You have kept a record of tens of thousands of girders over many years and have noticed that the strengths (in whatever units) follow a normal distribution with a certain mean and standard deviation.
You also notice that every girder has the name of the worker that made it stamped on one end. There have been several tens of workers names you've noticed over the years. So far you have never done any statistical tests to see whether the strength of the beams has any correlation with the worker that made it. Indeed it may be possible that any variation in strength is entirely due to variations in the raw materials and nothing to do with the worker at all... you never checked. Indeed you have no idea how the girder is made. Its a secret process the factory owners don't want to tell you.
Now you recently noticed that a new name appeared on the end of a handful of girders, say it's "Fred". You measure their strengths and find they appear stronger than average, but the sample size is small. The question I have now is "What is the likely strength of Fred's next girder?" Presumably it must be some combination of the two distributions. The large overall distribution combined somehow with the small sized distribution... and presumably the weighting given to the smaller distribution should be determined by the degree of correlation between the choice of worker and the girder strength (which we don't know yet - but presumably we can do some test to find out).
So I know the gist of what the answer might look like, but I don't know exactly how to do the math.
Can anyone either tell me the answer, or point me to some resource that solves an analogous problem?
EDIT: I thought of a kind of very inelegant, botched, solution that I shall describe below to suggest that an answer must be possible:
Let us assume that the answer is a combination of the small "Fred" distribution (size f) and the much larger overall distribution (size N). Say that we invent a weighting factor K which expresses how much attention we pay to Fred's small distribution as opposed to the overall distribution... something like:
Strength of Fred's next comes from a combined distribution of K x Fred's distribution + (1-K) x overall distribution.
Then we can do an experiment where for each other worker w we select f girders at random, and derive a small distribution. Then we note the mean and standard deviation of that small distribution. Then we select one new randomly chosen girder from the rest of that worker's production (analogous to choosing Fred's next girder) and note its strength Sw. Then we ask the question "what value of K" would result in the combined distribution most likely to produce Sw for each worker... armed with this K, we can now solve our original question.
I think this proves that there is a solution. Now I just want a more elegant one.
 A: I agree a Bayesian approach makes sense. Let me provide some detail. 
Forecasts for two models
Let $x_i$ denote the strength of one of Fred's girders and let $x_{1:n} = (x_1, \ldots, x_n)$ denote the current set of observations. We assume
\begin{equation}
p(x_i|\theta) = \textsf{N}(x_i|\mu,\sigma^2),
\end{equation} 
where $\theta = (\mu,\sigma^2)$. Our prediction for $x_i$ depends on what we know about $\theta$. 
We will distinguish between two sets of assumptions, $A_1$ and $A_2$, which may be thought of as two models. According to $A_1$ we know $\theta$ (from the the historical distribution), while according to $A_2$ the value of $\theta$ is not completely known. 
According to $A_1$, the forecast for $x_{n+1}$ is independent of $x_{1:n}$ because the value of $\theta$ is known. Let $\widetilde\theta = (\widetilde\mu,\widetilde\sigma^2)$ denote the known value. Then
\begin{equation}
p(x_{n+1}|x_{1:n},A_1) = p(x_{n+1}|\widetilde\theta) = \textsf{N}(x_{n+1}|\widetilde\mu,\widetilde\sigma^2).
\end{equation}
By contrast, according to $A_2$ the knowledge about $\theta$ is summarized by a probability distribution that incorporates what has been learned from the current set of observations (the posterior distribution), building on was surmised about $\theta$ before any such observations were available (the prior distribution).
According to $A_2$, the posterior distribution for $\theta$ is
\begin{equation}
p(\theta|x_{1:n},A_2) = \frac{p(x_{1:n}|\theta)\,p(\theta|A_2)}{p(x_{1:n}|A_2)},
\end{equation}
where $p(\theta|A_2)$ is the prior distribution. 
(We address how to formulate the prior distribution below.) 
The likelihood for $\theta$ embodies the independence of the observations conditional on $\theta$:
\begin{equation}
p(x_{1:n}|\theta) = \prod_{i=1}^n p(x_i|\theta).
\end{equation}
The expression in the denominator is the normalizing constant, also known as the marginal likelihood:
\begin{equation}
p(x_{1:n}|A_2) = \int p(x_{1:n}|\theta)\, p(\theta_2|A_2)\,d\theta.
\end{equation}
Having characterized the information regarding $\theta$ given the current set of observations, we can address the question of the prediction for the next observation. In particular, the predictive distribution is
\begin{equation}
p(x_{n+1}|x_{1:n},A_2) = \int p(x_{n+1}|\theta)\, p(\theta|x_{1:n},A_2)\,d\theta.
\end{equation}
Combining the forecasts
We can use Bayesian Model Averaging (BMA) to combine the predictive distributions from both sets of assumptions. The assumption behind BMA is that all of the observations come form one of the two models, but we do not know which. Therefore, we weight the forecasts by the probabilities we assign to each of the models. 
We use Bayes' rule to compute the posterior model probabilities:
\begin{equation}
p(A_j|x_{1:n}) = \frac{p(x_{1:n}|A_j)\,p(A_j)}{\sum_{j=1}^2 p(x_{1:n}|A_j)\,p(A_j)},
\end{equation}
where $p(A_j)$ is the prior model weight for model $j$. 
The likelihood of the observations $x_{1:n}$ according to $A_1$ is 
\begin{equation}
p(A_1|x_{1:n}) = \prod_{i=1}^n p(x_i|\theta_1)
\end{equation}
and the likelihood of the observations according to $A_2$ is given above. 
Given the posterior model probabilities, we average the forecasts to produce a combined predictive distribution:
\begin{equation}
p(x_{n+1}|x_{1:n}) = p(x_{n+1}|x_{1:n},A_1)\,p(A_1|x_{1:n}) + p(x_{n+1}|x_{1:n},A_2)\, p(A_2|x_{1:n}).
\end{equation}
In effect, we have "integrated" the two models out of the predictive distribution.
The posterior model probabilities can be expressed in ratio terms:
\begin{equation}
\frac{p(A_1|x_{1:n})}{p(A_2|x_{1:n})} = \frac{p(A_1)}{p(A_2)}\, \frac{p(x_{1:n}|A_1)}{p(x_{1:n}|A_2)},
\end{equation}
where the second factor on the right-hand side is called the Bayes factor. Since $A_1$ is nested within $A_2$ the Bayes factor can be expressed as the Savage-Dickey density ratio:
\begin{equation}
\frac{p(x_{1:n}|A_1)}{p(x_{1:n}|A_2)} = \frac{p(\widetilde\theta|x_{1:n},A_2)}{p(\widetilde \theta|A_2)}.
\end{equation}
We will have a closed-form expression for the combined predictive distribution if we have closed-form expressions for the posterior and the posterior predictive under $A_2$: $p(\theta|x_{1:n},A_2)$ and $p(x_{n+1}|x_{1:n},A_2)$. 
A prior for $\theta$
We now address the question as to how to formulate the prior for $\theta$ under the second set of assumptions. 
In the absence of information to the contrary, we adopt the conjugate prior. 
With this in mind, define the normal-inverse-gamma conjugate distribution:
\begin{equation}
\textsf{Normal-IG}(\mu,\sigma^2|m, \kappa, s^2, \nu) = \textsf{N}(\mu|m, \sigma^2/\kappa)\, \textsf{Inv-Gamma}(\sigma^2|\nu/2, s^2\nu/2),
\end{equation}
where
\begin{equation}
\textsf{Inv-Gamma}(\sigma^2|\nu/2, s^2\nu/2) = 
\frac{\left(s^2\nu/2\right)^{\nu /2}}{\Gamma(\nu/2)}\, \frac{e^{-\frac{s^2\nu/2}{\sigma^2}}}{\left(\sigma^2\right)^{\frac{\nu }{2}+1}}.
\end{equation}
The predictive distribution is 
\begin{equation}
\iint \textsf{N}(x|\mu,\sigma^2)\, \textsf{Normal-IG}(\mu,\sigma^2|m, \kappa, s^2, \nu)\,d\mu\,d\sigma^2 = \textsf{Student}(x|m,s^2/\kappa,\nu).
\end{equation}
Let the prior for $\theta$ be given by
\begin{equation}
p(\theta|A_2) = p(\mu,\sigma^2|\widetilde\mu, \underline\kappa, \widetilde\sigma^2, \underline s^2).
\end{equation}
This distribution is "centered" on $(\underline\mu,\underline\sigma^2) = (\widetilde\mu,\widetilde\sigma^2)$. There are two remaining free parameters: $\underline\kappa$ and $\underline\nu$. These parameters have the interpretation of the number of "prior" observations pertaining to $\mu$ and $\sigma^2$, respectively. The values of these parameters can affect the posterior model probabilities, and thus one should give some thought to the values chosen for them. To this end, note that
\begin{equation}
p(x_1|A_2) = \textsf{Student}(x_1|\widetilde\mu, \widetilde\sigma^2/\underline\kappa, \underline\nu).
\end{equation}
Ask yourself what values of $\underline\kappa$ and $\underline\nu$ characterize your uncertainty about the distribution for the strength of Fred's girders (before you observed any of Fred's output). 
One possible setting is $\underline\kappa = \underline \nu = 1$. This setting produces a Cauchy distribution with location $\widetilde\mu$ and scale $\widetilde\sigma$, where
\begin{equation}
\textsf{Cauchy}(x|\mu,\sigma^2) = \frac{1}{\pi\,  \sigma  \left(1+\frac{(x-\mu)^2}{\sigma^2}\right)}.
\end{equation}
(The Cauchy distribution does not have a mean or a variance.) In thinking about this prior, consider the density ratio
\begin{equation}
\frac{\textsf{Cauchy}(x|\mu,\sigma^2)}{\textsf{N}(x|\mu,\sigma^2)}.
\end{equation} 
The density ratio is greater that one unless $|x-\mu| \le 1.851\,\sigma$.
The posterior distribution for $\theta$ is given by
\begin{equation}
p(\theta|x_{1:n},A_2) = \textsf{Normal-IG}(\mu,\sigma^2|\overline\mu,\overline\kappa,\overline s^2, \overline \nu),
\end{equation}
and hence the posterior predictive distribution for $x_{n+1}$ is given by
\begin{equation}
p(x_{n+1}|x_{1:n},A_2) = \textsf{N}(x_{n+1}|\overline\mu, \overline\sigma^2/\overline\kappa,\overline\nu),
\end{equation}
where
\begin{align}
\overline\nu &= \underline\nu + n \\
\overline\kappa &= \underline\kappa + n \\
\overline\mu &= (\underline \kappa/\overline \kappa)\, \widetilde\mu + (n/\overline \kappa)\, \widehat\mu \\
\overline\sigma^2 &= (\underline\nu/\overline \nu)\,\widetilde\sigma^2 + (n/\overline\nu)\, \widehat\sigma^2 + (\underline\kappa/\overline\kappa)\,(n/\overline\nu)\,(\widehat\mu-\widetilde\mu)^2,
\end{align}
where
\begin{equation}
\widehat\mu = \frac{1}{n}\sum_{i=1}^n x_i 
\qquad\text{and}\qquad
\widehat\sigma^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \widehat\mu)^2.
\end{equation}
You now have closed-form expressions for all the parts of the combined predictive distribution. 
A: For a more general solution than mef's answer that allows you to make predictions for any worker who's already made a girder (not just Fred in particular), you could use a hierarchical model.
Let $x$ denote girder strength, let $i$ index girders, and let $j$ index workers. Then:
$$\begin{align}
x_{ij} &\sim N\left( \mu_i , \sigma^2 \right)\\
\mu_i &\sim N\left( \nu , \tau^2 \right)
\end{align}$$
where $\sigma^2$, $\nu$, and $\tau^2$ are unknown parameters.
There is no harm in assuming normality here because $x_i \sim N$, and it was already found that $p(x_i) = \sum_j p(x_{ij}) p(j)$. That is, the density of all girder strengths is the sum of the density of girder strengths for each worker, weighted by the probability of that worker making a given girder. It's a well-known result that a mixture of normal distributions is itself normal.
Then the complete posterior density for the parameters is
$$
p \left( \mu , \sigma^2 , \nu , \tau^2\ |\ x \right) = \prod_j \prod_i \phi \left( x_{ij}\ |\ \mu_i , \sigma^2 \right) \phi \left( \mu_j\ |\ \nu, \tau^2 \right) p \left( \nu, \sigma^2, \tau^2 \right)
$$
where $p \left( \nu, \sigma^2, \tau^2 \right)$ is the joint prior distribution; generally you'd assume that $p \left( \nu, \sigma^2, \tau^2 \right) = p \left( \nu \right) p \left( \sigma^2 \right) p \left( \tau^2 \right)$. For $\nu$, you'd be safe specifying a normal prior with a large variance. For $\sigma^2$ and $\tau^2$, you could look into the "folded noncentral t family" described in "Prior distributions for variance parameters in hierarchical models" by Gelman (2006).
For an example of how this model works in practice, take a look at "Multilevel (hierarchical) modeling: what it can and can't do", by Gelman (2005).
Disclaimer: I don't have any experience using these models outside of the classroom. But I did take a class and a half with Prof Gelman and another with a postdoc who works for him.
