Target Fluctuations

This notebook is not a complete effort, although it does contain several potentially useful plots which show how QSO, ELG, BGS, and LRG target density varies with Galactic reddening and imaging depth (based on DR5).

The bottom of the notebook also contains a simple regression of target density versus reddening, which we use in desitarget.mock.targets_truth.


In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt

In [2]:
import pandas as pd
from astropy.table import Table

In [3]:
from sklearn.mixture import GaussianMixture as GMM

In [4]:
import seaborn as sns
sns.set(style='white', font_scale=1.1, palette='deep')

In [5]:
%matplotlib inline

Read the heapix information table for DR5...

...take the log of the densities, and convert to a pandas dataframe.


In [6]:
def read_infotable():
    infofile = os.path.join(os.getenv('DESI_ROOT'), 'target', 'catalogs', 'hp-info-dr5-0.17.1.fits')
    info = Table.read(infofile)
    cols = ['{}DEPTH_{}_PERCENTILES'.format(prefix, band) for band in ('G', 'R', 'Z') for prefix in ('PSF', 'GAL')]
    for target in ('ELG', 'LRG', 'QSO', 'BGS'):
        with np.errstate(all='ignore'):
            info['DENSITY_{}'.format(target)] = np.log10(info['DENSITY_{}'.format(target)])
    info.remove_columns(cols)
    return info.to_pandas()

In [7]:
_info = read_infotable()
_info


