<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Building Histograms with Bayesian Priors
An Introduction to Bayesian Blocks
========
Version 0.1
By LM Walkowicz 2019 June 14
This notebook makes heavy use of Bayesian block implementations by Jeff Scargle, Jake VanderPlas, Jan Florjanczyk, and the Astropy team.
Before you begin, please download the dataset for this notebook.
Problem 1) Histograms Lie!
One of the most common and useful tools for data visualization can be incredibly misleading. Let's revisit how.
Problem 1a
First, let's make some histograms! Below, I provide some data; please make a histogram of it.
Step1:  Hey, nice histogram! 
But how do we know we have visualized all the relevant structure in our data? 
Play around with the binning and consider
Step2:  Problem 1b
What are some issues with histograms? 
Take a few min to discuss this with your partner
Solution 1b
write your solution here
Problem 1c
We have previously covered a few ways to make histograms better. What are some ways you could improve your histogram?
Take a few min to discuss this with your partner
Problem 1d
There are lots of ways to improve the previous histogram-- let's implement a KDE representation instead! As you have seen in previous sessions, we will borrow a bit of code from Jake VanderPlas to estimate the KDE. 
As a reminder, you have a number of choices of kernel in your KDE-- some we have used in the past
Step4:  Problem 1d
Which parameters most affected the shape of the final distribution?
What are some possible issues with using a KDE representation of the data?
Discuss with your partner
Solution 1d
Write your response here
Problem 2) Histograms Episode IV
Step5:  Problem 2a 
Let's visualize our data again, but this time we will use Bayesian Blocks.
Plot a standard histogram (as above), but now plot the Bayesian Blocks representation of the distribution over it.
Step6:  Problem 2b
How is the Bayesian Blocks representation different or similar? 
How might your choice of representation affect your scientific conclusions about your data? 
Take a few min to discuss this with your partner
If you are using histograms for analysis, you might infer physical meaning from the presence or absence of features in these distributions. As it happens, histograms of time-tagged event data are often used to characterize physical events in time domain astronomy, for example gamma ray bursts or stellar flares. 
Problem 3) Bayesian Blocks in the wild
Now we'll apply Bayesian Blocks to some real astronomical data, and explore how our visualization choices may affect our scientific conclusions.
First, let's get some data!
All data from NASA missions is hosted on the Mikulski Archive for Space Telescopes (aka MAST). As an aside, the M in MAST used to stand for "Multimission", but was changed to honor Sen. Barbara Mikulski (D-MD) for her tireless support of science. 
Some MAST data (mostly the original data products) can be directly accessed using astroquery (there's an extensive guide to interacting with MAST via astroquery here
Step7:  Problem 3c
Let's look at the distribution of the small planet radii in this table, which are given in units of Earth radii. Select the planets whose radii are Neptune-sized or smaller. 
Select the planet radii for all planets smaller than Neptune in the table, and visualize the distribution of planet radii using a standard histogram.
Step8:  Problem 3d
What features do you see in the histogram of planet radii? Which of these features are important? 
Discuss with your partner
Solution 3d
Write your answer here
Problem 3e
Now let's try visualizing these data using Bayesian Blocks. Please recreate the histogram you plotted above, and then plot the Bayesian Blocks version over it.
Step9:  Problem 3f
What features do you see in the histogram of planet radii? Which of these features are important? 
Discuss with your partner.
Hint
Step10:  Problem 3g
Please repeat the previous problem, but this time use the astropy implementation of Bayesian Blocks.
Step11:  Putting these results in context
Both standard histograms and KDEs can be useful for quickly visualizing data, and in some cases, getting an intuition for the underlying PDF of your data.
However, keep in mind that they both involve making parameter choices that are largely not motivated in any quantitative way. These choices can create wildly misleading representations of your data. 
In particular, your choices may lead you to make a physical interpretation that may or may not be correct (in our example, bear in mind that the observed distribution of exoplanetary radii informs models of planet formation).
Bayesian Blocks is more than just a variable-width histogram
While KDEs often do a better job of visualizing data than standard histograms do, they also create a loss of information. Philosophically speaking, what Bayesian Blocks do is posit that the "change points", also known as the bin edges, contain information that is interesting. When you apply a KDE, you are smoothing your data by creating an approximation, and that can mean you are losing potential insights by removing information. 
While Bayesian Blocks are useful as a replacement for histograms in general, their ability to identify change points makes them especially useful for time series analysis. 
Problem 4) Bayesian Blocks for Time Series Analysis
While Bayesian Blocks can be very useful as a simple replacement for histograms, one of its great strengths is in finding "change points" in time series data. Finding these change points can be useful for discovering interesting events in time series data you already have, and it can be used in real-time to detect changes that might trigger another action (for example, follow up observations for LSST). 
Let's take a look at a few examples of using Bayesian Blocks in the time series context. 
First and foremost, it's important to understand the different between various kinds of time series data. 
Event data come from photon counting instruments. In these data, the time series typically consists of photon arrival times, usually in a particular range of energies that the instrument is sensitive to. Event data are univariate, in that the time series is "how many photons at a given time", or "how many photons in a given chunk of time". 
Point measurements are measurements of a (typically) continuous source at a given moment in time, often with some uncertainty associated with the measurement. These data are multivariate, as your time series relates time, your measurement (e.g. flux, magnitude, etc) and its associated uncertainty to one another. 
Problem 4a
Let's look at some event data from BATSE, a high energy astrophysics experiment that flew on NASA's Compton Gamma-Ray Observatory. BATSE primarily studied gamma ray bursts (GRBs), capturing its detections in four energy channels
Step12:  Problem 4b
When you reach this point, you and your partner should pick a number between 1 and 4; your number is the channel whose data you will work with. 
Using the data for your channel, please visualize the photon events in both a standard histogram and using Bayesian Blocks.
Step13:  Problem 4c
Let's take a moment to reflect on the differences between these two representations of our data. 
Please discuss with your partner
Step14:  Cool, we have loaded in the FITS file. Let's look at what's in it
Step15:  We want the light curve, so let's check out what's in that part of the file!
Step16:  Problem 5b
Use a scatter plot to visualize the Kepler lightcurve. I strongly suggest you try displaying it at different scales by zooming in (or playing with the axis limits), so that you can get a better sense of the shape of the lightcurve.
Step17:  Problem 5c
These data consist of a variable background, with occasional bright points caused by stellar flares.
Brainstorm possible approaches to find the flare events in this data. Write down your ideas, and discuss their potential advantages, disadvantages, and any pitfalls that you think might arise. 
Discuss with your partner
Solution 5c
Write your notes here
There are lots of possible ways to approach this problem. In the literature, a very common traditional approach has been to fit the background variability while ignoring the outliers, then to subtract the background fit, and flag any point beyond some threshold value as belonging to a flare. More sophisticated approaches also exist, but they are often quite time consuming (and in many cases, detailed fits require a good starting estimate for the locations of flare events). 
Recall, however, that Bayesian Blocks is particularly effective at identifying change points in our data. Let's see if it can help us in this case!
Problem 5d
Use Bayesian Blocks to visualize the Kepler lightcurve. Note that you are now using data that consists of point measurements, rather than event data (as for the BATSE example).
Step18:  Concluding Remarks
As you can see, using Bayesian Blocks allowed us to represent the data, including both quiescent variability and flares, without having to smooth, clip, or otherwise alter the data. 
One potential drawback here is that the change points don't identify the flares themselves, or at least don't identify them as being different from the background variability-- the algorithm identifies change points, but does not know anything about what change points might be interesting to you in particular. 
Another potential drawback is that it is possible that the Bayesian Block representation may not catch all events, or at least may not, on its own, provide an unambiguous sign that a subtle signal of interest is in the data.
In these cases, it is sometimes instructive to use a hybrid approach, where one combines the bins determined from a traditional histogram with the Bayesian Blocks change points. Alternatively, if one has a good model for the background, one can compared the blocks representation of a background-only simulated data set with one containing both background and model signal. 
Further interesting examples (in the area of high energy physics and astrophysics) are provided by this paper
<ASSISTANT_TASK:>
Python Code:
# execute this cell
np.random.seed(0)
x = np.concatenate([stats.cauchy(-5, 1.8).rvs(500),
                    stats.cauchy(-4, 0.8).rvs(2000),
                    stats.cauchy(-1, 0.3).rvs(500),
                    stats.cauchy(2, 0.8).rvs(1000),
                    stats.cauchy(4, 1.5).rvs(500)])
