Figure 3.5

Asymmetric distributions of astronomical magnitudes

In [1]: 
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [3]: 
%load http://www.astroml.org/_downloads/fig_flux_errors.py
In [7]: 
"""
Flux Errors
-----------
Figure 3.5.

An example of Gaussian flux errors becoming non-Gaussian magnitude errors.
The dotted line shows the location of the mean flux; note that this is not
coincident with the peak of the magnitude distribution.
"""
# Author: Jake VanderPlas
# License: BSD
#   The figure produced by this code is published in the textbook
#   "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
#   For more information, see http://astroML.github.com
#   To report a bug or issue, use the following forum:
#    https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import norm

#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX.  This may
# result in an error if LaTeX is not installed on your system.  In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=16, usetex=True)

#------------------------------------------------------------
# Create our data

# generate 10000 normally distributed points
np.random.seed(1)
dist = norm(1, 0.25)
flux = dist.rvs(10000)
flux_fit = np.linspace(0.001, 2, 1000)
pdf_flux_fit = dist.pdf(flux_fit)

# transform this distribution into magnitude space
mag = -2.5 * np.log10(flux)
mag_fit = -2.5 * np.log10(flux_fit)
pdf_mag_fit = pdf_flux_fit.copy()
pdf_mag_fit[1:] /= abs(mag_fit[1:] - mag_fit[:-1])
pdf_mag_fit /= np.dot(pdf_mag_fit[1:], abs(mag_fit[1:] - mag_fit[:-1]))

#------------------------------------------------------------
# Plot the result
fig = plt.figure(figsize=(10, 5.0))
fig.subplots_adjust(bottom=0.17, top=0.9,
                    left=0.12, right=0.95, wspace=0.3)

# first plot the flux distribution
ax = fig.add_subplot(121)
ax.hist(flux, bins=np.linspace(0, 2, 50),
        histtype='stepfilled', fc='gray', alpha=0.5, normed=True)
ax.plot(flux_fit, pdf_flux_fit, '-k')
ax.plot([1, 1], [0, 2], ':k', lw=1)
ax.set_xlim(-0.1, 2.1)
ax.set_ylim(0, 1.8)

ax.set_xlabel(r'${\rm flux}$')
ax.set_ylabel(r'$p({\rm flux})$')
ax.yaxis.set_major_locator(plt.MultipleLocator(0.4))
ax.text(0.04, 0.98, r'${\rm 20\%\ flux\ error}$',
        ha='left', va='top', transform=ax.transAxes,
        bbox=dict(ec='none', fc='w'))

# next plot the magnitude distribution
ax = fig.add_subplot(122)
ax.hist(mag, bins=np.linspace(-1, 2, 50),
        histtype='stepfilled', fc='gray', alpha=0.5, normed=True)
ax.plot(mag_fit, pdf_mag_fit, '-k')
ax.plot([0, 0], [0, 2], ':k', lw=1)
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(0, 1.8)
ax.yaxis.set_major_locator(plt.MultipleLocator(0.4))
ax.text(0.04, 0.98, r'${\rm mag} = -2.5\log_{10}({\rm flux})$',
        ha='left', va='top', transform=ax.transAxes,
        bbox=dict(ec='none', fc='w'))

ax.set_xlabel(r'${\rm mag}$')
ax.set_ylabel(r'$p({\rm mag})$')

plt.show()

Understanding this figure

In [8]: 
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import norm
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=16, usetex=True)
In [12]: 
#------------------------------------------------------------
# Create our data

# generate 10000 normally distributed points
np.random.seed(1)
dist = norm(1, 0.25) 
flux = dist.rvs(10000)
flux_fit = np.linspace(0.001, 2, 1000)
pdf_flux_fit = dist.pdf(flux_fit)