Out[7]:
HPXID HPXAREA RA DEC EBV PSFDEPTH_G PSFDEPTH_R PSFDEPTH_Z GALDEPTH_G GALDEPTH_R ... NEXP_G NEXP_R NEXP_Z DENSITY_ALL DENSITY_ELG DENSITY_LRG DENSITY_QSO DENSITY_LYA DENSITY_BGS DENSITY_MWS
0 0 0.052456 45.000000 0.149208 0.088500 24.575624 24.286928 23.380188 24.343454 23.938946 ... 4 2 3 5261.567383 3.373628 2.771568 2.558960 0.0 3.337111 0.0
1 1 0.052456 45.175781 0.298417 0.106785 24.486000 24.268520 23.200266 24.241188 23.930630 ... 3 2 2 4899.357910 3.387416 2.742604 2.660417 0.0 3.243994 0.0
2 2 0.052456 44.824219 0.298417 0.083006 24.630730 24.333153 23.313522 24.393602 23.978317 ... 4 2 2 5985.985840 3.362991 2.581236 2.234449 0.0 3.513202 0.0
3 3 0.052456 45.000000 0.447628 0.090561 24.559509 24.311333 23.288239 24.317987 23.964497 ... 3 2 2 4823.103516 3.384010 2.771568 2.394149 0.0 3.234448 0.0
4 4 0.052456 45.351562 0.447628 0.106791 24.343983 24.076357 23.164227 24.101957 23.745398 ... 2 1 2 4251.193848 3.344664 2.711570 2.321599 0.0 3.155267 0.0
5 5 0.052456 45.527344 0.596842 0.097694 24.225916 23.954306 23.153233 23.985544 23.622919 ... 2 1 2 2573.592773 3.106281 2.581236 2.125304 0.0 2.970402 0.0
6 6 0.052456 45.175781 0.596842 0.095500 24.465931 24.220013 23.273710 24.224958 23.882454 ... 3 2 2 4098.684570 3.348392 2.641934 2.456297 0.0 3.119055 0.0
7 7 0.052456 45.351562 0.746060 0.093295 24.284094 24.047815 23.294277 24.030888 23.734314 ... 2 1 3 4308.384766 3.309590 2.798720 2.484326 0.0 3.194020 0.0
8 8 0.052456 44.648438 0.447628 0.081472 24.618589 24.345335 23.159859 24.374355 23.986643 ... 3 2 2 5032.803223 3.370111 2.359387 2.426334 0.0 3.373628 0.0
9 9 0.052456 44.824219 0.596842 0.083834 24.632193 24.349581 23.189117 24.392372 23.990971 ... 4 2 2 5490.331055 3.473331 2.660417 2.581236 0.0 3.284527 0.0
10 10 0.052456 44.472656 0.596842 0.077690 24.629478 24.435400 23.341959 24.382429 24.084297 ... 4 2 3 4422.766602 3.280206 2.394149 2.280206 0.0 3.355753 0.0
11 11 0.052456 44.648438 0.746060 0.077288 24.650406 24.372604 23.448380 24.410416 24.019161 ... 4 2 4 5128.121582 3.384010 2.581236 2.456297 0.0 3.348392 0.0
12 12 0.052456 45.000000 0.746060 0.088749 24.519896 24.255062 23.270868 24.282492 23.907808 ... 3 2 3 5242.503418 3.432494 2.757327 2.394149 0.0 3.248689 0.0
13 13 0.052456 45.175781 0.895283 0.089137 24.467789 24.168743 23.444002 24.229855 23.859684 ... 3 2 4 4823.103516 3.394149 2.678146 2.321599 0.0 3.253334 0.0
14 14 0.052456 44.824219 0.895283 0.075638 24.545727 24.260260 23.461494 24.310770 23.913776 ... 3 2 4 5795.349609 3.413745 2.952304 2.125304 0.0 3.352088 0.0
15 15 0.052456 45.000000 1.044512 0.079337 24.497442 24.182848 23.439434 24.258589 23.875919 ... 3 2 4 5757.222168 3.450468 2.871271 2.321599 0.0 3.333284 0.0
16 16 0.052456 45.703125 0.746060 0.084560 24.125685 23.887634 23.141653 23.904728 23.559713 ... 2 1 2 1315.391846 2.771568 2.058357 1.757327 0.0 2.785356 0.0
17 17 0.052456 45.878906 0.895283 0.085857 24.076340 23.700874 22.912132 23.849855 23.444963 ... 1 1 2 2211.383301 2.970402 2.558960 1.757327 0.0 2.961447 0.0
18 18 0.052456 45.527344 0.895283 0.096097 24.112364 24.024450 23.426197 23.903336 23.725445 ... 1 1 4 2916.738281 3.112715 2.558960 2.125304 0.0 3.065536 0.0
19 19 0.052456 45.703125 1.044512 0.095707 24.115780 23.738264 23.193213 23.886045 23.477760 ... 1 1 3 3507.711426 3.214705 2.359387 2.581236 0.0 3.149438 0.0
20 20 0.052456 46.054688 1.044512 0.085629 24.091810 23.714767 22.958872 23.873138 23.466526 ... 1 1 2 5642.840332 3.280206 2.456297 2.660417 0.0 3.484326 0.0
21 21 0.052456 46.230469 1.193748 0.078183 24.110388 23.720509 22.909740 23.890947 23.477365 ... 1 1 2 4441.830078 3.293043 2.602425 2.058357 0.0 3.309590 0.0
22 22 0.052456 45.878906 1.193748 0.093103 24.122772 23.732430 23.218002 23.895760 23.480539 ... 1 1 3 4518.084961 3.352088 2.660417 2.484326 0.0 3.199284 0.0
23 23 0.052456 46.054688 1.342993 0.080176 24.123001 23.723740 23.085318 23.899841 23.474148 ... 1 1 3 5471.267578 3.370111 2.942964 2.484326 0.0 3.325529 0.0
24 24 0.052456 45.351562 1.044512 0.098581 24.150642 23.769289 23.314091 23.921608 23.503830 ... 1 1 4 5566.585449 3.462050 2.785356 2.456297 0.0 3.340904 0.0
25 25 0.052456 45.527344 1.193748 0.102357 24.140881 23.750698 23.267275 23.909264 23.493347 ... 1 1 4 3812.729980 3.370111 2.602425 2.456297 0.0 3.004482 0.0
26 26 0.052456 45.175781 1.193748 0.082129 24.452795 24.104370 23.406866 24.216398 23.803457 ... 2 2 4 5814.413086 3.487032 2.695179 2.394149 0.0 3.340904 0.0
27 27 0.052456 45.351562 1.342993 0.082786 24.301718 23.978790 23.329519 24.074739 23.693916 ... 2 1 4 5051.867188 3.377116 2.785356 2.359387 0.0 3.284527 0.0
28 28 0.052456 45.703125 1.342993 0.105988 24.094769 23.709671 23.365061 23.868139 23.463650 ... 1 1 4 5242.503418 3.400780 2.859990 2.535479 0.0 3.262477 0.0
29 29 0.052456 45.878906 1.492246 0.094263 24.103170 23.724384 23.314394 23.875549 23.472290 ... 1 1 4 5395.012695 3.447523 2.871271 2.359387 0.0 3.266978 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
152819 786402 0.052456 314.121094 -1.492246 0.094376 24.018852 23.932585 23.210379 23.939960 23.832573 ... 3 2 4 4384.639160 3.137538 2.711570 2.359387 0.0 3.390796 0.0
152820 786403 0.052456 314.296875 -1.342993 0.102598 23.988201 23.959518 23.316519 23.910828 23.857100 ... 3 3 4 3793.666260 3.177833 2.558960 2.321599 0.0 3.243994 0.0
152821 786404 0.052456 314.648438 -1.342993 0.093823 23.903200 24.017292 23.544586 23.826229 23.906673 ... 2 3 4 4346.512207 3.166697 2.678146 2.280206 0.0 3.400780 0.0
152822 786405 0.052456 314.824219 -1.193748 0.087239 23.560488 23.814610 23.473936 23.484316 23.697956 ... 1 2 3 4155.875488 3.043634 2.535479 2.183296 0.0 3.450468 0.0
152823 786406 0.052456 314.472656 -1.193748 0.094280 23.957632 24.007113 23.469946 23.880722 23.898230 ... 3 3 4 3564.902344 3.099750 2.602425 1.979176 0.0 3.280206 0.0
152824 786407 0.052456 314.648438 -1.044512 0.082309 24.023838 24.056711 23.533920 23.934374 23.939646 ... 3 3 4 3259.884033 3.099750 2.711570 2.280206 0.0 3.172301 0.0
152825 786408 0.052456 313.945312 -1.342993 0.101245 23.993797 23.904636 23.094446 23.914877 23.805468 ... 3 3 4 4975.612305 3.155267 2.836509 2.859990 0.0 3.366566 0.0
152826 786409 0.052456 314.121094 -1.193748 0.094644 23.972733 23.903349 23.146341 23.894602 23.803793 ... 3 2 4 3564.902344 3.065536 2.484326 2.359387 0.0 3.301395 0.0
152827 786410 0.052456 313.769531 -1.193748 0.122275 23.963762 23.880398 23.072142 23.884373 23.781246 ... 3 2 3 4060.557373 2.970402 2.426334 2.510655 0.0 3.420085 0.0
152828 786411 0.052456 313.945312 -1.044512 0.099958 24.044668 23.971811 23.152668 23.966000 23.872461 ... 3 3 4 3965.239014 3.051058 2.727364 2.510655 0.0 3.329424 0.0
152829 786412 0.052456 314.296875 -1.044512 0.084715 23.996414 23.907940 23.040237 23.918442 23.808187 ... 3 3 3 3698.347900 3.072598 2.641934 2.602425 0.0 3.321599 0.0
152830 786413 0.052456 314.472656 -0.895283 0.085018 24.076115 23.951986 23.086266 23.979897 23.842068 ... 3 3 3 3736.475342 3.004482 2.581236 2.394149 0.0 3.333284 0.0
152831 786414 0.052456 314.121094 -0.895283 0.083892 24.040159 23.942726 23.130444 23.961306 23.843172 ... 3 3 4 4060.557373 3.043634 2.711570 2.359387 0.0 3.390796 0.0
152832 786415 0.052456 314.296875 -0.746060 0.085945 24.083170 23.956715 23.265049 23.987278 23.848045 ... 3 3 4 3469.584229 3.106281 2.602425 2.484326 0.0 3.243994 0.0
152833 786416 0.052456 315.000000 -1.044512 0.085760 23.121010 23.525394 23.435410 23.044760 23.388336 ... 1 1 2 3145.502197 2.695179 2.660417 1.757327 0.0 3.352088 0.0
152834 786417 0.052456 315.175781 -0.895283 0.080456 22.829222 23.379101 23.331753 22.554985 23.195923 ... 0 1 2 2783.292725 2.510655 2.785356 1.581236 0.0 3.271432 0.0
152835 786418 0.052456 314.824219 -0.895283 0.079681 23.907375 23.879704 23.463318 23.758825 23.726309 ... 2 3 3 2992.992920 2.923659 2.510655 2.234449 0.0 3.271432 0.0
152836 786419 0.052456 315.000000 -0.746060 0.072440 23.702925 23.666187 23.399021 23.432970 23.444141 ... 1 2 2 4213.066406 3.275841 2.558960 2.484326 0.0 3.257930 0.0
152837 786420 0.052456 315.351562 -0.746060 0.070227 23.728437 23.536110 23.122562 23.460102 23.295765 ... 1 1 2 3069.247559 3.112715 2.359387 2.321599 0.0 3.172301 0.0
152838 786421 0.052456 315.527344 -0.596842 0.073579 24.137239 23.938948 23.113367 23.913351 23.697338 ... 2 3 2 4232.130371 3.275841 2.602425 2.321599 0.0 3.288806 0.0
152839 786422 0.052456 315.175781 -0.596842 0.066928 23.828575 23.567972 23.166485 23.550447 23.315256 ... 1 1 2 3946.175293 3.166697 2.394149 2.321599 0.0 3.313630 0.0
152840 786423 0.052456 315.351562 -0.447628 0.065692 24.131567 23.851118 22.914812 23.895559 23.592812 ... 2 2 2 4289.321289 3.183296 2.798720 2.426334 0.0 3.321599 0.0
152841 786424 0.052456 314.648438 -0.746060 0.084029 24.230057 24.062405 23.303020 24.086124 23.919247 ... 3 3 5 3488.647949 3.099750 2.660417 2.280206 0.0 3.224689 0.0
152842 786425 0.052456 314.824219 -0.596842 0.076643 23.882961 23.576998 23.149717 23.616783 23.316378 ... 1 1 2 4327.448242 3.293043 2.456297 2.394149 0.0 3.284527 0.0
152843 786426 0.052456 314.472656 -0.596842 0.101459 24.240847 24.050329 23.277225 24.093578 23.907990 ... 3 3 4 3984.302734 3.072598 2.581236 2.678146 0.0 3.337111 0.0
152844 786427 0.052456 314.648438 -0.447628 0.099815 23.880224 23.556881 22.971525 23.623327 23.298447 ... 1 1 2 4708.721191 3.305512 2.660417 2.359387 0.0 3.348392 0.0
152845 786428 0.052456 315.000000 -0.447628 0.073731 23.851231 23.474743 22.992561 23.564070 23.185822 ... 1 1 2 4155.875488 3.248689 2.558960 2.535479 0.0 3.262477 0.0
152846 786429 0.052456 315.175781 -0.298417 0.072409 24.054617 23.747658 23.056339 23.796898 23.480560 ... 2 2 3 4060.557373 3.209625 2.622629 1.882266 0.0 3.313630 0.0
152847 786430 0.052456 314.824219 -0.298417 0.084263 23.884518 23.503691 22.919777 23.593447 23.214861 ... 1 1 2 4289.321289 3.149438 2.660417 2.359387 0.0 3.359387 0.0
152848 786431 0.052456 315.000000 -0.149208 0.083533 23.785103 23.464945 23.069107 23.499277 23.186174 ... 1 1 3 4537.148438 3.370111 2.484326 2.535479 0.0 3.224689 0.0

