Combine three equally likely observed PDFs describing carrot sizes in a carrot farm I observe things so that, during an observational session, I develop a probability density function $f(x)$ describing the likelihood that $x$ is equal to different real numbers. I have $n$ of these observed probability function $[f_1,f_2,\dots,f_n]$. All of these observations seem equally valid to me. How can I combine these into a single probability function, describing the likelihood that $x$ is really different numbers given all of these data?
Another way to think about it is with a somewhat contrived example. Imagine I want to know the distribution of the sizes of carrots currently in my organic carrot farm, so I can advertise statistics about how big they are. In fact, I want to know the exact shape of the carrot size distribution. My farm is huge, so I can't measure them all. I choose three parcels of land and measure 1000 carrots in each parcel, only recording the density functions of carrot sizes found during each measuring session. Since this is all the information I'm willing to invest in, I have to estimate a probability distribution for the entire farm using only this information. I think all three of the measured density functions are equally likely to describe the true carrot size distribution.
Here is a coded example for $n=3$. Some commenters suggested I try averaging; it is not clear to me if that is the right approach. If I choose to average them, I have two options. I can divide the $x$ axis into small rectangles and find the average $f$ value for every $\Delta x$, or divide the $f$ axis into small rectangles and find the average $x$ for each $\Delta f$. I want to rigorously approach this problem, and it is not clear to me if one of these methods is better than the other, if that is at all a valid approach, or if there is a better way. My reasoning in exploring convolution is that the average carrot distribution would be $f_1$ + $f_2$ + $f_3$ normalized such that the area under the curve is one. Based on the community's feedback, I now think that this calls for bayesian updating math that I am not sure how to do.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# Make observations. We only have PDFs, not data. We loose some information
#    by assuming they are normally distributed.
f1 = np.random.normal(0, 2, 100) + np.random.normal (1, 4, 100)
f2 = np.random.normal(0, 2, 100) + np.random.normal (1, 4, 100)
f3 = np.random.normal(0, 2, 100) + np.random.normal (1, 4, 100)
mu1, std1 = norm.fit(f1)
mu2, std2 = norm.fit(f2)
mu3, std3 = norm.fit(f3)
del f1
del f2
del f3

# plot distributions
x = np.arange(-10, 10, .1)
f1 = norm.pdf(x, mu1, std1)
f2 = norm.pdf(x, mu2, std2)
f3 = norm.pdf(x, mu3, std3)
plt.plot(x, f1, label='f1')
plt.plot(x, f2, label='f2')
plt.plot(x, f3, label='f3')
plt.legend()
plt.show()


What I want is to estimate the "true" distribution that I took these samples from. Here's what that should look like:
import numpy as np
import matplotlib.pyplot as plt
nsamps = 20000
truth = np.random.normal(0, 2, nsamps) + np.random.normal (1, 4, nsamps)
plt.hist(truth, 400, normed=1)
plt.show()


 A: In your example you're separately fitting normal distributions to three samples (resulting from the same data-generating process), & asking how to combine these fits into a single fitted distribution; as if you'd had all the observations in one large sample to start with.
The mean is doubtless estimated from the $n_j$ observations of the $j$th sample, the $x_j$s, as
$$\hat\mu_j = \frac {\sum_{i=1}^{n_j} x_{ij}}{n_j}$$
& the standard deviation as
$$\hat\sigma_j = \sqrt{\frac{\sum_i^{n_j} x_{ij}^2 - \frac{\left(\sum_i^{n_j} x_{ij}\right)^2}{n_j}}{d(n_j)}}$$
where $d(n_j)$ is some function of the sample size (commonly $n-1$, equating to an unbiased estimator of the variance). If you've kept all the observations  you can estimate the mean, by the same procedure, as
$$\hat\mu = \frac {\sum_{j=1}^3\sum_{i=1}^{n_j} x_{ij}}{\sum_{j=1}^3 n_j}$$
& the standard deviation as 
$$\hat\sigma = \sqrt{\frac{\sum_{j=1}^3\sum_{i=1}^{n_j} x_{ij}^2 - \frac{\left(\sum_{j=1}^3\sum_i^{n_j} x_{ij}\right)^2}{\sum_{j=1}^3 n_j}}{d\left(\sum_{j=1}^3 n_j\right)}}$$
But if you haven't it doesn't matter, because it's straightforward to recover $\sum_{i=1}^{n_j} x_{ij}$ & $\sum_{i=1}^{n_j} x_{ij}^2$ from the parameter estimates for each sample.
The above goes for method-of-moments estimation in general, not just for fitting normal distributions.
It's interesting you mention information loss—in fact this is a case where you're not losing information about the parameters by reducing the data to a vector of parameter estimates. When the parameter estimates are jointly sufficient you can get a sufficient estimate for a pooled sample using just the individual sample estimates. So, for example, with a gamma distribution you can  recover $\sum x$ & $\sum \log x$ from the maximum-likelihood estimates of scale & shape. But if you were estimating the location parameter of a Laplace distribution by maximum likelihood, you'd have to keep the order statistic at least because the maximum-likelihood estimate, the median, is not sufficient.
