Chandra ACA Summer Project

Now With PyStan Bayesian Tools


In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from sklearn import linear_model
import astropy
import functions as f
import pystan

from numpy import genfromtxt
from astropy.time import Time

import warnings
warnings.filterwarnings('ignore')

Data import


In [2]:
#Setup
#Loading data from file
acq_data = np.load('data/acq_table.npy')

#Adding fields required for analysis
acq_data = f.add_column(acq_data, 'tstart_jyear' , np.zeros(len(acq_data)))
acq_data = f.add_column(acq_data, 'tstart_quarter' , np.zeros(len(acq_data)))
acq_data = f.add_column(acq_data, 'mag_floor' , np.zeros(len(acq_data)))
acq_data = f.add_column(acq_data, 'year' , np.zeros(len(acq_data)))
acq_data = f.add_column(acq_data, 'failure' , np.zeros(len(acq_data)))
acq_data['tstart_jyear'] = Time(acq_data['tstart'], format='cxcsec').jyear
acq_data['year'] = np.floor(acq_data.tstart_jyear)
acq_data['mag_floor'] = np.floor(acq_data['mag'])
acq_data['failure'][acq_data['obc_id']=="NOID"] = 1.

for acq in acq_data:
    acq.tstart_quarter = f.quarter_bin(acq.tstart_jyear)

New subset functions

Created to come up with similar analyses that I was performing in R.


In [3]:
def subset_range_tstart_jyear(dset, start, end):
    return dset[(dset['tstart_jyear']>=start) & (dset['tstart_jyear']<=end)]

In [4]:
def subset_range_warmpix(dset, wpf):
    return dset[(dset['warm_pix']>=(wpf - 0.01)) & (dset['warm_pix']<=(wpf + 0.01))]

Binomial proportion confidence interval estimation

We also need a function that estimates error bars for each bin that we'll eventually put together. Normal statistics textbooks would set up the interval as $\hat p \pm \sqrt{\frac{\hat p(1-\hat p)}{n}}$ where $n$ is the total number of trials and $\hat p$ in our case would be the sum of the failures over the total number of observations. This represents the results from the Central Limit Theorem. Unfortunately this does not work well for small $n$ or extreme probabilities, which is a majority of dataset. Thus in order to calculate our intervals, we'll use Wilson's method with a continuity correction which is better for small $n$ and extreme probabilities.

http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval


In [5]:
def wilsoninterval(failures, total_count):
    p = failures.astype(float)/total_count
    q = 1 - p
    z = 1.96
    n = total_count
    L = (2*n*p + z**2 - 1 - z*np.sqrt(z**2 - 2 - (1/n) + 4*p*(n*q + 1))) / (2*(n + z**2))
    U = (2*n*p + z**2 + 1 + z*np.sqrt(z**2 - 2 - (1/n) + 4*p*(n*q + 1))) / (2*(n + z**2)) 
    return U, L

Setting up the training subset of data

Below we set up the training data for the problem. Previously I had been using a subset by year and then using that year to test the next year. I've changed the code to be more variable. Right now a years worth of data takes approximately 5-10 mins to run for the model. The entire dataset takes approximately an hour to run.

Currently, I am using data from the last three quarters of 2013 and the first quarter of 2014 to estimate the newest observations in the dataset. Right now I don't believe that we need to use the entire dataset to create our estimates. We have a lot of data in the last year that is very "similar" to our present observations, while the data at the beginning of the mission is not "similar". I can show a few plots next week to describe what I'm talking about and give a little bit more of a visual justification.


In [6]:
#Subsetting up the training data
train_start = 2013.25
train_end = 2014.25

traindata = subset_range_tstart_jyear(acq_data, train_start, train_end)
train_success = f.smlset(traindata, 'obc_id', 'ID')
train_fail = f.smlset(traindata, 'obc_id', 'NOID')

Setting up the test subset of data

Below you can chance the Warm Pixel Fraction that we're testing. I'm subsetting the data this far up in the notebook so that we can get a visualization of what the training data looks like for the samples we're trying to estimate. Right now I'm using the newest data as out of sample test data to try and see how we'll perform in the future.


In [7]:
#Subsetting the Test Data
wp_test = 0.15

test_start = 2014.25
test_end = 2014.5

test_data = subset_range_tstart_jyear(acq_data, test_start, test_end)
wp_range = subset_range_warmpix(test_data, wp_test)

Visualization of the Training Data

Included are the following plots:

  • Histograms of the successes/failures for the training range.
  • Histogram of the observed warm pixel fractions for the training data.
    • Motivation: Are we using similar data points to what we eventually want to test?
  • Line chart of the proportion of failed stars binned by star magnitude.
    • Overall for the training period
  • Histograms of a subset of successes/failures for the training range for a given warm pixel fraction
    • For a given warm pixel fraction, I bin all successes and failures in a +/- 0.01 range
  • Histograms of the window of eventually tested warm pixel fraction
    • Motivation: have we observed a lot of events in warm pixel range we want to eventually test or are we extrapolating from the dataset.
  • Summary Table
    • Warm Pixel Fraction Range bins
    • Average Warm Pixel Fraction in from those bins
    • Columns: Mag Bin, Percent Failed, Total Observed, Failure Count
      • Shows specific failure statistics for the given warm pixel fraction range
      • Idea is to get a relative magnitude for the number of failures as opposed to the line chart which is proportions
  • Line chart of the proportion of failed stars binned by star magnitude for the given warm pixel fraction range

In [8]:
# Visualizing the Training Data
nbins = np.arange(6,11.75,0.25)

wp_train_range = subset_range_warmpix(traindata, wp_test)

F = plt.figure()
p1 = plt.hist(train_success['mag'], bins=nbins, color='b', alpha=0.5, label="Successes")
p2 = plt.hist(train_fail['mag'], bins=nbins, color='r', alpha=0.5, label="Failures")
lgd = plt.legend()
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {} (All Data)'.format(train_start, train_end))
plt.grid()
plt.show()

F = plt.figure()
p1 = plt.hist(traindata['warm_pix'])
F.set_size_inches(10,3)
plt.xlabel('Warm Pixel Fraction')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {} (All Data)'.format(train_start, train_end))
plt.grid()
plt.show()

failure_counts, bin_edges = np.histogram(train_fail['mag'], bins=nbins)
success_counts, bin_edges = np.histogram(train_success['mag'], bins=nbins)
total_counts, bin_edges = np.histogram(traindata['mag'], bins=nbins)
percent_fail = failure_counts.astype(float) / total_counts
pupper, plower = wilsoninterval(failure_counts, total_counts)