# truncate values to a reasonable range
x = x[(x > -15) & (x < 15)]
# complete
# plt.hist( 
Explanation: Building Histograms with Bayesian Priors
An Introduction to Bayesian Blocks
========
Version 0.1
By LM Walkowicz 2019 June 14
This notebook makes heavy use of Bayesian block implementations by Jeff Scargle, Jake VanderPlas, Jan Florjanczyk, and the Astropy team.
Before you begin, please download the dataset for this notebook.
Problem 1) Histograms Lie!
One of the most common and useful tools for data visualization can be incredibly misleading. Let's revisit how.
Problem 1a
First, let's make some histograms! Below, I provide some data; please make a histogram of it.
End of explanation
# complete
# plt.hist 
Explanation: Hey, nice histogram! 
But how do we know we have visualized all the relevant structure in our data? 
Play around with the binning and consider: 
What features do you see in this data? Which of these features are important?
End of explanation
# execute this cell
from sklearn.neighbors import KernelDensity
def kde_sklearn(data, grid, bandwidth = 1.0, **kwargs):
    kde_skl = KernelDensity(bandwidth = bandwidth, **kwargs)
    kde_skl.fit(data[:, np.newaxis])
    log_pdf = kde_skl.score_samples(grid[:, np.newaxis]) # sklearn returns log(density)
    return np.exp(log_pdf)