152849 rows × 21 columns

Select healpixels with three or more exposures in all three (grz) bands.

The caps on density are to toss out a handful of couple outliers.


In [8]:
minexp = 3

In [9]:
info = _info.loc[(_info['NEXP_G'] >= minexp) & (_info['NEXP_R'] >= minexp) & 
                 (_info['NEXP_Z'] >= minexp) & 
                 np.isfinite(_info['DENSITY_ELG']) & 
                 np.isfinite(_info['DENSITY_LRG']) & 
                 np.isfinite(_info['DENSITY_BGS']) & 
                 np.isfinite(_info['DENSITY_QSO'])
                 #(_info['DENSITY_ELG'] < 5000) & (_info['DENSITY_ELG'] > 0) & 
                 #(_info['DENSITY_LRG'] < 1500) & (_info['DENSITY_LRG'] > 0) & 
                 #(_info['DENSITY_BGS'] < 6000) & (_info['DENSITY_BGS'] > 0) & 
                 #(_info['DENSITY_QSO'] < 1000) & (_info['DENSITY_QSO'] > 0) 
                ]
print('There are {} / {} healpixels with >={} exposures in grz.'.format(len(info), len(_info), minexp))


There are 42741 / 152849 healpixels with >=3 exposures in grz.

Visualize the relationships between density, depth, and reddening for each target class.