F = plt.figure()
p1, = plt.plot(bin_edges[1:], percent_fail)
p2, = plt.plot(bin_edges[1:], plower, linestyle='--', color='red')
p3, = plt.plot(bin_edges[1:], pupper, linestyle='--', color='red')
plt.legend([p1, p2],
           ["Observed Data", "95% Interval"], loc=2)
plt.ylim((0.,1.))
plt.xlim((6.,11.))
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Proportion Failed')
plt.title('Star Acquisitions for Time Range: {} - {} (All Data)'.format(train_start, train_end))
plt.grid()
plt.show()


failure_counts, bin_edges = np.histogram(wp_train_range[wp_train_range['obc_id']=="NOID"]['mag'], bins=nbins)
success_counts, bin_edges = np.histogram(wp_train_range[wp_train_range['obc_id']=="ID"]['mag'], bins=nbins)
total_counts, bin_edges = np.histogram(wp_train_range['mag'], bins=nbins)
percent_fail = failure_counts.astype(float) / total_counts
pupper, plower = wilsoninterval(failure_counts, total_counts)

F = plt.figure()
plt.hist(wp_train_range[wp_train_range['obc_id']=="ID"]['mag'], bins=nbins, color='b', alpha=0.5, label="Successes")
plt.hist(wp_train_range[wp_train_range['obc_id']=="NOID"]['mag'], bins=nbins, color='r', alpha=0.5, label="Failures")
lgd = plt.legend()
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {} (Warm Pix Subset)'.format(train_start, train_end))
plt.grid()
plt.show()

F = plt.figure()
plt.hist(wp_train_range['warm_pix'])
F.set_size_inches(10,3)
plt.xlabel('Warm Pixel Fraction')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {} '.format(train_start, train_end))
plt.grid()
plt.show()


print """Warm Pixel Fraction Range: {0:<3.2f} - {1:3.2f}""".format(wp_test - 0.01, wp_test + 0.01)
print """Average WPF              : {}\n""".format(np.mean(wp_train_range['warm_pix']))

print """{0:^15} | {1:^15} | {2:^15} | {3:^15}""".format("Mag Bin", "Percent Failed", "Total Observed", "Failure Count")
print """{}""".format("-"*70)
for i in np.arange(0, len(total_counts)):
    print """{0:<15} | {1:^15.3f} | {2:^15} | {3:^15}""".format("{0:<5.2f} - {1:>5.2f}".format(nbins[i+1] - 0.25, nbins[i+1]),
                                                    100*percent_fail[i], 
                                                    total_counts[i], 
                                                    failure_counts[i])  

F = plt.figure()
p1, = plt.plot(bin_edges[1:], percent_fail)
p2, = plt.plot(bin_edges[1:], plower, linestyle='--', color='red')
p3, = plt.plot(bin_edges[1:], pupper, linestyle='--', color='red')
plt.xlim((8,11))
plt.ylim((0,1))
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Proportion Failed')
plt.title('Star Acquisitions for Time Range: {} - {}'.format(train_start, train_end))
plt.grid()
plt.show()


Warm Pixel Fraction Range: 0.14 - 0.16
Average WPF              : 0.148038133042

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count 
----------------------------------------------------------------------
6.00  -  6.25   |      0.000      |       22        |        0       
6.25  -  6.50   |      0.000      |       21        |        0       
6.50  -  6.75   |      0.000      |       28        |        0       
6.75  -  7.00   |      0.000      |       29        |        0       
7.00  -  7.25   |      0.000      |       44        |        0       
7.25  -  7.50   |      5.882      |       51        |        3       
7.50  -  7.75   |      1.449      |       69        |        1       
7.75  -  8.00   |      1.802      |       111       |        2       
8.00  -  8.25   |      0.658      |       152       |        1       
8.25  -  8.50   |      0.000      |       174       |        0       
8.50  -  8.75   |      0.467      |       214       |        1       
8.75  -  9.00   |      1.045      |       287       |        3       
9.00  -  9.25   |      3.226      |       341       |       11       
9.25  -  9.50   |      5.391      |       371       |       20       
9.50  -  9.75   |      8.949      |       257       |       23       
9.75  - 10.00   |     16.923      |       260       |       44       
10.00 - 10.25   |     27.206      |       136       |       37       
10.25 - 10.50   |     55.814      |       43        |       24       
10.50 - 10.75   |     66.667      |       24        |       16       
10.75 - 11.00   |     75.000      |        4        |        3       
11.00 - 11.25   |       nan       |        0        |        0       
11.25 - 11.50   |       nan       |        0        |        0       

Visualization of the Training Data

Included are the following plots:

  • Histograms of the successes/failures for observations in the observed test range.
  • Histogram of the observed warm pixel fractions for the test data.
    • Motivation: Is this data similar to what we used to test?
  • Line chart of the proportion of failed stars binned by star magnitude.
    • Overall for the test period
  • Histograms of a subset of successes/failures for the test range for a given warm pixel fraction
    • For a given warm pixel fraction, I bin all successes and failures in a +/- 0.01 range
  • Histograms of the window of tested warm pixel fraction
    • Motivation: How many more are we observing now compared with then...
  • Summary Table
    • Warm Pixel Fraction Range bins
    • Average Warm Pixel Fraction in from those bins
    • Columns: Mag Bin, Percent Failed, Total Observed, Failure Count
      • Shows specific failure statistics for the given warm pixel fraction range
      • Idea is to get a relative magnitude for the number of failures as opposed to the line chart which is proportions
  • Line chart of the proportion of failed stars binned by star magnitude for the given warm pixel fraction range

In [9]:
# Test Data Visualization
nbins = np.arange(6,11.75,0.25)

test_data = subset_range_tstart_jyear(acq_data, test_start, test_end)
wp_range = subset_range_warmpix(test_data, wp_test)

failure_counts, bin_edges = np.histogram(test_data[test_data['obc_id']=="NOID"]['mag'], bins=nbins)
success_counts, bin_edges = np.histogram(test_data[test_data['obc_id']=="ID"]['mag'], bins=nbins)
total_counts, bin_edges = np.histogram(test_data['mag'], bins=nbins)
percent_fail = failure_counts.astype(float) / total_counts
pupper, plower = wilsoninterval(failure_counts, total_counts)

F = plt.figure()
plt.hist(test_data[test_data['obc_id']=="ID"]['mag'], bins=nbins, color='b', alpha=0.5, label="Successes")
plt.hist(test_data[test_data['obc_id']=="NOID"]['mag'], bins=nbins, color='r', alpha=0.5, label="Failures")
lgd = plt.legend()
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {} (All Data)'.format(test_start, test_end))
plt.grid()
plt.show()

F = plt.figure()
p1 = plt.hist(test_data['warm_pix'])
F.set_size_inches(10,3)
plt.xlabel('Warm Pixel Fraction')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {} (All Data)'.format(test_start, test_end))
plt.grid()
plt.show()