Ok, this next part is the transformation from flux to magnitude. Let's sketch out the math here: $$ p(y) = p(x) (dy / dx)^{-1} $$ $$ x \sim \mathcal{N}(1, \sigma) $$ $$y(x) = -2.5 \log_{10}(x)$$ where y is the magnitude and x is the flux. $$dy / dx = -2.5 \frac{1}{x \ln{10}} $$ so: $$p(y)= \mathcal{N}(1, \sigma) \times \mid -\frac{\ln{10}}{2.5} x \mid $$

In [32]: 
x_sample = flux
x = np.linspace(0.001, 2, 1000)
Px = dist.pdf(x) 

y_sample = -2.5 * np.log10(x_sample)
y = -2.5 * np.log10(x)

dy_dx= -2.5 / ( x * np.log(10)) #natural log

Py= Px / np.abs(dy_dx) #must do absolute value, see Eq. 3.20

mag = -2.5 * np.log10(flux) 
mag_fit = -2.5 * np.log10(flux_fit) 

#pdf_mag_fit = pdf_flux_fit.copy()
#pdf_mag_fit[1:] /= abs(mag_fit[1:] - mag_fit[:-1])
#pdf_mag_fit /= np.dot(pdf_mag_fit[1:], abs(mag_fit[1:] - mag_fit[:-1]))
#pdf_mag_fit = pdf_flux_fit / dy_dx
In [31]: 
#------------------------------------------------------------
# Plot the result
fig = plt.figure(figsize=(10, 5.0), facecolor='#f3f3f3')
fig.subplots_adjust(bottom=0.17, top=0.9,
                    left=0.12, right=0.95, wspace=0.3)

# first plot the flux distribution
ax = fig.add_subplot(121, axis_bgcolor='#f3f3f3')
ax.hist(flux, bins=np.linspace(0, 2, 50),
        histtype='stepfilled', fc='#00ffff', alpha=0.3, normed=True)
ax.plot(x, Px, '-k') #Boom, x and y work...
ax.plot([1, 1], [0, 2], ':k', lw=1)
ax.set_xlim(-0.1, 2.1)
ax.set_ylim(0, 1.8)
ax.set_xlabel(r'${\rm flux}$')
ax.set_ylabel(r'$p({\rm flux})$')
ax.yaxis.set_major_locator(plt.MultipleLocator(0.4))
ax.text(0.04, 0.98, r'${\rm 20\%\ flux\ error}$',  #not 25% flux error?
        ha='left', va='top', transform=ax.transAxes,
        bbox=dict(ec='none', fc='#f3f3f3'))


# next plot the magnitude distribution
ax = fig.add_subplot(122, axis_bgcolor='#f3f3f3')
ax.hist(y_sample, bins=np.linspace(-1, 2, 50),
        histtype='stepfilled', fc='#00FFFF', alpha=0.3, normed=True)

#ax.plot(mag_fit, pdf_mag_fit, '-k')

ax.plot(y, Py, '-k') #Boom shakalaka.  Applying (modified) Equation 3.20 worked, as expected.

ax.plot([0, 0], [0, 2], ':k', lw=1)
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(0, 1.8)
ax.yaxis.set_major_locator(plt.MultipleLocator(0.4))
ax.text(0.04, 0.98, r'${\rm mag} = -2.5\log_{10}({\rm flux})$',
        ha='left', va='top', transform=ax.transAxes,
        bbox=dict(ec='none', fc='#f3f3f3'))

ax.set_xlabel(r'${\rm mag}$')
ax.set_ylabel(r'$p({\rm mag})$')

plt.show()

Recap

OK- here are the key ideas and lessons learned from this figure:

  • Equation 3.20 still works. It's for calculating the scaled pdf. It seems pretty easy, as long as you have a functional form of the transformation, which you should have in hand.
  • This figure used a different approach to calculate $p(y)$. No pencil and paper derivative necessary- I suppose it was a numerical derivative.
  • It really is stunning that the median is no longer the mode of the transformed distribution. It makes sense in this figure. Do people know this?