In [10]:
def pairplot(target, only_dust=False):
    if only_dust:
        cols = ['DENSITY_{}'.format(target), 'EBV']
    else:
        cols = ['DENSITY_{}'.format(target), 'GALDEPTH_G', 'GALDEPTH_R', 'GALDEPTH_Z', 'EBV']
    print(info[cols].describe())
    sns.pairplot(info, vars=cols, diag_kind='hist', plot_kws={'s': 3}, size=2)

In [11]:
pairplot('ELG')


        DENSITY_ELG    GALDEPTH_G    GALDEPTH_R    GALDEPTH_Z           EBV
count  42741.000000  42741.000000  42741.000000  42741.000000  42741.000000
mean       3.374060     24.505133     24.110043     23.032774      0.041722
std        0.086989      0.308172      0.312294      0.251312      0.022888
min        1.581236     22.579695     21.880945     22.227142      0.008564
25%        3.325529     24.326187     23.896286     22.867973      0.026394
50%        3.380577     24.463923     24.125931     23.013962      0.035180
75%        3.429425     24.668438     24.283543     23.179966      0.049815
max        3.757327     25.882170     25.414646     24.174570      0.244828

In [12]:
pairplot('LRG')


        DENSITY_LRG    GALDEPTH_G    GALDEPTH_R    GALDEPTH_Z           EBV