F = plt.figure()
p1, = plt.plot(bin_edges[1:], percent_fail)
p2, = plt.plot(bin_edges[1:], plower, linestyle='--', color='red')
p3, = plt.plot(bin_edges[1:], pupper, linestyle='--', color='red')
plt.xlim((8,11))
plt.ylim((0,1))
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Proportion Failed')
plt.title('Star Acquisitions for Time Range: {} - {}'.format(test_start, test_end))
plt.grid()
plt.show()

failure_counts, bin_edges = np.histogram(wp_range[wp_range['obc_id']=="NOID"]['mag'], bins=nbins)
success_counts, bin_edges = np.histogram(wp_range[wp_range['obc_id']=="ID"]['mag'], bins=nbins)
total_counts, bin_edges = np.histogram(wp_range['mag'], bins=nbins)
percent_fail = failure_counts.astype(float) / total_counts
pupper, plower = wilsoninterval(failure_counts, total_counts)

F = plt.figure()
plt.hist(wp_range[wp_range['obc_id']=="ID"]['mag'], bins=nbins, color='b', alpha=0.5, label="Successes")
plt.hist(wp_range[wp_range['obc_id']=="NOID"]['mag'], bins=nbins, color='r', alpha=0.5, label="Failures")
lgd = plt.legend()
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {}'.format(test_start, test_end))
plt.grid()
plt.show()

F = plt.figure()
plt.hist(wp_range['warm_pix'])
F.set_size_inches(10,3)
plt.xlabel('Warm Pixel Fraction')
plt.ylabel('Frequency')
plt.title('Star Acquisitions for Time Range: {} - {}'.format(test_start, test_end))
plt.grid()
plt.show()


print """Warm Pixel Fraction Range: {0:<3.2f} - {1:3.2f}""".format(wp_test - 0.01, wp_test + 0.01)
print """Average WPF              : {}\n""".format(np.mean(wp_range['warm_pix']))

print """{0:^15} | {1:^15} | {2:^15} | {3:^15}""".format("Mag Bin", "Percent Failed", "Total Observed", "Failure Count")
print """{}""".format("-"*70)
for i in np.arange(0, len(total_counts)):
    print """{0:<15} | {1:^15.3f} | {2:^15} | {3:^15}""".format("{0:<5.2f} - {1:>5.2f}".format(nbins[i+1] - 0.25, nbins[i+1]),
                                                    100*percent_fail[i], 
                                                    total_counts[i], 
                                                    failure_counts[i])  

F = plt.figure()
p1, = plt.plot(bin_edges[1:], percent_fail)
p2, = plt.plot(bin_edges[1:], plower, linestyle='--', color='red')
p3, = plt.plot(bin_edges[1:], pupper, linestyle='--', color='red')
plt.xlim((8,11))
plt.ylim((0,1))
F.set_size_inches(10,3)
plt.xlabel('Star Magnitude')
plt.ylabel('Proportion Failed')
plt.title('Star Acquisitions for Time Range: {} - {}'.format(test_start, test_end))
plt.grid()
plt.show()


Warm Pixel Fraction Range: 0.14 - 0.16
Average WPF              : 0.14961704778

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count 
----------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0       
6.25  -  6.50   |      0.000      |        6        |        0       
6.50  -  6.75   |      0.000      |       14        |        0       
6.75  -  7.00   |      0.000      |        9        |        0       
7.00  -  7.25   |      0.000      |        8        |        0       
7.25  -  7.50   |      0.000      |       21        |        0       
7.50  -  7.75   |      0.000      |       15        |        0       
7.75  -  8.00   |      0.000      |       41        |        0       
8.00  -  8.25   |      0.000      |       34        |        0       
8.25  -  8.50   |      1.961      |       51        |        1       
8.50  -  8.75   |      1.887      |       53        |        1       
8.75  -  9.00   |      1.220      |       82        |        1       
9.00  -  9.25   |      0.000      |       94        |        0       
9.25  -  9.50   |      4.505      |       111       |        5       
9.50  -  9.75   |      9.474      |       95        |        9       
9.75  - 10.00   |     18.182      |       66        |       12       
10.00 - 10.25   |     34.667      |       75        |       26       
10.25 - 10.50   |     61.111      |       18        |       11       
10.50 - 10.75   |     90.000      |       10        |        9       
10.75 - 11.00   |     50.000      |        2        |        1       
11.00 - 11.25   |       nan       |        0        |        0       
11.25 - 11.50   |       nan       |        0        |        0       

Welcome to PyStan

Below is the PyStan Model that we'll be using. Documentation can be found at the following locations:

To understand the "model" below, use the modeling documentation at the second link. The first link will allow you to understand how to access elements from the model fit as well as what's available.

Model Description

NOTE: All variables are centered for computational purposes. Without centering a years worth of data goes from completing in 5 mins to completing in 70 minutes. Do NOT uncenter and rerun, you won't be happy.

Response:

  • Star Acquisition Failure (1 corresponds to a failure).
    • Model Variable: y

Covariates:

  • Star Magnitude.
    • Model Variable: Mag
  • Warm Pixel Fraction
    • Model Variable: WarmPix

Model:

\begin{eqnarray} y_i &\sim& Bernoulli(\Phi(\eta_i)) \\ \eta_i &=& \alpha + \beta_1 \times Mag + \beta_2 \times Mag^2 + \beta_3 \times WarmPix \\ & & + \beta_4 \times Mag \times WarmPix + \beta_5 \times Mag^2 \times WarmPix \\ \beta &\sim& N(0, 5I) \end{eqnarray}

Coefficient Descriptions:

\begin{eqnarray} \alpha &:& \mbox{Intercept term (Included in beta vector)} \\ \beta_1 &:& \mbox{First Order Star Magnitude Coefficient} \\ \beta_2 &:& \mbox{Second Order Star Magnitude Coefficient} \\ \beta_3 &:& \mbox{Warm Pixel Fraction Coefficient} \\ \beta_4 &:& \mbox{First Order Interaction Coefficient} \\ \beta_5 &:& \mbox{Second Order Interaction Coefficient} \\ \end{eqnarray}

Model Runtime:

Factor of number of observations and number of iterations. Currently 3000 iterations are used on 4 chains. The Stan modeling language guidelines suggest NOT using 1 chain as reference. After burn in, and merging we have about 6000 samples to play with.

  • 1 years data approximately 7 mins
  • 14 years data approximately 1 hour