# complete
# plt.hist( 
# grid = 
# PDF = 
# plt.plot(
Explanation: Problem 1b
What are some issues with histograms? 
Take a few min to discuss this with your partner
Solution 1b
write your solution here
Problem 1c
We have previously covered a few ways to make histograms better. What are some ways you could improve your histogram?
Take a few min to discuss this with your partner
Problem 1d
There are lots of ways to improve the previous histogram-- let's implement a KDE representation instead! As you have seen in previous sessions, we will borrow a bit of code from Jake VanderPlas to estimate the KDE. 
As a reminder, you have a number of choices of kernel in your KDE-- some we have used in the past: tophat, Epanechnikov, Gaussian. Please plot your original histogram, and then overplot a few example KDEs on top of it.
End of explanation
# execute this cell
def bayesian_blocks(t):
    Bayesian Blocks Implementation
    By Jake Vanderplas.  License: BSD
    Based on algorithm outlined in http://adsabs.harvard.edu/abs/2012arXiv1207.5578S
    Parameters
    ----------
    t : ndarray, length N
        data to be histogrammed
    Returns
    -------
    bins : ndarray
        array containing the (N+1) bin edges
    Notes
    -----
    This is an incomplete implementation: it may fail for some
    datasets.  Alternate fitness functions and prior forms can
    be found in the paper listed above.
    
    # copy and sort the array
    t = np.sort(t)
    N = t.size
    # create length-(N + 1) array of cell edges
    edges = np.concatenate([t[:1],
                            0.5 * (t[1:] + t[:-1]),
                            t[-1:]])
    block_length = t[-1] - edges
    # arrays needed for the iteration
    nn_vec = np.ones(N)
    best = np.zeros(N, dtype=float)
    last = np.zeros(N, dtype=int)
    #-----------------------------------------------------------------
    # Start with first data cell; add one cell at each iteration
    #-----------------------------------------------------------------
    for K in range(N):
        # Compute the width and count of the final bin for all possible
        # locations of the K^th changepoint
        width = block_length[:K + 1] - block_length[K + 1]
        count_vec = np.cumsum(nn_vec[:K + 1][::-1])[::-1]
        # evaluate fitness function for these possibilities
        fit_vec = count_vec * (np.log(count_vec) - np.log(width))
        fit_vec -= 4  # 4 comes from the prior on the number of changepoints
        fit_vec[1:] += best[:K]
        # find the max of the fitness: this is the K^th changepoint
        i_max = np.argmax(fit_vec)
        last[K] = i_max
        best[K] = fit_vec[i_max]
    
    #-----------------------------------------------------------------
    # Recover changepoints by iteratively peeling off the last block
    #-----------------------------------------------------------------
    change_points =  np.zeros(N, dtype=int)
    i_cp = N
    ind = N
    while True:
        i_cp -= 1
        change_points[i_cp] = ind
        if ind == 0:
            break
        ind = last[ind - 1]
    change_points = change_points[i_cp:]
    return edges[change_points]
