# get a bunch of imports out of the way
import matplotlib.pyplot as plt
plt.rc('text', usetex=True)
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'
import numpy as np
import scipy.stats as st
%matplotlib inline

M = st.uniform.rvs(1.0, 100.0, size=10)
F = np.sqrt(M)
mu = 2*F
N = st.poisson.rvs(mu)

Mgrid = np.linspace(M.min(), M.max())
mu_of_Mgrid = 2*np.sqrt(Mgrid)
plt.plot(Mgrid, mu_of_Mgrid, label=r'$\mu(M)$');
plt.plot(M, N, 'o', label='data');
plt.xlabel(r'$M$', fontsize='x-large');
plt.ylabel(r'$N$', fontsize='x-large');
plt.legend(fontsize='x-large');

"""
Explanation: Notes: So, about those error bars...
In which we will make sense of the common delusion of error bars.
The Poisson example we've used so far when thinking about generative models and Bayes' Law is one where the idea of an error bar doesn't really apply. We measure a certain number of counts, period. That's it. It doesn't make sense to say that we measured, as the case may be, 5 $\pm$ something (even more so if that something is not an integer), because literally what we measured was the number 5.
But error bars are enough of a thing that even the most pedantic among us will talk about them fairly often. So what's the deal? Let's look at a contrived example to sketch it out.
Specifically, let's expand on the previous Poisson example by imagining we (a) measure a number of counts (with the same telescope and exposure time) from many stars with different masses, and (b) live in magical universe where stars have name tags that tell us their precise masses. We're interested in how luminosity, which is directly proportional to the expectation value for number of counts, relates to this magically known mass. In generative terms,
* masses $M_i$ are known somehow,
* fluxes $F_i \Leftarrow$ $M_i$,
* expected counts $\mu_i \Leftarrow F_i$,
* measured counts $N_i \sim \mathrm{Poisson}(\mu_i)$.
For simplicity, let's implement the details as
* Choose $M_i$,
* Let $F_i=M_i^{0.5}$,
* Let $\mu_i = 2F_i$.
Here's a made-up data set.
End of explanation
"""


Nlower = st.poisson.ppf(0.16, mu_of_Mgrid)
Nupper = st.poisson.ppf(0.84, mu_of_Mgrid)
plt.plot(Mgrid, mu_of_Mgrid, label=r'$\mu(M)$');
plt.plot(M, N, 'o', label='data');
plt.plot(Mgrid, Nlower, '--', color='C0');
plt.plot(Mgrid, Nupper, '--', color='C0');
plt.xlabel(r'$M$', fontsize='x-large');
plt.ylabel(r'$N$', fontsize='x-large');
plt.legend(fontsize='x-large');

"""
Explanation: Here the curve shows the Poisson mean as a function of $M$. Clearly, the data don't sit on the curve, nor should they. But it would be nice to represent the width of the sampling distribution somehow, so we can see how compatible the data are with the model. Among other possibilities, we might do that by also showing curves reflecting the 16th and 84th percentiles of the sampling distribution, in addition to the mean, as a function of $M$ (so that the probability between the curves, at fixed $M$, is 68%). That would look like this:
End of explanation
"""


M = st.uniform.rvs(1000.0, 10000.0, size=10)
F = np.sqrt(M)
mu = 2*F
N = st.poisson.rvs(mu)

Mgrid = np.linspace(M.min(), M.max())
mu_of_Mgrid = 2*np.sqrt(Mgrid)
Nlower = st.poisson.ppf(0.16, mu_of_Mgrid)
Nupper = st.poisson.ppf(0.84, mu_of_Mgrid)
plt.plot(Mgrid, mu_of_Mgrid, label=r'$\mu(M)$');
plt.plot(M, N, 'o', label='data');
plt.plot(Mgrid, Nlower, '--', color='C0');
plt.plot(Mgrid, Nupper, '--', color='C0');
plt.xlabel(r'$M$', fontsize='x-large');
plt.ylabel(r'$N$', fontsize='x-large');
plt.legend(fontsize='x-large');

"""
Explanation: You can see that about $2/3$ of the data lie within these limits. The jaggedness of the lines is due to the fact that, of course, we can only get integers out of the Poisson distribution.
Let's now take this to a regime where there would be many more counts in each measurement.
End of explanation
"""


F_of_Mgrid = mu_of_Mgrid / 2.
Flower = (mu_of_Mgrid - np.sqrt(mu_of_Mgrid)) / 2.
Fupper = (mu_of_Mgrid + np.sqrt(mu_of_Mgrid)) / 2.
Fhat = N / 2.

plt.plot(Mgrid, F_of_Mgrid, label=r'$\hat{F}(M)$');
plt.plot(M, Fhat, 'o', label='data');
plt.plot(Mgrid, Flower, '--', color='C0');
plt.plot(Mgrid, Fupper, '--', color='C0');
plt.xlabel(r'$M$', fontsize='x-large');
plt.ylabel(r'$\hat{F}$', fontsize='x-large');
plt.legend(fontsize='x-large');

"""
Explanation: We're now well into the limit where the Poisson (mean $\mu$) distribution begins to resemble the Gaussian distribution (mean $\mu$ and standard deviation $\sqrt \mu$). As we're looking at many more counts, the discreteness of the dashed lines is also harder to see.
Now, we might think about displaying this with $F$ rather than the actual measurement, $N$, on the $y$ axis. Of course, we know that the plotted points don't actually correspond to the true flux of each star, only the naive extimate we would make by plugging through the equations at the top of the notebook and ignoring any uncertainty in the model. So let's call it "estimated flux", $\hat{F}$, instead. We'll use the fact that we're now in the Gaussian limit to simply compute the dashed limits on the model prediction as the mean $\pm$ the standard deviation.
End of explanation
"""


Fhat_err = np.sqrt(mu) / 2.
plt.plot(Mgrid, F_of_Mgrid, label=r'$\hat{F}(M)$');
plt.errorbar(M, Fhat, yerr=Fhat_err, fmt='o', label='data');
plt.xlabel(r'$M$', fontsize='x-large');
plt.ylabel(r'$\hat{F}$', fontsize='x-large');
plt.legend(fontsize='x-large');

"""
Explanation: In this simple example, apart from the tiny change to the dashed lines, this just corresponds to a rescaling of the $y$ axis.
Finally, we could (if we wanted to) remember the form of the Gaussian distribution,
$\mathrm{Normal}(x|\mu,\sigma) = \frac{1}{\sqrt{2\pi}\sigma} \exp\left[ -\frac{(x-\mu)^2}{2\sigma^2} \right]$.
Because this density depends on $x$ and $\mu$ only through the distance between them, if we're using $\sigma$ (or some other scale) to represent model uncertainty, we could choose to position that visual cue with the data points rather than with the model. The distance between the data and model in units of that shown uncertainty would, of course, be the same either way. Note that the same statement could not be made when we were in the small-$N$ regime, or at least not made as simply.
End of explanation
"""