In [10]:
chandra_model = """
data {
  int N;             // number of observations
  int y[N];          // binary outcome for star acquisition event 'n'
  real Mag[N];       // predictive feature - star magnitude
  real WarmPix[N];   // predictive feature - warm pixel fraction
}
parameters {
  real alpha;  // intercept coefficient
  real beta1;  // Star Magnitude coefficient
  real beta2;
  real beta3;  // Warm Pixel Fraction coefficient
  real beta4;  // Magnitude and WP interaction coefficient
  real beta5;
}
model {
  alpha ~ normal(0,5);  // weakly informative
  beta1 ~ normal(0,5);  // Star Magnitude coefficient
  beta2 ~ normal(0,5);  // Star Magnitude Second Order
  beta3 ~ normal(0,5);  // Warm Pixel Fraction coefficient
  beta4 ~ normal(0,5);  // Magnitude and WP interaction coefficient
  beta5 ~ normal(0,5);  // M & WP Second Order Interaction
  
  for (n in 1:N)
    y[n] ~ bernoulli(Phi(alpha + 
                            beta1 * Mag[n] +
                            beta2 * Mag[n] * Mag[n] +
                            beta3 * WarmPix[n] + 
                            beta4 * Mag[n] * WarmPix[n] + 
                            beta5 * Mag[n] * Mag[n] * WarmPix[n])); //  
}
"""

chandra_data = {'N': int(len(traindata)),
                'y':traindata["failure"].astype(int),
                'Mag':traindata["mag"] - np.mean(traindata["mag"]),
                'WarmPix':traindata["warm_pix"] - np.mean(traindata["warm_pix"])
               }

chandra_fit = pystan.stan(model_code=chandra_model, data=chandra_data, iter=3000, chains=4)


INFO:pystan:
NOW ON CHAIN 0

INFO:pystan:
NOW ON CHAIN 1

INFO:pystan:
NOW ON CHAIN 2

INFO:pystan:
NOW ON CHAIN 3

Model Summary


In [11]:
chandra_fit


Out[11]:
Inference for Stan model: anon_model_c3826222a135b76d5333abf787a17622.
4 chains, each with iter=3000; warmup=1500; thin=1; 
post-warmup draws per chain=1500, total post-warmup draws=6000.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha  -2.16  8.7e-4   0.03  -2.23  -2.19  -2.16  -2.14   -2.1 1469.0    1.0
beta1   0.75  7.6e-4   0.03    0.7   0.73   0.75   0.77    0.8 1249.0    1.0
beta2   0.32  4.5e-4   0.02   0.29   0.31   0.32   0.34   0.36 1791.0    1.0
beta3  -2.89    0.03   1.67  -6.24  -4.01  -2.86  -1.74    0.3 4014.0    1.0
beta4    7.6    0.04   1.41   4.84   6.65   7.62   8.56  10.39 1587.0    1.0
beta5   3.16    0.02    1.0   1.23   2.49   3.14   3.84   5.15 2784.0    1.0
lp__   -2221    0.04    1.8  -2226  -2222  -2221  -2220  -2219 1905.0    1.0

Samples were drawn using NUTS(diag_e) at Tue Sep  2 13:00:16 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

PyStan Supplied Plots

The important thing to look at is the plots on the right showing the traces


In [12]:
chandra_fit.plot().set_size_inches(15,15)


Extracting the Simulated Draws for each Covariate


In [13]:
beta_draws = chandra_fit.extract()

Covariate Plots

Used to assess correlation between covariates


In [33]:
def plotdiagnostic(draws):
    n_coefs = len(draws) - 1
    nplots = n_coefs**2
    names = []
    for d in draws:
        names.append(d)
    
    F = plt.figure()
    spot = 1
    for n in names[:-1]:
        for i, d in enumerate(draws):
            if i == n_coefs:
                break
            else:
                if n == d:
                    plt.subplot(n_coefs,n_coefs,spot)
                    plt.hist(draws[d], alpha=0.5)
                    plt.grid()
                    plt.xticks(rotation=45)
                    if spot <= n_coefs:
                        plt.title(d)
                    spot += 1
                else:
                    plt.subplot(n_coefs,n_coefs,spot)
                    plt.plot(draws[n],draws[d],marker='.',linestyle='', alpha=0.3)
                    plt.grid()
                    plt.xticks(rotation=45)
                    if spot <= n_coefs:
                        plt.title(d)
                    spot += 1
    F.set_size_inches(15,15)
    plt.show()
plotdiagnostic(beta_draws)


Setting up Coefficient Vectors for use later.

Taking the posterior mean for each coefficient. Also computing variables to center the data.


In [15]:
a = np.mean(beta_draws['alpha'])
b1 = np.mean(beta_draws['beta1'])
b2 = np.mean(beta_draws['beta2'])
b3 = np.mean(beta_draws['beta3'])
b4 = np.mean(beta_draws['beta4'])
b5 = np.mean(beta_draws['beta5'])
betas = np.array([a, b1, b2, b3, b4, b5]) #, b3
m_center = np.mean(traindata['mag'])
wp_center = np.mean(traindata['warm_pix'])

In [16]:
def build_matrix_line(m, wp, m_center, wp_center):
    mag = m - m_center
    warm = wp - wp_center
    return np.array([1, mag, mag**2, warm, mag*warm, mag*mag*warm]) #

In [17]:
def build_matrix_fixedWP(wp, m_center, wp_center):
    mag = np.arange(5.5,13.25,0.25)
    Xline = [build_matrix_line(mag[0], wp, m_center, wp_center)]
    for i, m in enumerate(mag[1:],start=1):
        Xline = np.vstack((Xline, build_matrix_line(mag[i], wp, m_center, wp_center)))
    return Xline

Creating a function to plot observed values v. estmated values for a given warm pixel fraction