count  42741.000000  42741.000000  42741.000000  42741.000000  42741.000000
mean       2.644160     24.505133     24.110043     23.032774      0.041722
std        0.146595      0.308172      0.312294      0.251312      0.022888
min        1.280206     22.579695     21.880945     22.227142      0.008564
25%        2.558960     24.326187     23.896286     22.867973      0.026394
50%        2.660417     24.463923     24.125931     23.013962      0.035180
75%        2.742604     24.668438     24.283543     23.179966      0.049815
max        3.256887     25.882170     25.414646     24.174570      0.244828

In [13]:
pairplot('QSO')


        DENSITY_QSO    GALDEPTH_G    GALDEPTH_R    GALDEPTH_Z           EBV
count  42741.000000  42741.000000  42741.000000  42741.000000  42741.000000
mean       2.370539     24.505133     24.110043     23.032774      0.041722
std        0.157054      0.308172      0.312294      0.251312      0.022888
min        1.280206     22.579695     21.880945     22.227142      0.008564
25%        2.280206     24.326187     23.896286     22.867973      0.026394
50%        2.394149     24.463923     24.125931     23.013962      0.035180
75%        2.484326     24.668438     24.283543     23.179966      0.049815
max        3.591960     25.882170     25.414646     24.174570      0.244828