Explanation: Problem 1d
Which parameters most affected the shape of the final distribution?
What are some possible issues with using a KDE representation of the data?
Discuss with your partner
Solution 1d
Write your response here
Problem 2) Histograms Episode IV: A New Hope
How can we create representations of our data that are robust against the known issues with histograms and KDEs?
Introducing: Bayesian Blocks
We want to represent our data in the most general possible way, a method that
* avoids assumptions about smoothness or shape of the signal (which might place limitations on scales and resolution)
* is nonparametric (doesn't fit some model)
* finds and characterizes local structure in our time series (in contrast to periodicities)
(continued)
handles arbitrary sampling (i.e. doesn't require evenly spaced samples, doesn't care about sparse samples)
is as hands-off as possible-- user interventions should be minimal or non-existent 
is applicable to multivariate data
can both analyze data after they are collected, and in real time 
Bayesian Blocks works by creating a super-simple representation of the data, essentially a piecewise fit that segments our time series. 
In the implementations we will use today, the model is a piecewise linear fit in time across each individual bin, or "block". 
one modeling the signal as linear in time across the block:
$$x(t) = λ(1 + a(t − t_{fid}))$$ 
where $\lambda$ is the signal strength at the fiducial time $t_{fid}$, and the coefficient $a$ determines the rate of change over the block. 
Scargle et al. (2012) point out that using a linear fit is good because it makes calculating the fit really easy, but you could potentially use something more complicated (they provide some details for using an exponential model, $x(t) = λe^{a(t−t_{fid}})$, in their Appendix C.
The Fitness Function
The insight in Bayesian Blocks is that you can use a Bayesian likelihood framework to compute a "fitness function" that depends only on the number and size of the blocks. 
In every block, you are trying to maximize some goodness-of-fit measure for data in that individual block. This fit depends only on the data contained in its block, and is independent of all other data.
The optimal segmentation of the time series, then, is the segmentation that maximizes fitness-- the total goodness-of-fit over all the blocks (so for example, you could use the sum over all blocks of whatever your quantitative expression is for the goodness of fit in individual blocks). 
The entire time series is then represented by a series of segments (blocks) characterized by very few parameters:
$N_{cp}$: the number of change-points
$t_{k}^{cp}$: the change-point starting block k
$X_k$: the signal amplitude in block k
for k = 1, 2, ... $N_{cp}$
When using Bayesian Blocks (particularly on time series data) we often speak of "change points" rather than segmentation, as the block edges essentially tell us the discrete times at which a signal’s statistical properties change discontinuously, though the segments themselves are constant between these points. 
You looking at KDEs right now:
In some cases (such as some of the examples below), the Bayesian Block representation may look kind of clunky. HOWEVER: remember that histograms and KDEs may sometimes look nicer, but can be really misleading! If you want to derive physical insight from these representations of your data, Bayesian Blocks can provide a means of deriving physically interesting quantities (for example, better estimates of event locations, lags, amplitudes, widths, rise and decay times, etc). 
On top of that, you can do all of the above without losing or hiding information via smoothing or other model assumptions. 
HOW MANY BLOCKS, THO?!
We began this lesson by bemoaning that histograms force us to choose a number of bins, and that KDEs require us to choose a bandwidth. Furthermore, one of the requirements we had for a better way forward was that the user interaction be minimal or non-existent. What to do?
Bayesian Blocks works by defining a prior distribution for the number of blocks, such that a single parameter controls the steepness of this prior (in other words, the relative probability for smaller or larger numbers of blocks. 
Once this prior is defined, the size, number, and locations of the blocks are determined solely and uniquely by the data. 
So, what does the prior look like?
In most cases, $N_{blocks}$ << N (you are, after all, still binning your data-- if $N_{blocks}$ was close to N, you wouldn't really be doing much). Scargle et al. (2012) adopts a geometric prior (Coram 2002), which assigns smaller probability to a large number of blocks:
$$P(N_{blocks}) = P_{0}\gamma N_{blocks}$$
for $0 ≤ N_{blocks} ≤ N$, and zero otherwise since $N_{blocks}$ cannot be negative or larger than the number of data cells. 
Substituting in the normalization constant $P_{0}$ gives
$$P(N_{blocks}) = \frac{1−\gamma}{1-\gamma^{N+1}}\gamma^{N_{blocks}}$$
<sub>Essentially, this prior says that finding k + 1 blocks is less likely than finding k blocks by the constant factor $\gamma$. Scargle (2012) also provides a nice intuitive way of thinking about $\gamma$: $\gamma$ is adjusting the amount of structure in the resulting representation.</sub>
The Magic Part
At this point, you may be wondering about how the algorithm is capable of finding an optimal number of blocks. As Scargle et al (2012) admits
the number of possible partitions (i.e. the number of ways N cells can be arranged in blocks) is $2^N$. This number is exponentially large, rendering an explicit exhaustive search of partition space utterly impossible for all but very small N.
In his blog post on Bayesian Blocks, Jake VdP compares the algorithm's use of dynamic programming to mathematical induction. For example: how could you prove that  
$$1 + 2 + \cdots + n = \frac{n(n+1)}{2}$$
is true for all positive integers $n$? An inductive proof of this formula proceeds in the following fashion:
Base Case: We can easily show that the formula holds for $n = 1$.
Inductive Step: For some value $k$, assume that $1 + 2 + \cdots + k = \frac{k(k+1)}{2}$ holds. 
Adding $(k + 1)$ to each side and rearranging the result yields 
$$1 + 2 + \cdots + k + (k + 1) = \frac{(k + 1)(k + 2)}{2}$$
Looking closely at this, we see that we have shown the following: if our formula is true for $k$, then it must be true for $k + 1$.
By 1 and 2, we can show that the formula is true for any positive integer $n$, simply by starting at $n=1$ and repeating the inductive step $n - 1$ times.
In the Bayesian Blocks algorithm, one can find the optimal binning for a single data point; so by analogy with our example above (full details are given in the Appendix of Scargle et al. 2012), if you can find the optimal binning for $k$ points, it's a short step to the optimal binning for $k + 1$ points. 
So, rather than performing an exhaustive search of all possible bins, which would scale as $2^N$, the time to find the optimal binning instead scales as $N^2$.
Playing with Blocks
We will begin playing with (Bayesian) blocks with a simple implementation, outlined by Jake VanderPlas in this blog: https://jakevdp.github.io/blog/2012/09/12/dynamic-programming-in-python/
End of explanation
# complete 
plt.hist(
Explanation: Problem 2a 
Let's visualize our data again, but this time we will use Bayesian Blocks.
Plot a standard histogram (as above), but now plot the Bayesian Blocks representation of the distribution over it.
End of explanation
# complete
Explanation: Problem 2b
How is the Bayesian Blocks representation different or similar? 
How might your choice of representation affect your scientific conclusions about your data? 
Take a few min to discuss this with your partner
If you are using histograms for analysis, you might infer physical meaning from the presence or absence of features in these distributions. As it happens, histograms of time-tagged event data are often used to characterize physical events in time domain astronomy, for example gamma ray bursts or stellar flares. 
Problem 3) Bayesian Blocks in the wild
Now we'll apply Bayesian Blocks to some real astronomical data, and explore how our visualization choices may affect our scientific conclusions.
First, let's get some data!
All data from NASA missions is hosted on the Mikulski Archive for Space Telescopes (aka MAST). As an aside, the M in MAST used to stand for "Multimission", but was changed to honor Sen. Barbara Mikulski (D-MD) for her tireless support of science. 
Some MAST data (mostly the original data products) can be directly accessed using astroquery (there's an extensive guide to interacting with MAST via astroquery here: https://astroquery.readthedocs.io/en/latest/mast/mast.html). 
In addition, MAST also hosts what are called "Higher Level Science Products", or HLSPs, which are data derived by science teams in the course of doing their analyses. You can see a full list of HLSPs here: https://archive.stsci.edu/hlsp/hlsp-table
These data tend to be more heterogeneous, and so are not currently accessible through astroquery (for the most part). They will be added in the future. But never fear! You can also submit SQL queries via MAST's CasJobs interface.
Go to the MAST CasJobs http://mastweb.stsci.edu/mcasjobs/home.aspx
If I have properly remembered to tell you to create a MAST CasJobs login, you can login now (and if not, just go ahead and sign up now, it's fast)!
We will be working with the table of new planet radii by Berger et al. (2019). 
If you like, you can check out the paper here! https://arxiv.org/pdf/1805.00231.pdf
From the "Query" tab, select "HLSP_KG_RADII" from the Context drop-down menu. 
You can then enter your query. In this example, we are doing a simple query to get all the KG-RADII radii and fluxes from the exoplanets catalog, which you could use to reproduce the first figure, above. For short queries that can execute in less than 60 seconds, you can hit the "Quick" button and the results of your query will be displayed below, where you can export them as needed. For longer queries like this one, you can select into an output table (otherwise a default like MyDB.MyTable will be used), hit the "Submit" button, and when finished your output table will be available in the MyDB tab.
Problem 3a
Write a SQL query to fetch this table from MAST using CasJobs.
Your possible variables are KIC_ID, KOI_ID, Planet_Radius, Planet_Radius_err_upper, Planet_Radius_err_lower, Incident_Flux, Incident_Flux_err_upper, Incident_Flux_err_lower, AO_Binary_Flag
For very short queries you can use "Quick" for your query; this table is large enough that you should use "Submit".
Hint: You will want to SELECT some stuff FROM a table called exoplanet_parameters
Solution 3a
Write your SQL query here
Once your query has completed, you will go to the MyDB tab to see the tables you have generated (in the menu at left). From here, you can click on a table, and select Download. I would recommend downloading your file as a CSV (comma-separated value) file, as CSV are simple, and can easily read into python via a variety of methods.
Problem 3b
Time to read in the data! There are several ways of importing a csv into python... choose your favorite and load in the table you downloaded.
End of explanation
# complete
Explanation: Problem 3c
Let's look at the distribution of the small planet radii in this table, which are given in units of Earth radii. Select the planets whose radii are Neptune-sized or smaller. 
Select the planet radii for all planets smaller than Neptune in the table, and visualize the distribution of planet radii using a standard histogram.
End of explanation
# complete
Explanation: Problem 3d
What features do you see in the histogram of planet radii? Which of these features are important? 
Discuss with your partner
Solution 3d
Write your answer here
Problem 3e
Now let's try visualizing these data using Bayesian Blocks. Please recreate the histogram you plotted above, and then plot the Bayesian Blocks version over it.
End of explanation
import astropy.stats.bayesian_blocks as bb
Explanation: Problem 3f
What features do you see in the histogram of planet radii? Which of these features are important? 
Discuss with your partner.
Hint: maybe you should look at some of the comments in the implementation of Bayesian Blocks we are using
Solution 3f
Write your answer here
OK, so in this case, the Bayesian Blocks representation of the data looks fairly different. A couple things might stand out to you:
* There are large spikes in the Bayesian Blocks representation that are not present in the standard histogram
* If we're just looking at the Bayesian Block representation, it's not totally clear whether one should believe that there are two peaks in the distribution. 
HMMMMMmmmm....
Wait! Jake VDP told us to watch out for this implementation. Maybe we should use something a little more official instead...
GOOD NEWS~! There is a Bayesian Blocks implementation included in astropy. Let's try that. 
http://docs.astropy.org/en/stable/api/astropy.stats.bayesian_blocks.html
Important note
There is a known issue in the astropy implementation of Bayesian Blocks; see: https://github.com/astropy/astropy/issues/8317
It is possible this issue will be fixed in a future release, but in the event this problem arises for you, you will need to edit bayesian_blocks.py to include the following else statement (see issue link above for exact edit):
if self.ncp_prior is None:
            ncp_prior = self.compute_ncp_prior(N)
 else:
            ncp_prior = self.ncp_prior
End of explanation
# complete
Explanation: Problem 3g
Please repeat the previous problem, but this time use the astropy implementation of Bayesian Blocks.
End of explanation
# complete
Explanation: Putting these results in context
Both standard histograms and KDEs can be useful for quickly visualizing data, and in some cases, getting an intuition for the underlying PDF of your data.
However, keep in mind that they both involve making parameter choices that are largely not motivated in any quantitative way. These choices can create wildly misleading representations of your data. 
In particular, your choices may lead you to make a physical interpretation that may or may not be correct (in our example, bear in mind that the observed distribution of exoplanetary radii informs models of planet formation).
Bayesian Blocks is more than just a variable-width histogram
While KDEs often do a better job of visualizing data than standard histograms do, they also create a loss of information. Philosophically speaking, what Bayesian Blocks do is posit that the "change points", also known as the bin edges, contain information that is interesting. When you apply a KDE, you are smoothing your data by creating an approximation, and that can mean you are losing potential insights by removing information. 
While Bayesian Blocks are useful as a replacement for histograms in general, their ability to identify change points makes them especially useful for time series analysis. 
Problem 4) Bayesian Blocks for Time Series Analysis
While Bayesian Blocks can be very useful as a simple replacement for histograms, one of its great strengths is in finding "change points" in time series data. Finding these change points can be useful for discovering interesting events in time series data you already have, and it can be used in real-time to detect changes that might trigger another action (for example, follow up observations for LSST). 
Let's take a look at a few examples of using Bayesian Blocks in the time series context. 
First and foremost, it's important to understand the different between various kinds of time series data. 
Event data come from photon counting instruments. In these data, the time series typically consists of photon arrival times, usually in a particular range of energies that the instrument is sensitive to. Event data are univariate, in that the time series is "how many photons at a given time", or "how many photons in a given chunk of time". 
Point measurements are measurements of a (typically) continuous source at a given moment in time, often with some uncertainty associated with the measurement. These data are multivariate, as your time series relates time, your measurement (e.g. flux, magnitude, etc) and its associated uncertainty to one another. 
Problem 4a
Let's look at some event data from BATSE, a high energy astrophysics experiment that flew on NASA's Compton Gamma-Ray Observatory. BATSE primarily studied gamma ray bursts (GRBs), capturing its detections in four energy channels: ~25-55 keV, 55-110 keV, 110-320 keV, and >320 keV.
You have been given four text files that record one of the BATSE GRB detections. Please read these data in.
End of explanation
# complete
Explanation: Problem 4b
When you reach this point, you and your partner should pick a number between 1 and 4; your number is the channel whose data you will work with. 
Using the data for your channel, please visualize the photon events in both a standard histogram and using Bayesian Blocks.
End of explanation
# execute this cell
kplr_hdul = fits.open('./data/kplr009726699-2011271113734_llc.fits')
Explanation: Problem 4c
Let's take a moment to reflect on the differences between these two representations of our data. 
Please discuss with your partner:
* How many bursts are present in these two representations? 
* How accurately would you be able to identify the time of the burst(s) from these representations? What about other quantities?
For the groups who worked with Channel 3:
You may have noticed very sharp features in your blocks representation of the data. Are they real?
To quote Jake VdP:
Simply put, there are spikes because the piecewise constant likelihood model says that spikes are favored. By saying that the spikes seem unphysical, you are effectively adding a prior on the model based on your intuition of what it should look like.
To quote Jeff Scargle:
Trust the algorithm!
Problem 5: Finding flares
As we have just seen, Bayesian Blocks can be very useful for finding transient events-- it worked great on our photon counts from BATSE! Let's try it on some slightly more complicated data: lightcurves from NASA's Kepler mission. Kepler's data consists primarily of point measures (rather than events)-- a Kepler lightcurve is just the change in the brightness of the star over time (with associated uncertainties).
Problem 5a
People often speak of transients (like the GRB we worked with above) and variables (like RR Lyrae or Cepheid stars) as being two completely different categories of changeable astronomical objects. However, some objects exhibit both variability and transient events. Magnetically active stars are one example of these objects: many of them have starspots that rotate into and out of view, creating periodic (or semi-periodic) variability, but they also have flares, magnetic reconnection events that create sudden, rapid changes in the stellar brightness. 
A challenge in identifying flares is that they often appear against a background that is itself variable. While their are many approaches to fitting both quiescent and flare variability (Gaussian processes, which you saw earlier this week, are often used for exactly this purpose!), they can be very time consuming. 
Let's read in some data, and see whether Bayesian Blocks can help us here.
End of explanation
# execute this cell
kplr_hdul.info()
Explanation: Cool, we have loaded in the FITS file. Let's look at what's in it:
End of explanation
# execute this cell
lcdata = kplr_hdul[1].data
lcdata.columns
# execute this cell 
t = lcdata['TIME']
f = lcdata['PDCSAP_FLUX']
e = lcdata['PDCSAP_FLUX_ERR']
t = t[~np.isnan(f)]
e = e[~np.isnan(f)]
f = f[~np.isnan(f)]
nf = f / np.median(f)
ne = e / np.median(f)
Explanation: We want the light curve, so let's check out what's in that part of the file!
End of explanation
# complete 
Explanation: Problem 5b
Use a scatter plot to visualize the Kepler lightcurve. I strongly suggest you try displaying it at different scales by zooming in (or playing with the axis limits), so that you can get a better sense of the shape of the lightcurve.
End of explanation
# complete
edges = 
#
#
#
plt.step(
Explanation: Problem 5c
These data consist of a variable background, with occasional bright points caused by stellar flares.
Brainstorm possible approaches to find the flare events in this data. Write down your ideas, and discuss their potential advantages, disadvantages, and any pitfalls that you think might arise. 
Discuss with your partner
Solution 5c
Write your notes here
There are lots of possible ways to approach this problem. In the literature, a very common traditional approach has been to fit the background variability while ignoring the outliers, then to subtract the background fit, and flag any point beyond some threshold value as belonging to a flare. More sophisticated approaches also exist, but they are often quite time consuming (and in many cases, detailed fits require a good starting estimate for the locations of flare events). 
Recall, however, that Bayesian Blocks is particularly effective at identifying change points in our data. Let's see if it can help us in this case!
Problem 5d
Use Bayesian Blocks to visualize the Kepler lightcurve. Note that you are now using data that consists of point measurements, rather than event data (as for the BATSE example).
End of explanation
# complete
your_data = 
# complete
def stratified_bayesian_blocks(
Explanation: Concluding Remarks
As you can see, using Bayesian Blocks allowed us to represent the data, including both quiescent variability and flares, without having to smooth, clip, or otherwise alter the data. 
One potential drawback here is that the change points don't identify the flares themselves, or at least don't identify them as being different from the background variability-- the algorithm identifies change points, but does not know anything about what change points might be interesting to you in particular. 
Another potential drawback is that it is possible that the Bayesian Block representation may not catch all events, or at least may not, on its own, provide an unambiguous sign that a subtle signal of interest is in the data.
In these cases, it is sometimes instructive to use a hybrid approach, where one combines the bins determined from a traditional histogram with the Bayesian Blocks change points. Alternatively, if one has a good model for the background, one can compared the blocks representation of a background-only simulated data set with one containing both background and model signal. 
Further interesting examples (in the area of high energy physics and astrophysics) are provided by this paper:
https://arxiv.org/abs/1708.00810
Challenge Problem
Bayesian Blocks is so great! However, there are occasions in which it doesn't perform all that well-- particularly when there are large numbers of repeated values in the data. Jan Florjanczyk, a senior data scientist at Netflix, has written up a description of the problem, and implemented a version of Bayesian Blocks that does a better job on data with repeating values. Read his blog post on "Stratfied Bayesian Blocks" and try implementing it! 
https://medium.com/@janplus/stratified-bayesian-blocks-2bd77c1e6cc7
You can either apply this to your favorite data set, or you can use it on the stellar parameters table associated with the data set we pulled from MAST earlier (recall that you previously worked with the exoplanet parameters table). You can check out the fields you can search on using MAST CasJobs, or just download a FITS file of the full table, here: https://archive.stsci.edu/prepds/kg-radii/
<sub>This GIF is of Dan Shiffman, who has a Youtube channel called Coding Train</sub>
End of explanation
<END_TASK>