In [34]:
def compare_obs_2_calc(wp_test, betas):
    nbins = np.arange(6,11.75,0.25)
    
    test_data = subset_range_tstart_jyear(acq_data, test_start, test_end)
    wp_range = subset_range_warmpix(test_data, wp_test)
    
    test_matrix = build_matrix_fixedWP(wp_test, m_center, wp_center)

    y = stats.norm.cdf(np.dot(test_matrix, betas))
    x = np.arange(5.5,13.25,0.25)

    failure_counts, bin_edges = np.histogram(wp_range[wp_range['obc_id']=="NOID"]['mag'], bins=nbins)
    success_counts, bin_edges = np.histogram(wp_range[wp_range['obc_id']=="ID"]['mag'], bins=nbins)
    total_counts, bin_edges = np.histogram(wp_range['mag'], bins=nbins)
    percent_fail = failure_counts.astype(float) / total_counts
    pupper, plower = wilsoninterval(failure_counts, total_counts)
    
    exp_fails = []
    exp_sd = []
    for i in np.arange(0, len(total_counts)):
        mag = ((nbins[i+1] - 0.25) + nbins[i+1])/2.0
        pi = stats.norm.cdf(np.dot(build_matrix_line(mag, wp_test, m_center, wp_center), betas))
        exp_fail = total_counts[i] * pi
        std_dev = np.sqrt(total_counts[i] * pi * (1 - pi))    
        exp_fails.append(exp_fail)
        exp_sd.append(std_dev)

    upper = np.array(exp_fails) + 2 * np.array(exp_sd)
    lower = np.array(exp_fails) - 2 * np.array(exp_sd)    
    
    pltmax = np.max((np.max(upper), np.max(failure_counts)))
    
    F = plt.figure()
    plt.subplot(121)
    p1, = plt.plot(x,y)
    p2, = plt.plot(bin_edges[1:], percent_fail)
    p3, = plt.plot(bin_edges[1:], plower, linestyle='--', color='red')
    p4, = plt.plot(bin_edges[1:], pupper, linestyle='--', color='red')
    plt.grid()
    plt.legend([p2,p3,p1], ["Observed: WPF={}".format(wp_test), "Observed Interval", "Est. Probability Curve"], loc=2)
    plt.xlabel("Star Magnitude")
    plt.ylabel("Proportion Failed")
    plt.xlim((8,11))
    plt.subplot(122)
    p1, = plt.plot(bin_edges[1:], failure_counts)
    p2, = plt.plot(bin_edges[1:], exp_fails, color='r', alpha=0.5)
    p3, = plt.plot(bin_edges[1:], upper, color='r', linestyle="--", alpha=0.5)
    p4, = plt.plot(bin_edges[1:], lower, color='r', linestyle="--", alpha=0.5)
    plt.ylim((0,pltmax))
    plt.legend([p1,p2, p3], ["Observed Failures: WPF={}".format(wp_test), "Est. # Failures", "95% Interval"], loc=2)
    plt.xlabel("Star Magnitude")
    plt.ylabel("Number Failed")
    F.set_size_inches(15,5)
    plt.grid()
    plt.show()

In [35]:
compare_obs_2_calc(wp_test, betas)


Creating a function to show simulated draws with observed. As well as estimates

  • The figure plotted shows for each simulated draw the estimated probability curve.
  • Summary Table Description
    • Warm Pixel Fraction Range bins
    • Average Warm Pixel Fraction in from those bins
    • Columns: Mag Bin, Percent Failed, Total Observed, Failure Count, Expected Fails, 1-Std Dev
      • Shows specific failure statistics for the given warm pixel fraction range
      • Two new columns: Expected Fails, 1-Std Dev
        • Idea is when binned, data is binomial.
        • Binomial Expect Value: $E[X] = np$
        • Binomial Std Dev: $SD[X] = \sqrt{np(1-p)}$
      • These new columns give an idea of what we should have seen, given the probabilities that are estimated with errors quantified for those estimates. This is a rough idea for a goodness of fit measure.

In [20]:
def compare_show_errors(wp_test, beta_draws, betas):
    nbins = np.arange(6,11.75,0.25)
    
    test_data = subset_range_tstart_jyear(acq_data, test_start, test_end)
    wp_range = subset_range_warmpix(test_data, wp_test)    
    
    test_matrix = build_matrix_fixedWP(wp_test, m_center, wp_center)    
    
    failure_counts, bin_edges = np.histogram(wp_range[wp_range['obc_id']=="NOID"]['mag'], bins=nbins)
    success_counts, bin_edges = np.histogram(wp_range[wp_range['obc_id']=="ID"]['mag'], bins=nbins)
    total_counts, bin_edges = np.histogram(wp_range['mag'], bins=nbins)
    percent_fail = failure_counts.astype(float) / total_counts
    plower, pupper = wilsoninterval(failure_counts, total_counts)
    
    F = plt.figure()
    for n in np.arange(1,len(beta_draws['alpha']),1):
        betas = beta_draws['alpha'][n], beta_draws['beta1'][n], beta_draws['beta2'][n], beta_draws['beta3'][n], beta_draws['beta4'][n], beta_draws['beta5'][n]
        y = stats.norm.cdf(np.dot(test_matrix, betas))
        x = np.arange(5.5,13.25,0.25)
        plt.plot(x,y, alpha=0.5)
    p1, = plt.plot(bin_edges[1:], percent_fail, linewidth=5)
    p2, = plt.plot(bin_edges[1:], plower, linestyle='--', color='red', linewidth=3)
    p3, = plt.plot(bin_edges[1:], pupper, linestyle='--', color='red', linewidth=3)
    plt.legend([p1], ["Observed: WPF={}".format(wp_test)], loc=2)
    plt.xlabel("Star Magnitude")
    plt.ylabel("Proportion Failed")
    plt.xlim((8,11))
    F.set_size_inches(15,5)
    plt.grid()
    plt.show()

    print """Warm Pixel Fraction Range: {0:<3.2f} - {1:3.2f}""".format(wp_test - 0.01, wp_test + 0.01)
    print """Average WPF              : {}\n""".format(np.mean(wp_range['warm_pix']))

    print """{0:^15} | {1:^15} | {2:^15} | {3:^15} | {4:^15} | {5:^15}""".format("Mag Bin", 
                                                                                 "Percent Failed", 
                                                                                 "Total Observed", 
                                                                                 "Failure Count", 
                                                                                 "Expected Fails",
                                                                                 "1-Std Dev")
    print """{}""".format("-"*105)
    for i in np.arange(0, len(total_counts)):
        mag = ((nbins[i+1] - 0.25) + nbins[i+1])/2.0
        pi = stats.norm.cdf(np.dot(build_matrix_line(mag, wp_test, m_center, wp_center), betas))
        exp_fail = total_counts[i] * pi
        std_dev = np.sqrt(total_counts[i] * pi * (1 - pi))
        print """{0:<15} | {1:^15.3f} | {2:^15} | {3:^15} | {4:^15.1f} | {5:^15}""".format("{0:<5.2f} - {1:>5.2f}".format(nbins[i+1] - 0.25, nbins[i+1]),
                                                        100*percent_fail[i], 
                                                        total_counts[i], 
                                                        failure_counts[i],
                                                        exp_fail,
                                                        "+/- {:<4.1f}".format(std_dev))

In [21]:
compare_show_errors(wp_test, beta_draws, betas)