In [14]:
pairplot('BGS')


        DENSITY_BGS    GALDEPTH_G    GALDEPTH_R    GALDEPTH_Z           EBV
count  42741.000000  42741.000000  42741.000000  42741.000000  42741.000000
mean       3.278477     24.505133     24.110043     23.032774      0.041722
std        0.114281      0.308172      0.312294      0.251312      0.022888
min        1.882266     22.579695     21.880945     22.227142      0.008564
25%        3.204485     24.326187     23.896286     22.867973      0.026394
50%        3.275841     24.463923     24.125931     23.013962      0.035180
75%        3.348392     24.668438     24.283543     23.179966      0.049815
max        4.036081     25.882170     25.414646     24.174570      0.244828

Determine the number of Gaussian components needed for each target class.

Regress just the QSOs against the dust reddening for now, to demonstrate the method. Unfortunately, a simple mixture of Gaussians does not appear to be a good representation of the data (i.e., the Bayesian information criterion, BIC, does not converge), so this approach requires more work.


In [15]:
ncomparray = np.arange(10, 50)

In [16]:
def qa_bic(ncomp, bic, title):
    fig, ax = plt.subplots(figsize=(8, 5))
    ax.plot(ncomp, bic, marker='s', ls='-')
    ax.set_xlabel('Number of Gaussian Components')
    ax.set_ylabel('Bayesian Information Criterion')
    ax.set_title(title)

In [17]:
def getbic(target, ncomp=[3]):
    cols = ['DENSITY_{}'.format(target), 'EBV']
    X = info[cols]
    bic = [GMM(n_components=nc).fit(X).bic(X) for nc in ncomp]
    return bic

In [18]:
%time bic_qso = getbic('QSO', ncomparray)


CPU times: user 5min 27s, sys: 11.4 s, total: 5min 39s
Wall time: 1min 25s

In [19]:
qa_bic(ncomparray, bic_qso, title='QSO')



In [20]:
qsocols = ['DENSITY_QSO', 'EBV']

In [21]:
mog = GMM(n_components=25).fit(info[qsocols])

QSO pairplot and data


In [22]:
print(info[qsocols].describe())
sns.pairplot(info[qsocols], plot_kws={'s': 10})


        DENSITY_QSO           EBV
count  42741.000000  42741.000000
mean       2.370539      0.041722
std        0.157054      0.022888
min        1.280206      0.008564
25%        2.280206      0.026394
50%        2.394149      0.035180
75%        2.484326      0.049815
max        3.591960      0.244828
Out[22]:
<seaborn.axisgrid.PairGrid at 0x117d16748>

QSO pairplot and model

Note how the model does not capture the correlation and spread in the data, especially at high reddening.


In [24]:
samp = pd.DataFrame(np.vstack(mog.sample(1000)[0]), columns=qsocols)

In [25]:
print(samp.describe())
sns.pairplot(samp, plot_kws={'s': 10})


       DENSITY_QSO          EBV
count  1000.000000  1000.000000
mean      2.365006     0.041066
std       0.150573     0.022535
min       1.579644    -0.014273
25%       2.279986     0.026766
50%       2.361255     0.038653
75%       2.460854     0.051553
max       2.836147     0.212188
Out[25]:
<seaborn.axisgrid.PairGrid at 0x115f54048>

Model (log) density vs E(B-V), for each target class.


