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