Warm Pixel Fraction Range: 0.14 - 0.16
Average WPF              : 0.14961704778

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.1       |    +/- 0.2     
6.25  -  6.50   |      0.000      |        6        |        0        |       0.1       |    +/- 0.3     
6.50  -  6.75   |      0.000      |       14        |        0        |       0.1       |    +/- 0.3     
6.75  -  7.00   |      0.000      |        9        |        0        |       0.0       |    +/- 0.2     
7.00  -  7.25   |      0.000      |        8        |        0        |       0.0       |    +/- 0.2     
7.25  -  7.50   |      0.000      |       21        |        0        |       0.1       |    +/- 0.2     
7.50  -  7.75   |      0.000      |       15        |        0        |       0.0       |    +/- 0.2     
7.75  -  8.00   |      0.000      |       41        |        0        |       0.1       |    +/- 0.3     
8.00  -  8.25   |      0.000      |       34        |        0        |       0.1       |    +/- 0.4     
8.25  -  8.50   |      1.961      |       51        |        1        |       0.3       |    +/- 0.5     
8.50  -  8.75   |      1.887      |       53        |        1        |       0.4       |    +/- 0.6     
8.75  -  9.00   |      1.220      |       82        |        1        |       1.1       |    +/- 1.0     
9.00  -  9.25   |      0.000      |       94        |        0        |       2.3       |    +/- 1.5     
9.25  -  9.50   |      4.505      |       111       |        5        |       5.3       |    +/- 2.2     
9.50  -  9.75   |      9.474      |       95        |        9        |       8.7       |    +/- 2.8     
9.75  - 10.00   |     18.182      |       66        |       12        |      11.4       |    +/- 3.1     
10.00 - 10.25   |     34.667      |       75        |       26        |      22.7       |    +/- 4.0     
10.25 - 10.50   |     61.111      |       18        |       11        |       8.7       |    +/- 2.1     
10.50 - 10.75   |     90.000      |       10        |        9        |       6.8       |    +/- 1.5     
10.75 - 11.00   |     50.000      |        2        |        1        |       1.7       |    +/- 0.5     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

Now what does it do if we give it a stars magnitude & warm pixel fraction


In [22]:
allbetas = np.vstack((beta_draws['alpha'], beta_draws['beta1'], beta_draws['beta2'], beta_draws['beta3'], beta_draws['beta4'], beta_draws['beta5']))

F = plt.figure()
plt.hist(stats.norm.cdf(np.dot(build_matrix_line(10.6,wp_test, m_center, wp_center), allbetas)), bins=25, normed=True)
plt.xlim((0,1))
plt.grid()
plt.show()


What does it look like for other values of Warm Pixel Fraction?


In [37]:
compare_obs_2_calc(0.11, betas)
compare_show_errors(0.11, beta_draws, betas)


Warm Pixel Fraction Range: 0.10 - 0.12
Average WPF              : 0.114558692465

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.1       |    +/- 0.3     
6.25  -  6.50   |      0.000      |        1        |        0        |       0.0       |    +/- 0.1     
6.50  -  6.75   |      0.000      |        8        |        0        |       0.1       |    +/- 0.3     
6.75  -  7.00   |      0.000      |        7        |        0        |       0.1       |    +/- 0.3     
7.00  -  7.25   |      0.000      |        6        |        0        |       0.0       |    +/- 0.2     
7.25  -  7.50   |      0.000      |        8        |        0        |       0.1       |    +/- 0.2     
7.50  -  7.75   |      0.000      |        6        |        0        |       0.0       |    +/- 0.2     
7.75  -  8.00   |      0.000      |       16        |        0        |       0.1       |    +/- 0.3     
8.00  -  8.25   |      0.000      |       23        |        0        |       0.2       |    +/- 0.4     
8.25  -  8.50   |      0.000      |       31        |        0        |       0.3       |    +/- 0.5     
8.50  -  8.75   |      0.000      |       20        |        0        |       0.2       |    +/- 0.5     
8.75  -  9.00   |      0.000      |       39        |        0        |       0.6       |    +/- 0.8     
9.00  -  9.25   |      0.000      |       29        |        0        |       0.7       |    +/- 0.8     
9.25  -  9.50   |      0.000      |       40        |        0        |       1.4       |    +/- 1.2     
9.50  -  9.75   |     23.529      |       34        |        8        |       2.0       |    +/- 1.4     
9.75  - 10.00   |      3.030      |       33        |        1        |       3.2       |    +/- 1.7     
10.00 - 10.25   |      6.667      |       15        |        1        |       2.4       |    +/- 1.4     
10.25 - 10.50   |     16.667      |        6        |        1        |       1.5       |    +/- 1.1     
10.50 - 10.75   |     25.000      |        4        |        1        |       1.5       |    +/- 1.0     
10.75 - 11.00   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [38]:
compare_obs_2_calc(0.12, betas)
compare_show_errors(0.12, beta_draws, betas)


Warm Pixel Fraction Range: 0.11 - 0.13
Average WPF              : 0.123955742305

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.1       |    +/- 0.3     
6.25  -  6.50   |      0.000      |        4        |        0        |       0.1       |    +/- 0.3     
6.50  -  6.75   |      0.000      |       10        |        0        |       0.1       |    +/- 0.4     
6.75  -  7.00   |      0.000      |       10        |        0        |       0.1       |    +/- 0.3     
7.00  -  7.25   |      0.000      |        6        |        0        |       0.0       |    +/- 0.2     
7.25  -  7.50   |      0.000      |       10        |        0        |       0.1       |    +/- 0.2     
7.50  -  7.75   |      0.000      |       13        |        0        |       0.1       |    +/- 0.3     
7.75  -  8.00   |      0.000      |       31        |        0        |       0.2       |    +/- 0.4     
8.00  -  8.25   |      0.000      |       47        |        0        |       0.3       |    +/- 0.5     
8.25  -  8.50   |      0.000      |       42        |        0        |       0.3       |    +/- 0.6     
8.50  -  8.75   |      0.000      |       46        |        0        |       0.5       |    +/- 0.7     
8.75  -  9.00   |      0.000      |       92        |        0        |       1.4       |    +/- 1.2     
9.00  -  9.25   |      2.899      |       69        |        2        |       1.6       |    +/- 1.3     
9.25  -  9.50   |      3.896      |       77        |        3        |       3.0       |    +/- 1.7     
9.50  -  9.75   |      8.065      |       62        |        5        |       4.1       |    +/- 2.0     
9.75  - 10.00   |      3.279      |       61        |        2        |       6.9       |    +/- 2.5     
10.00 - 10.25   |     14.035      |       57        |        8        |      10.8       |    +/- 3.0     
10.25 - 10.50   |     38.462      |       26        |       10        |       7.9       |    +/- 2.3     
10.50 - 10.75   |     28.571      |        7        |        2        |       3.2       |    +/- 1.3     
10.75 - 11.00   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [39]:
compare_obs_2_calc(0.13, betas)
compare_show_errors(0.13, beta_draws, betas)