In [26]:
def medxbin(x, y, binsize, minpts=20, xmin=None, xmax=None):
    """
    Compute the median (and other statistics) in fixed bins along the x-axis.
    """
    import numpy as np
    from scipy import ptp

    if xmin == None:
        xmin = x.min()
    if xmax == None:
        xmax = x.max()

    nbin = int(ptp(x) / binsize)
    bins = np.linspace(xmin, xmax, nbin)
    idx  = np.digitize(x, bins)

    stats = np.zeros(nbin,[('med', 'f4'), ('sig', 'f4'), 
                           ('q25', 'f4'), ('q75', 'f4') ])
    for kk in np.arange(nbin):
        npts = len(y[idx == kk])
        if npts > minpts:
            stats['med'][kk] = np.median(y[idx == kk])
            stats['sig'][kk] = np.std(y[idx == kk])
            stats['q25'][kk] = np.percentile(y[idx == kk], 25)
            stats['q75'][kk] = np.percentile(y[idx == kk], 75)

    good = np.nonzero(stats['med'])
    stats = stats[good]

    return bins[good], stats

In [27]:
def density_vs_ebv(info, target, order=2, doplot=True, ylim=None):

    xx, yy = info['EBV'].values, info['DENSITY_{}'.format(target)].values
    
    if ylim:
        good = (yy > ylim[0]) * (yy < ylim[1])
    else:
        good = np.ones(len(yy), dtype='bool')

    bins, stats = medxbin(xx[good], yy[good], 0.01, minpts=10, xmin=0)

    coeff = np.polyfit(bins, stats['med'], order)
    
    #from desiutil.funcfits import iter_fit
    #fit, mask = iter_fit(xx, yy, 'polynomial', order, sig_rej=10)
    #good = np.where(mask == 0)[0]
    #coeff = fit['coeff']
    #scatter = np.std(yy[good] - np.polyval(coeff, xx[good]))

    med = 10**np.mean(stats['med'])
    sig = np.log(10) * med * np.mean(stats['sig'])
    scatter = np.std(yy[good] - np.polyval(coeff, xx[good]))
    print('{}: median={:.2f}+/-{:.2f} deg2'.format(target, med, sig))
    print('  Slope, Intercept, scatter = {:.5f}, {:.3f}, {:.3f}'.format(coeff[0], coeff[1],
                                                                        scatter))
    
    if doplot:
        fig, ax = plt.subplots()
        ax.scatter(xx, yy, s=5, alpha=0.5)
        ax.plot(bins, stats['med'], lw=3, ls='-', color='orange')
        ax.plot(bins, stats['med'] + stats['sig'], lw=3, ls='--', color='orange')
        ax.plot(bins, stats['med'] - stats['sig'], lw=3, ls='--', color='orange')
        #ax.plot(bins, stats['q25'], lw=3, ls='--', color='orange')
        #ax.plot(bins, stats['q75'], lw=3, ls='--', color='orange')

        xaxis = np.linspace(xx.min(), xx.max(), 100)
        ax.plot(xaxis, np.polyval(coeff, xaxis), color='red', lw=2, alpha=0.5)
        #ax.plot(xaxis, np.polyval(coeff, xaxis) + scatter)
        #ax.plot(xaxis, np.polyval(coeff, xaxis) - scatter)
        
        ax.set_xlabel(r'$E(B-V)$')
        ax.set_ylabel(r'$\log_{{10}}\ (N\ {}\ /\ {{\rm deg}}^2)$'.format(target))
        
        if ylim:
            ax.set_ylim(ylim)
                                                                 
    return coeff

In [28]:
coeff = density_vs_ebv(info, 'QSO', order=1)


QSO: median=246.57+/-94.91 deg2
  Slope, Intercept, scatter = 0.02972, 2.389, 0.157

In [32]:
coeff = density_vs_ebv(info, 'ELG', order=1, ylim=(2.5, 4))


ELG: median=2092.90+/-441.48 deg2
  Slope, Intercept, scatter = -0.55792, 3.380, 0.081

In [30]:
coeff = density_vs_ebv(info, 'LRG', order=1, ylim=(1.7, 3.3))


LRG: median=457.16+/-161.89 deg2
  Slope, Intercept, scatter = 0.27216, 2.631, 0.145

In [31]:
coeff = density_vs_ebv(info, 'BGS', order=1, ylim=(2.5, 4))


BGS: median=1925.66+/-519.12 deg2
  Slope, Intercept, scatter = 0.33321, 3.249, 0.112