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()


• Your question is unclear in several ways, not least that you seem to conflate density with probability (they're different). You seem to be observing densities. Can you explain how they arise? Why would you convolve them rather than take an equal-weight mixture distribution (i.e. why not actually average the densities)? Your question might be more readable with slightly smaller plots (placing m directly before .png will make them half as big each way but perhaps that's too much smaller) – Glen_b -Reinstate Monica May 25 '17 at 6:05
• I think you're confused about what a convolution is - it's the distribution of a sum of random variables. A mixture on the other hand has a density function that's the weighted sum of other density functions - an average. Which is the correct one to use depends on what you're trying to do, which I don't quite get from your question. – Scortchi - Reinstate Monica May 25 '17 at 19:11
• @Scortchi I added an explanatory example where I own a carrot farm and want to know the probability distribution of carrot sizes – kilojoules May 25 '17 at 19:41
• 1. The Q is improved but I think will require some further edits. I'll reopen as amoeba asks. 2. You say "I spend [...] sessions creating probability density functions of what I think the likelihood of x being different values is" -- it's precisely this that I wanted explained. If you're giving probabilities of values, you don't have a density, so I want to know what exactly you are doing and it may be necessary to know more about the process of doing it in order to arrive at a probability model for it. 3. while the Q treats density and probability as the same it remains unclear in that sense – Glen_b -Reinstate Monica May 25 '17 at 20:56
• @Glen_b I think your questions show that you suspect that OP does not well appreciate & describe what their real application might be about, and this is fair enough. But OP's question as illustrated by the carrot farm example (and by the code) is pretty clear, no? In each "session" some data are collected, empirical PDF estimated (by some kind of smoothing method), then data are discarded. Assuming that the true PDF did not change between sessions, the question is, how to estimate it by combining session-wise PDFs in some manner. – amoeba May 25 '17 at 21:31

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