Warm Pixel Fraction Range: 0.12 - 0.14
Average WPF              : 0.131210095991

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.1       |    +/- 0.2     
6.25  -  6.50   |      0.000      |        4        |        0        |       0.1       |    +/- 0.3     
6.50  -  6.75   |      0.000      |       11        |        0        |       0.1       |    +/- 0.3     
6.75  -  7.00   |      0.000      |       12        |        0        |       0.1       |    +/- 0.3     
7.00  -  7.25   |      0.000      |       13        |        0        |       0.1       |    +/- 0.3     
7.25  -  7.50   |      0.000      |       12        |        0        |       0.1       |    +/- 0.2     
7.50  -  7.75   |      0.000      |       22        |        0        |       0.1       |    +/- 0.3     
7.75  -  8.00   |      0.000      |       38        |        0        |       0.2       |    +/- 0.4     
8.00  -  8.25   |      0.000      |       56        |        0        |       0.3       |    +/- 0.5     
8.25  -  8.50   |      0.000      |       55        |        0        |       0.4       |    +/- 0.6     
8.50  -  8.75   |      0.000      |       72        |        0        |       0.7       |    +/- 0.8     
8.75  -  9.00   |      0.794      |       126       |        1        |       1.8       |    +/- 1.3     
9.00  -  9.25   |      1.639      |       122       |        2        |       2.9       |    +/- 1.7     
9.25  -  9.50   |      3.509      |       114       |        4        |       4.7       |    +/- 2.1     
9.50  -  9.75   |      2.128      |       94        |        2        |       6.9       |    +/- 2.5     
9.75  - 10.00   |      8.333      |       72        |        6        |       9.4       |    +/- 2.9     
10.00 - 10.25   |     18.447      |       103       |       19        |      23.0       |    +/- 4.2     
10.25 - 10.50   |     34.286      |       35        |       12        |      12.6       |    +/- 2.8     
10.50 - 10.75   |     20.000      |        5        |        1        |       2.7       |    +/- 1.1     
10.75 - 11.00   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [40]:
compare_obs_2_calc(0.14, betas)
compare_show_errors(0.14, beta_draws, betas)


Warm Pixel Fraction Range: 0.13 - 0.15
Average WPF              : 0.139901701288

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        1        |        0        |       0.0       |    +/- 0.2     
6.25  -  6.50   |      0.000      |        3        |        0        |       0.0       |    +/- 0.2     
6.50  -  6.75   |      0.000      |       15        |        0        |       0.1       |    +/- 0.4     
6.75  -  7.00   |      0.000      |       10        |        0        |       0.1       |    +/- 0.2     
7.00  -  7.25   |      0.000      |       11        |        0        |       0.0       |    +/- 0.2     
7.25  -  7.50   |      0.000      |       18        |        0        |       0.1       |    +/- 0.3     
7.50  -  7.75   |      0.000      |       18        |        0        |       0.1       |    +/- 0.2     
7.75  -  8.00   |      0.000      |       42        |        0        |       0.2       |    +/- 0.4     
8.00  -  8.25   |      0.000      |       45        |        0        |       0.2       |    +/- 0.4     
8.25  -  8.50   |      0.000      |       58        |        0        |       0.3       |    +/- 0.6     
8.50  -  8.75   |      0.000      |       67        |        0        |       0.6       |    +/- 0.8     
8.75  -  9.00   |      0.935      |       107       |        1        |       1.5       |    +/- 1.2     
9.00  -  9.25   |      0.000      |       114       |        0        |       2.8       |    +/- 1.6     
9.25  -  9.50   |      4.348      |       115       |        5        |       5.1       |    +/- 2.2     
9.50  -  9.75   |      6.140      |       114       |        7        |       9.4       |    +/- 2.9     
9.75  - 10.00   |     16.901      |       71        |       12        |      10.7       |    +/- 3.0     
10.00 - 10.25   |     27.184      |       103       |       28        |      27.0       |    +/- 4.5     
10.25 - 10.50   |     39.130      |       23        |        9        |       9.7       |    +/- 2.4     
10.50 - 10.75   |     62.500      |        8        |        5        |       4.9       |    +/- 1.4     
10.75 - 11.00   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [41]:
compare_obs_2_calc(0.15, betas)
compare_show_errors(0.15, beta_draws, betas)


Warm Pixel Fraction Range: 0.14 - 0.16
Average WPF              : 0.14961704778

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.1       |    +/- 0.2     
6.25  -  6.50   |      0.000      |        6        |        0        |       0.1       |    +/- 0.3     
6.50  -  6.75   |      0.000      |       14        |        0        |       0.1       |    +/- 0.3     
6.75  -  7.00   |      0.000      |        9        |        0        |       0.0       |    +/- 0.2     
7.00  -  7.25   |      0.000      |        8        |        0        |       0.0       |    +/- 0.2     
7.25  -  7.50   |      0.000      |       21        |        0        |       0.1       |    +/- 0.2     
7.50  -  7.75   |      0.000      |       15        |        0        |       0.0       |    +/- 0.2     
7.75  -  8.00   |      0.000      |       41        |        0        |       0.1       |    +/- 0.3     
8.00  -  8.25   |      0.000      |       34        |        0        |       0.1       |    +/- 0.4     
8.25  -  8.50   |      1.961      |       51        |        1        |       0.3       |    +/- 0.5     
8.50  -  8.75   |      1.887      |       53        |        1        |       0.4       |    +/- 0.6     
8.75  -  9.00   |      1.220      |       82        |        1        |       1.1       |    +/- 1.0     
9.00  -  9.25   |      0.000      |       94        |        0        |       2.3       |    +/- 1.5     
9.25  -  9.50   |      4.505      |       111       |        5        |       5.3       |    +/- 2.2     
9.50  -  9.75   |      9.474      |       95        |        9        |       8.7       |    +/- 2.8     
9.75  - 10.00   |     18.182      |       66        |       12        |      11.4       |    +/- 3.1     
10.00 - 10.25   |     34.667      |       75        |       26        |      22.7       |    +/- 4.0     
10.25 - 10.50   |     61.111      |       18        |       11        |       8.7       |    +/- 2.1     
10.50 - 10.75   |     90.000      |       10        |        9        |       6.8       |    +/- 1.5     
10.75 - 11.00   |     50.000      |        2        |        1        |       1.7       |    +/- 0.5     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [42]:
compare_obs_2_calc(0.16, betas)
compare_show_errors(0.16, beta_draws, betas)


Warm Pixel Fraction Range: 0.15 - 0.17
Average WPF              : 0.157103115504

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.0       |    +/- 0.2     
6.25  -  6.50   |     20.000      |        5        |        1        |       0.1       |    +/- 0.2     
6.50  -  6.75   |      0.000      |        8        |        0        |       0.1       |    +/- 0.2     
6.75  -  7.00   |      0.000      |        6        |        0        |       0.0       |    +/- 0.2     
7.00  -  7.25   |      0.000      |        7        |        0        |       0.0       |    +/- 0.1     
7.25  -  7.50   |      0.000      |       16        |        0        |       0.0       |    +/- 0.2     
7.50  -  7.75   |      0.000      |        9        |        0        |       0.0       |    +/- 0.1     
7.75  -  8.00   |      0.000      |       22        |        0        |       0.1       |    +/- 0.2     
8.00  -  8.25   |      0.000      |       22        |        0        |       0.1       |    +/- 0.3     
8.25  -  8.50   |      3.125      |       32        |        1        |       0.1       |    +/- 0.4     
8.50  -  8.75   |      2.778      |       36        |        1        |       0.3       |    +/- 0.5     
8.75  -  9.00   |      1.667      |       60        |        1        |       0.8       |    +/- 0.9     
9.00  -  9.25   |      0.000      |       61        |        0        |       1.5       |    +/- 1.2     
9.25  -  9.50   |      1.282      |       78        |        1        |       4.0       |    +/- 1.9     
9.50  -  9.75   |     11.765      |       51        |        6        |       5.2       |    +/- 2.2     
9.75  - 10.00   |     16.981      |       53        |        9        |      10.4       |    +/- 2.9     
10.00 - 10.25   |     34.211      |       38        |       13        |      13.2       |    +/- 2.9     
10.25 - 10.50   |     58.333      |       12        |        7        |       6.6       |    +/- 1.7     
10.50 - 10.75   |     100.000     |        6        |        6        |       4.5       |    +/- 1.1     
10.75 - 11.00   |     75.000      |        4        |        3        |       3.6       |    +/- 0.6     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [43]:
compare_obs_2_calc(0.17, betas)
compare_show_errors(0.17, beta_draws, betas)


Warm Pixel Fraction Range: 0.16 - 0.18
Average WPF              : 0.167272273587

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.0       |    +/- 0.2     
6.25  -  6.50   |     100.000     |        1        |        1        |       0.0       |    +/- 0.1     
6.50  -  6.75   |      0.000      |        2        |        0        |       0.0       |    +/- 0.1     
6.75  -  7.00   |      0.000      |        3        |        0        |       0.0       |    +/- 0.1     
7.00  -  7.25   |      0.000      |        4        |        0        |       0.0       |    +/- 0.1     
7.25  -  7.50   |      0.000      |        9        |        0        |       0.0       |    +/- 0.1     
7.50  -  7.75   |      0.000      |        4        |        0        |       0.0       |    +/- 0.1     
7.75  -  8.00   |      0.000      |       16        |        0        |       0.0       |    +/- 0.2     
8.00  -  8.25   |      0.000      |       15        |        0        |       0.0       |    +/- 0.2     
8.25  -  8.50   |      0.000      |       21        |        0        |       0.1       |    +/- 0.3     
8.50  -  8.75   |      0.000      |       21        |        0        |       0.1       |    +/- 0.4     
8.75  -  9.00   |      0.000      |       32        |        0        |       0.4       |    +/- 0.6     
9.00  -  9.25   |      0.000      |       25        |        0        |       0.6       |    +/- 0.8     
9.25  -  9.50   |      6.250      |       32        |        2        |       1.7       |    +/- 1.3     
9.50  -  9.75   |     11.429      |       35        |        4        |       4.0       |    +/- 1.9     
9.75  - 10.00   |     11.765      |       34        |        4        |       7.5       |    +/- 2.4     
10.00 - 10.25   |     27.778      |       18        |        5        |       7.1       |    +/- 2.1     
10.25 - 10.50   |     62.500      |        8        |        5        |       4.9       |    +/- 1.4     
10.50 - 10.75   |     75.000      |        4        |        3        |       3.2       |    +/- 0.8     
10.75 - 11.00   |     100.000     |        2        |        2        |       1.9       |    +/- 0.3     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

In [44]:
compare_obs_2_calc(0.18, betas)
compare_show_errors(0.18, beta_draws, betas)


Warm Pixel Fraction Range: 0.17 - 0.19
Average WPF              : 0.17046768592

    Mag Bin     | Percent Failed  | Total Observed  |  Failure Count  | Expected Fails  |    1-Std Dev   
---------------------------------------------------------------------------------------------------------
6.00  -  6.25   |      0.000      |        2        |        0        |       0.0       |    +/- 0.2     
6.25  -  6.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
6.50  -  6.75   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
6.75  -  7.00   |      0.000      |        1        |        0        |       0.0       |    +/- 0.1     
7.00  -  7.25   |      0.000      |        3        |        0        |       0.0       |    +/- 0.1     
7.25  -  7.50   |      0.000      |        4        |        0        |       0.0       |    +/- 0.1     
7.50  -  7.75   |      0.000      |        3        |        0        |       0.0       |    +/- 0.1     
7.75  -  8.00   |      0.000      |       11        |        0        |       0.0       |    +/- 0.1     
8.00  -  8.25   |      0.000      |        7        |        0        |       0.0       |    +/- 0.1     
8.25  -  8.50   |      0.000      |       12        |        0        |       0.0       |    +/- 0.2     
8.50  -  8.75   |      0.000      |       11        |        0        |       0.1       |    +/- 0.3     
8.75  -  9.00   |      0.000      |       11        |        0        |       0.1       |    +/- 0.4     
9.00  -  9.25   |      0.000      |       13        |        0        |       0.3       |    +/- 0.6     
9.25  -  9.50   |     15.385      |       13        |        2        |       0.8       |    +/- 0.8     
9.50  -  9.75   |     10.000      |       20        |        2        |       2.5       |    +/- 1.5     
9.75  - 10.00   |      0.000      |       13        |        0        |       3.2       |    +/- 1.6     
10.00 - 10.25   |     20.000      |        5        |        1        |       2.2       |    +/- 1.1     
10.25 - 10.50   |     60.000      |        5        |        3        |       3.3       |    +/- 1.1     
10.50 - 10.75   |     50.000      |        2        |        1        |       1.7       |    +/- 0.5     
10.75 - 11.00   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.00 - 11.25   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     
11.25 - 11.50   |       nan       |        0        |        0        |       0.0       |    +/- 0.0     

Let's gather the things we'll need later

Pickling the beta_draws and means


In [45]:
beta_results = {'names':np.array(['alpha', 'beta1', 'beta2', 'beta3', 'beta4', 'beta5']), 
                'means':betas,
                'draws':beta_draws,
                'betamatrix':np.vstack((beta_draws['alpha'],
                                        beta_draws['beta1'], 
                                        beta_draws['beta2'], 
                                        beta_draws['beta3'], 
                                        beta_draws['beta4'], 
                                        beta_draws['beta5'])),
                'centers':{'mag':m_center, 
                           'warmpix':wp_center}}

In [46]:
import pickle
with open('chandra_dump.pkl', 'wb') as f:
    pickle.dump(beta_results, f)