\section{Introduction}
\label{sec:intro}
Astronomical sources that exhibit temporal variability or periodicity are among the most scientifically rich in the observed Universe.  Variation in source intensity can result from rotations, eclipses, magnetic activity cycles, or the turbulent flows of matter into deep gravitational wells.  Massive stellar explosions known as {\it super novae}, for example,  appear as abrupt peaks in the electromagnetic radiation emitted from a source. They sometime produce a narrow beam of intense high-energy radiation that appears as $\gamma$-ray bursts from Earth---although there are other sources of $\gamma$-ray bursts. These bursts originate in distant galaxies, can last from milliseconds to minutes, and are among the most energetic events that occur in the Universe. There are many less dramatic objects that produce temporal variability; they include signals with repeated peaks and troughs that may or may not follow a predictable periodic pattern. A {\it variable star}, for example, exhibits fluctuating emission that can be due to a natural expansion and contraction in its radius as it evolves (a so-called {\it pulsating variable star}) or to eruptions in its coronae such as flares or mass ejections. In extreme cases, superflares can erupt producing millions of times more energy than a typical solar flare; if such a flare occured on the Sun it could destroy the Earth's ozone layer and cause mass planetary extinction.  Although the mechanism is not well understood, super flares are most likely to occur on fast rotating stars \citep{mcke:12}.

In this article we develop statistical methods that enable us to identify changes in observed emission from astronomical sources, such as those associated with massive energetic events in the coronae of variable stars. Dramatic shifts in intensity of these sources are typically accompanied by changes in the distribution of the wavelength of the photons emitted, known as the {\it spectrum} of the source. Thus, our methods aim to take advantage of changes both in the intensity of the emission and in its wavelength distribution.  We illustrate our methods by applying them to high-energy observations of the flaring star known as  FK\,Com.

 \paragraph{\it Data collection in high-energy astrophysics}
High-energy astrophysics is the study of electromagnetic radiation in the ultraviolet, X-ray, and $\gamma$-ray bands. The production of such high-energy photons requires temperatures of millions of degrees and signals the release of deep wells of stored energy such as those in very strong magnetic fields, extreme gravity, explosive nuclear forces, and shock waves in hot plasmas. High-energy detectors have much in common in terms of their statistical properties (and much that distinguishes them from detectors in other energy regimes). In this article we focus on statistical methods that are appropriate for high-energy astrophysics.

Astronomical data from high-energy observatories are usually obtained as a list of photons;  the list records the two-dimensional direction from whence each photon arrived, the time at which it was recorded, and its wavelength (and hence its energy).  \changes{Owing to the intrinsic resolutions of high-energy photon telescope/detector systems, each of these four attributes is inherently discrete and the data can be compiled into a four-way table of photon counts.}  Although each attribute is subject to uncertainty (also inherent in the telescope/detector systems, e.g., the true direction of arrival is distorted by the point spread function), the quality of the resulting datasets is unprecedented in the history of astronomical observations.  Nonetheless the remarkable four-way tables of photon attributes are largely an untapped resource.  This is mostly due to a lack of principled statistical methods that can be used to detect and characterize astronomical sources in high-dimensional spaces.  This is especially true of high-spectral-resolution grating data (e.g., data obtained using the HETGS+ACIS-S configuration of the {\it Chandra X-ray Observatory}), where a self-selected sample of {\sl interesting} objects are observed for long durations (often for $>100$~ksec) with Doppler resolutions of ${\sim}100$~km~s$^{-1}$, and with an intrinsic temporal resolution that varies from one millisecond to three seconds.
Because the locations of the objects we study in this article do not change,  we ignore directional information and focus on the two-way table of photon counts indexed by time and wavelength. The one-way marginal table of wavelengths is called the (observed) {\it spectrum} and that of the times is the (observed) {\it light curve} of the source.

\paragraph{\it Change points in Poisson processes}
From a statistical point of view, sudden changes in an underlying physical process are modeled via change points. In a parametric model, for example, the value of the parameter is allowed to change at one or more time points during the observation. The periods of time where the parameter is fixed are called the {\it regimes} of the process. We might consider testing for a single change point, fitting the model with a given number of change points, or estimating the number of change points.  Given that the data in high-energy astrophysics are photon counts, we focus on change points in non-homogeneous Poisson processes. Because the wavelength (and direction) \changes{of} each photon are recorded, we can either model them using a marked Poisson process or count the events in each of a number of wavelength bins and model them as a multivariate Poisson process. Typical detectors have discrete temporal and spectral resolution so it is natural to use photon counts in time cross wavelength bins, i.e., to use discrete-time processes.
\changes{The methods that we propose are designed specifically for two-way tables of Poisson counts of this sort.}

The complexity of any statistical model naturally increases with the number of change points. Thus to avoid overfitting, methods must penalize complexity. Bayesian prior distributions and marginalization techniques are natural avenues for encouraging parsimony, for example in that they inherently embody Occam's Razor in model selection  \citep[e.g.,][chapter 28]{mack:03}. It is not surprising then that there is an extensive literature on Bayesian methods for change point models. \citet{raft:akma:86} is an early contribution for count data generated according to a homogeneous Poisson process until some unknown time when \changes{the} intensity of the process changes.  They propose computing the joint posterior distribution of the two Poisson rates and the single change point along with a Bayes Factor to test for existence of the change point; see \citet{carl:etal:92}  and \citet{more:etal:05} for derivation of hierarchal and intrinsic prior distributions, respectively, in the single change point Poisson setting.  \citet{gree:95} generalizes this approach by using a continuous time model, allowing for multiple regimes, putting a prior distribution on the number of regimes, and fitting this number.  \citet{lai:xing:11}, on the other hand, use a discrete time model with a Bernoulli process governing the change in regimes. This approach follows \citet{chib:98} who embedded a latent discrete-time discrete-state Markov process for the regime index into a multi-level Bayesian model. This formulation has been applied for inhomogeneous Poisson processes, for example, to model change points in the parameters of a log-linear model \citep{park:10}, see \citet{park:etal:12} for a related approach.

Perhaps the most popular change-point method for astronomical count data is known as Bayesian Blocks \citep{scar:98}. It starts by splitting the time interval into two, assuming a homogeneous Poisson process on each, and computing the posterior odds of this model relative to a homogeneous process on the entire interval.  The change point is chosen by maximizing the posterior odds and is accepted if the odds favor this model. This process continues recursively on each subinterval until no more change points are added. \citet{scar:etal:13} improved the method by achieving global optimization over the change points while adding the capacity to handle gaps in the data,  variable exposure, piecewise linear and piecewise exponential Poisson intensities, multivariate data, the analysis of variance,  data on the circle, etc. Bayesian Blocks does not, however, address changes in the spectrum at the change points.

Non-Bayesian methods for change-point detection have also been proposed for non-homogeneous Poisson processes.  \citet{Akman-Raftery86} and \citet{Worsley86} are early contributions.  Both papers developed methods for testing for the presence of a change-point in a Poisson process.  These methods follow classical frequentist arguments and also provide interval estimates for any detected change-points.  \citet{Loader92} considered using log-linear models for non-homogeneous Poisson processes.  In particular likelihood ratio tests are derived to choose a ``best'' model amongst different change-point and log-linear models for the observed data.  In \citet{Mei-et-al11} methods are proposed for detecting changes in Poisson rate after the effects of population size are taken into account.  These methods are based on a generalized likelihood ratio, weighted likelihood ratio, and adaptive thresholding.  More recently, \citet{shen:zhan:12} applied the modified BIC (mBIC) of \citet{Zhang-Siegmund07} to detect change-points in non-homogeneous Poisson processes.
\changes{
The mBIC can be viewed as a penalized likelihood model selection criterion that is tailored to estimating the number and locations of the change points.
}

\paragraph{\it Change points in marked Poisson processes}
Despite the extensive literature on change points in Poisson processes,
there is little on change points in {\it marked} Poisson processes.
\changes{Methods  have been developed to test for non-homogeneity in marked Poisson processes, for example via a scan statistic that is computed in a sliding time interval \citep{chan:zhan:07}.
}
 In high-energy astrophysics a common strategy is to compute a simple summary statistic describing the wavelength distribution in each time bin and to visually inspect a plot of the statistic as a function of time. A typical choice of statistics is the so-called {\it hardness ratio} which is the ratio of the photon count observed above a given wavelength threshold to that below the threshold; see \citet{park:etal:06} for a small-sample statistical treatment of hardness ratios. Plotting hardness ratios against time is an informal exploratory technique. The primary goal of this article is to provide a change point method for marked Poisson processes that is inspired by and tailored to a specific problem in high-energy astrophysics.

\paragraph{\it \changes{Methods for detecting change points}}
\changes{
There is a wide range of  statistical approaches for detecting statistical change points.  For the classical simple case of piecewise constant signals with additive noise, for example, \citet{Yao88} studied the use of the Schwartz criterion and showed it to be consistent.  This criterion can be interpreted as the classical Bayesian information criterion (BIC) if the location of a change point is treated as a free parameter (or as a dimension of the model). In this case, however, the likelihood functions violate the regularity conditions necessary for the classical BIC, so both criteria can likely be improved. To circumvent this issue, \citet{Zhang-Siegmund07} provide an ingenious method to approximate the Bayes factor and developed the above-mentioned mBIC for change point detection which has been shown to be consistent in various settings.  More recently, \citet{Ninomiya14} derived a specialized Akaike information criterion (AIC) for change point estimation.  Results from \citet[][Section~3.3]{Aue-Lee11}, however, can be used to show that this AIC method is not consistent.

In this article we focus on the minimum description length (MDL) principle \citep{Rissanen89,Rissanen07}, which can also be viewed as a penalized likelihood approach to model selection. A brief introduction of MDL is given in Section \ref{sec:firstMDL}. One of its earliest applications to change point estimation is the image segmentation work of \citet{Leclerc89}, where the corresponding MDL criterion can be easily reduced to the classical one-dimensional change point problem setting; see also \citet{Lee98:segcor} for an extension to correlated noise.  Other successful MDL applications to change point problems include discontinuity detection in nonparametric regression \citep{Lee02}, structural break detection in nonstationary time series models \citep{Davis-Yau13,Yau-et-al15}, and change point detection in quantile modeling \citep{Aue-et-al14b}.  Under mild regularity assumptions, these MDL solutions can be shown to be consistent.

For many change-point problems simpler than the one considered in this paper, the major penalty terms of mBIC and MDL are remarkably similar; e.g., compare Equation~(8) of \citet{Zhang-Siegmund07} and Equation~(6.2) of \citet{Lee97c}.  This provides reassuring evidence that both mBIC and MDL are reliable methods for change point detection.}

\paragraph{\it \changes{Strategy and outline}} We pursue a non-parametric penalized-likelihood strategy, using MDL to optimize the tuning parameters of the penalty function and to fit the number of change points.  Between the change points, we assume a homogeneous Poisson process where the spectrum includes a broad smooth term that is modeled using a radial basis expansion.  In particular the bases we use are cubic polynomials, which provide a good balance between model flexibility and ease of implementation.  We also include a number of single-bin local deviations from this extended smooth model. Physically, these deviations correspond to spectral emission lines that are the result of the production of photons of a particular wavelength as electrons fall to the lower energy shells of a specific ion. Both the radical basis expansion and the local deviations are over parameterized in that we expect most of their coefficients and intensities (respectively) to be zero. Thus, we employ an $\ell_1$ penalty on both.
An R package ``Automark" implementing our proposed methodology
is available to interested readers at
{\tt https://github.com/astrostat/Automark}.

Some simple spectral models such as power laws and exponential absorption can be formulated as log-linear models \citep{vand:etal:01} and thus in principle could be embedded into a temporal model with change points for their parameters, similar to the proposal of \citet{park:10}.  We take a different approach by using flexible non-parametric spectral models. This allows us to account for both blurring of the recorded energies and background contamination of the photon counts, and to apply a single model irrespective of the particular shape of the source spectra. (See \citet{vand:kang:04} for a review of the source and data-distortion models used in high-energy spectral analysis.)  This strategy is contrary to the typical approach in high-energy spectral analysis where physics-based parametric models are
preferred because their parameters have direct scientific meaning. Our primary goal here, however, is to identify the change points. The spectra corresponding to each regime within the process can be fit using physics-based models in a secondary analysis.

This article is organized into seven sections. Section~\ref{sec:back} describes the detectors used for data collection, the relevant astrophysical models, and basic notation for our statistical approach. We formalize our model  for the homogeneous Poisson processes within each regime of the overall process in Section~\ref{sec:hom} and describe how we embed this into the change-point model in Section~\ref{sec:break}. Numerical results including a set of simulation studies and an application to the yellow giant FK\,Com\ appear in Sections~\ref{sec:sim} and \ref{sec:real}, respectively. Discussion appears in Section~\ref{sec:conc} and details of the derivation of our MDL criterion in an Appendix.

\section{Astrophysical Background and Statistical Notation}
\label{sec:back}

Suppose that we observe photon counts in an $N\times J$ rectangular array of equally-sized wavelength-cross-time bins, where $y_{ij}$ is the photon count in the bin with wavelength range $[w_i-\delta_w/2, w_i+ \delta_w/2)$ and
temporal range $[t_j-\delta_t/2, t_j+\delta_t/2)$, for $i=1,\ldots, N$ and $j=1,\ldots, J$. \changes{Although this array may be compiled using the natural resolution of the detector, with high-resolution data it may be computationally advantageous to combine the natural bins in order to obtain a lower resolution array of photon counts. This strategy is illustrated and discussed in Section~\ref{sec:real}.}
In any case, the expected count in each bin is ideally proportional to the brightness of the astronomical source integrated over its time and wavelength ranges. Because of detector effects, however, a photon with wavelength $w$ has a distribution of possible recorded energies, $w'$. This distribution is called the {\it detector response function} and denoted $R(w,w')$. We apply our methods to data collected using a  {\it grating} that, like a prism,  spatially separates photons as a function of their wavelength. With grating data, wavelength is measured very accurately and we can ignore measurement uncertainty and treat $R(w,w')$ as a Dirac delta function. Another detector effect arises because the sensitivity of the detector varies with photon wavelength, an effect that is quantified via the so-called {\it effective area} of the detector and denoted by $A(w)$. The effective area quantifies characteristics of the telescope mirrors that do not relate to the nature of the source and is typically treated as known. (See, however, \citet{lee:etal:11} and \citet{xu:etal:14} for development of Bayesian statistical methods that account for uncertainty in $A(w)$.)

Using a Poisson model for the photon counts and letting $\lambda(\tm, \wl)$ be the expected count per unit time and per unit wavelength averaged over the bin centered at $(t_j,w_i)$,

\begin{equation}
y_{ij} \buildrel{\rm ind}\over\sim \mathrm{Poisson}\{ \delta_t \cdot  \delta_w  \cdot \lambda(\tm, \wl) \cdot A(w_i)\}.
\label{eq:pois-mdl}
\end{equation}
In practice, we may combine data from $K$ detectors, each with its own effective area. Letting  $A_1(w), \ldots, A_K(w)$ denote the $K$ effective areas, the total counts across the $K$ detectors can be modeled
\[
Y_{ij} = \sum_{k=1}^K y_{ijk},
\]
where
\[
y_{ijk} \buildrel{\rm ind}\over\sim \mathrm{Poisson}\{ \delta_t \cdot \delta_w \cdot \lambda(\tm, \wl) \cdot A_k(w_i)\}
\]
are the counts \changes{recorded with} detector $k$.  Because the $y_{ijk}$ are independent  across $k$, we have
\[
Y_{ij} \sim \mathrm{Poisson}\left\{  \delta_t \cdot \delta_w \cdot \lambda(\tm, \wl) \sum_{k=1}^K A_k(w_i)\right\}.
\]

The goal of this article is to infer the properties of $\lambda(\tm, \wl)$ using the observed photon counts, $\{Y_{ij}\}$.  We are particularly interested in detecting statistically significant changes in $\{Y_{ij}\}$ over time.  In other words, we would like to determine if there are any change point, $\pi$, such that $\lambda(\tm, \wl)\big|_{\{t_j\leq \pi\}}$ \changes{differs} from $\lambda(\tm, \wl)\big|_{\{t_j>\pi\}}$.  If there are, we aim to estimate the number of change points and their values.

For ease of presentation, we develop our estimation procedure in two stages.  First  we develop a flexible model for a time-homogeneous source spectrum. That is we assume that there are no change points and that $\lambda(\tm, \wl)$ is constant in $t_j$.
In the second stage we  introduce change points and allow $\lambda(\tm, \wl)$ to change over time.

\section{A Model for the Homogeneous Poisson Processes}
\label{sec:hom}

\subsection{A nonparametric spectral model}
Because there are no change points in our time-homogeneous Poisson model, we
drop $t_j$ from the argument of $\lambda(\tm, \wl)$ and write $\lambda(\tm, \wl)=\lambda(\wl)$
in this section. Also, we denote the expectation of $Y_{ij}$ by $\mu(w_i)$.
That is, $\mu(w_i) = s(w_i) \lambda(w_i)$ with
$s(w_i) = \delta_t\cdot\delta_w\sum^K_{k=1}A_k(w_i)$ being a completely specified function.
We represent the expected counts due to the underlying physical process by $\lambda$ and  expected detector counts by $\mu$; $\mu$ is adjusted for varying bin sizes and effective areas.

The parameter of scientific interest is $\lambda(\wl)$. For the astronomical sources \changes{that} we study $\lambda(\wl)$ is mostly smooth with a few emission lines overlaid and we adopt the model
\begin{equation}
\lambda(\wl) = f(w_i)+\sum^M_{m=1}\alpha_m g_m(w_i,\tau_m), \label{fun:lambdaw}
\end{equation}
where $f(w)$ is a nonnegative smooth function, $M$ is the number of emission lines, $\alpha_m > 0$ is the magnitude of emission line $m$, and $g_m(w,\tau_m)$ models the shape of the emission line centered at $\tau_m$.
Common choices for $g_m(w,\tau_m)$ include Gaussian density functions with small
variances, Lorentzian density functions, and delta functions
\citep[e.g.,][and the references therein]{Dyk-Connors-Esch06}.
Here we assume that each of the emission lines locates completely within a single bin, although our framework can be easily modified for cases when some emission lines span multiple bins.
For astrophysical reasons, we also assume that each $\alpha_m >0$, that is each emission line constitutes a positive deviation from $f(w)$, although in principle  the $\alpha_m$ may also be negative, so long as $\lambda(\wl)>0$.

If an individual spectral line is located in a single bin, the line profile, $g_m(w, \tau_m)$, is an indicator function for the interval of width $\delta_w$ centered at $\tau_m\in\{w_1,\dots,w_N\}$, that is, an indicator function for one of the wavelength bins. The natural resolution of the bin counts does not allow modeling of  $g_m(w, \tau_m)$ at any greater level of detail than this.
Under this model for the emission lines,  $\lambda(\wl)$ in (\ref{fun:lambdaw}) becomes
\[
\lambda(w_i) =
 f(w_i) + \sum^M_{m=1}\alpha_m I\left(\tau_m-\frac{\delta_w}{2} \le w_i < \tau_m+\frac{\delta_w}{2}\right),
\]
where
 $I(\mathcal{A})$ denotes the indicator function for $\mathcal{A}$ and  $\tau_m\in\{w_1,\dots,w_N\}$.
Because of the nonnegativity constraint on $\lambda(w_i)$, we introduce a log-link function,
\begin{equation}
  \log\lambda(w_i) =  \log f(w_i) + \sum^M_{m=1}\alpha_m' I\left(\tau_m-\frac{\delta_w}{2} \le w_i < \tau_m+\frac{\delta_w}{2}\right)
 \label{eq:log-link}
\end{equation}
 where $\alpha_m'$ is the transformed  $\alpha_m$ that ensures the equality in (\ref{eq:log-link}) holds.

To provide flexible modeling, we use a radial basis expansion \changes{\citep[see, e.g.,][]{Buhmann03, Ruppert-Wand-Carroll03}} to nonparametrically model $\log f(w)$, the smooth component of $\log\lambda(w)$.  Specifically, we use the $P$ basis functions, $(\xi_1(w), \ldots, \xi_P(w))$, which are polynomials of power 3 with \changes{equally spaced} knots, $\kappa_1, \ldots, \kappa_{P-4}$.\footnote{\changes{The gap between the first knot the left endpoint of the spectral range and the gap between the last knot and the right endpoint of the range are both slightly larger than the gaps between consecutive knots.}
}
\changes{Specifically, the basis functions are}
\begin{align*}
&\xi_1(w) = 1,\quad \xi_2(w) = w, \quad \xi_3(w)=w^2,\quad \xi_4(w)=w^3,\quad \mbox{and}\\
&\xi_{p+4}(w) = |w -\kappa_p|^3, \quad\quad p=1,\dots,P-4.
\end{align*}
\changes{
Substituting these into (\ref{eq:log-link}) yields}
\begin{align}
\log\lambda(\wl)  & =
  \sum_{p=1}^{P}\beta_p \xi_p(w_i) + \sum^M_{m=1}\alpha_m' I\left(\tau_m-\frac{\delta_w}{2} \le w_i < \tau_m+\frac{\delta_w}{2}\right),
\label{eqn:gensemi}
\end{align}
where
$\beta_p$ is the coefficient of basis function $p$.

Model~(\ref{eqn:gensemi}) is an example of a generalized semi-parametric model. Under this formulation, estimating $\lambda(\wl)$ is equivalent to estimating $P$, $(\beta_1, \ldots, \beta_P)$, $M$, and $(\alpha_1', \tau_1, \ldots, \alpha_M', \tau_M)$.

\subsection{A lasso model for the emission lines}
A difficulty in fitting model~(\ref{eqn:gensemi}) is the estimation of $(\tau_1, \ldots, \tau_M)$, the locations of the emission lines.  With $M$  unknown, there is a total of $2^N$ possibilities for $(\tau_1,\ldots,  \tau_M)$.  However, as emission lines are relatively rare (i.e., sparse), we can take advantage of the celebrated $\ell_1$ penalty techniques \citep[e.g.,][]{Tibshirani96} to provide a fast searching algorithm for $(\tau_1, \ldots, \tau_M)$.

We begin by rewriting~(\ref{eqn:gensemi}) as
\begin{equation}
  \log\lambda(\wl)
  = \sum^P_{p=1}\beta_p \xi_p(w_i) + \sum^N_{k=1} \eta_k I_{k}(w_i)
  = \sum^P_{p=1}\beta_p \xi_p(w_i) +\eta_i,
\label{eqn:pmodel}
\end{equation}
with $I_{k}(w) = I(w_k-\delta_w/2 \le w< w_k+\delta_w/2)$.  In this formulation $\eta_i=0$ implies there is no emission line in the wavelength bin centered at $w_i$.  There is a one-to-one correspondence between models~(\ref{eqn:gensemi}) and~(\ref{eqn:pmodel}) via
\begin{align*}
M = \sum^N_{i=1} I(\eta_i\neq 0) \quad \mbox{and} \quad
\Big\{(\tau_m,\alpha_m')\Big\}^M_{m=1} = \Big\{(w_i, \eta_i): \eta_i\neq0\Big\}.
\end{align*}
A major advantage of re-expressing model~(\ref{eqn:gensemi}) as (\ref{eqn:pmodel}) is that, many of the $\bm{\eta}=(\eta_1,\ldots, \eta_N)^\intercal$ are zero in~(\ref{eqn:pmodel}). Therefore by imposing an $\ell_1$ penalty, we can achieve fast simultaneous selection and estimation of the nonzero elements among $\bm{\eta}$.

We employ a common strategy for the smooth component, $\log f(w_i)=\sum_{p=1}^P \beta_p \xi_p(w_i)$, in~(\ref{eqn:pmodel}).  Namely, we specify a large value of $P$ {\it a priori} and estimate $\bm{\beta}=(\beta_1, \ldots, \beta_P)^\intercal$ under a penalty term that aims to avoid over-fitting.  In this way, the coefficients of the the redundant functions among $\xi_1(w), \ldots, \xi_P(w)$ are shrunk toward zero.  This enables us to assume that $P$ is known (and large).

Penalized maximum likelihood can be used to simultaneously estimate all of the unknown parameters in~(\ref{eqn:pmodel}), including those describing the smooth component and those describing the emission lines.

For a Poisson random variable with expectation $a>0$, we write the corresponding
probability mass function as
\begin{equation}
q(x;a)=-a + x\log a -\log(x!)
\label{eq:pois-loglike}
\end{equation}
for any
nonnegative integer $x$.
The log-likelihood of a bin count $Y_{ij}$ is
$$
l_{\mathrm{one}}(Y_{ij}; \bm{\theta}) = q\{Y_{ij}; \mu(w_i)\} = q\{Y_{ij}; s(w_i) \lambda(w_i)\},
$$
where $\lambda(\wl)$ is parameterized in terms of $\bm{\theta}=(\bm{\beta},\bm{\eta})$ as in (\ref{eqn:pmodel}).
(This likelihood will be used within one homogenous time segement in Section~\ref{sec:break}; hence the subscript ``one".)
Adding a penalty, we define the estimate of $\bm{\theta}$ as the maximizer of
\begin{equation}
\sum^N_{i=1}\sum^J_{j=1}l_{\mathrm{one}}(Y_{ij}; \bm{\theta}) - \gamma \{\rho\|\brm{D}\bm{\beta}\|_1 +
(1-\rho) \|\bm{\eta}\|_1\}, \label{eqn:penlike}
\end{equation}
where $\brm{D}=\mathrm{diag}\{(0,0,0,0,1,\dots,1)\}$ is a diagonal matrix, and $\gamma>0$ and $0\le\rho\le1$ are tuning parameters that determine the penalty on $\bm{\theta}$. {The first four diagonal elements of $\brm D$ are set to zero because there is no penalty for the inclusion of $(\beta_1, \ldots, \beta_4)$ in the model.}   If $\gamma$ and $\rho$ are specified, maximization of~(\ref{eqn:penlike}) is equivalent to lasso fitting of a generalized linear model. There are fast algorithms for this optimization problem \citep[e.g., coordinate descent,][]{Friedman-Hastie-Tibshirani10} and software is widely available; e.g., the $R$ package {\ttfamily glmnet}.

\subsection{Tuning parameter selection using MDL}
\label{sec:firstMDL}
The success of~(\ref{eqn:penlike}) depends heavily on our ability to choose good values of $\gamma$ and $\rho$, which determine the number of emission lines and the number of basis functions used in the smooth component of  the fitted model. To select $\gamma$ and $\rho$, we use the minimum description length (MDL) principle \citep{Rissanen89,Rissanen07}.  In short, MDL defines the best model as the one that achieves the highest compression rate of the data, $\mathcal{D}$.  In other words, the best model allows us to store  $\mathcal{D}$  with the shortest codelength.  A good statistical model and a good compression method share a common feature: they should be able to capture regularities and patterns hidden in the data.  Therefore it is reasonable to expect that a model chosen by MDL  to be a good statistical model;
\changes{some successful examples for change point problems were provided in Section~\ref{sec:intro}, while for other applications can be found in the later chapters of \citet{grunwald2005advances}.  There are various versions of MDL.}
We  use the so-called two-part form of MDL to derive a model selection criterion for choosing $\gamma$ and $\rho$.

Let $\mathsf{CL}(\mathcal{Z})$ be the codelength of $\mathcal{Z}$: one can treat $\mathsf{CL}(\mathcal{Z})$ as the amount of memory needed to store $\mathcal{Z}$.
Similarly, let $\mathsf{CL}(\mathcal{Z}|\mathcal{X})$ be the codelength of $\mathcal{Z}$
  conditional on $\mathcal{X}$, { that is, the codelength of $\mathcal{Z}$ with knowledge of $\mathcal{X}$.}
Generally speaking, the two-part form of MDL decomposes the codelength of the data $\mathcal{D}$ into two parts
\begin{equation}
  \mathsf{CL}(\mathcal{D}) = \mathsf{CL}(\widehat{\mathcal{M}}) +
  \mathsf{CL}(\mathcal{D}|\widehat{\mathcal{M}}), \label{eqn:cl}
\end{equation}
where the best fitting model $\widehat{\mathcal{M}}$ is defined to be  the minimizer of $\mathsf{CL}(\mathcal{D})$.
In Appendix~\ref{app:mdl1} it is shown that when fitting model~(\ref{eqn:penlike}), \changes{for large $N$} the two-part MDL (\ref{eqn:cl})
\changes{is asymptotically equal to}
\begin{equation}
  \mathsf{mdl}_{\mathrm{null}}(\gamma, \rho) =
  -\sum^N_{i=1}\sum^J_{j=1}l_{\mathrm{one}}(Y_{ij};\hat{\bm{\theta}}) +
  \frac{1}{2}\|\hat{\bm{\theta}}\|_0\log (NJ) +
\log {N \choose \|\hat{\bm{\eta}}\|_0}, \label{eqn:mdl_hom}
\end{equation}
where $\hat{\bm{\theta}}=(\hat{\bm{\beta}}, \hat{\bm{\eta}})$ is the maximizer of~(\ref{eqn:penlike}) given $\gamma$ and $\rho$, and the $\ell_0$ norm, $\|\bm{z}\|_0$, denotes the number of nonzero elements in the vector $\bm{z}$.  \changes{We choose the values of $\gamma$ and $\rho$ that jointly maximize (\ref{eqn:mdl_hom}) by evaluating $\mathsf{mdl}_{\mathrm{null}}(\gamma, \rho)$ on a fine grid, using {\tt glmnet} to compute $\hat{\bm{\theta}}(\gamma,\rho)$ on a grid of values of $\gamma$ for each value of $\rho$.}
In many MDL applications, a term similar to the last one in~(\ref{eqn:mdl_hom}) is  negligible and hence omitted.  Here, however, the number of unknown parameters ($N+P$) is not ignorable when compared to the number of observations and hence the last term  in~(\ref{eqn:mdl_hom}) is required.

In summary, for the time-homogeneous model with no change points, the parameter estimate $\hat{\bm{\theta}}=(\hat{\bm{\beta}}, \hat{\bm{\eta}})$ is obtained by maximizing~(\ref{eqn:penlike}) with $\gamma$ and $\rho$ chosen as the joint minimizer of~(\ref{eqn:mdl_hom}).  Finally, the estimate for $\lambda(w_i)$ is calculated as $\hat{\lambda}(w_i)=\exp\{\sum^P_{p=1}\hat{\beta}_p\xi_p(w_i)+\hat{\eta}_i\}$.

\section{Modeling with Change Points}
\label{sec:break}

\subsection{Adding change points to the spectral model}
\label{sec:addchange}
In this section, we allow the spectrum, $\lambda(\tm, \wl)$,  to change over time. Hence we reinstate the notation that emphasized the dependence on $t_j$. We model $\lambda(\tm, \wl)$ as piecewise constant as a function of $t_j$.  This involves partitioning the entire time interval $[0,T)$ into several subintervals and independently modeling $\lambda(\tm, \wl)$ as constant in $t_j$ in each subinterval using a model of the form given in~(\ref{eqn:pmodel}).
Again let $Y_{ij}$ be the photon count in the bin centered at $(w_i, t_j)$ and summed over the $K$ detectors. We can formalize the change-point model for $\lambda(\tm, \wl)$ as
\begin{equation}
\lambda(\tm, \wl) = \sum^B_{b=1} I_b(t_j) \lambda_b(\wl),
\label{eqn:piece}
\end{equation}
where
$B\ge 1$ is the number of time (sub)intervals,
$I_b(t_j)=I( \pi_{b-1} \le t_j < \pi_{b} )$ identifies the time bins in interval $b$,
$t_1-\delta_t/2=\pi_0 <\dots<\pi_B=t_J + \delta_t/2$ are the change points,
and $\lambda_b(\wl)$ is the spectrum in time interval $b$. {For ease of discussion, we refer to $\pi_1,\dots, \pi_{B-1}$ as the {\it change points}, excluding the endpoints $\pi_0$ and $\pi_B$\changes{.}  The within time-interval spectra are modeled as in~(\ref{eqn:pmodel}), i.e.,
\[
  \log\lambda_b(\wl)  = \sum^P_{p=1}\beta_{bp} \xi_p(w_i) +\eta_{bi},
\]
where $\bm{\beta}_b=(\beta_{b1}, \ldots, \beta_{bP})^\intercal$ and $\bm{\eta}_b=(\eta_{b1}, \ldots, \eta_{bN})^\intercal$ are the  parameters for the time-interval-specific spectra.

There is no gain in allowing the change points, $(\pi_1,\ldots, \pi_{B-1})$, to take values other than
bin breaks if there is no prior knowledge about the change points.
Further, to maintain sufficient data within each time interval for acceptable estimation, $|\pi_i-\pi_j|$ cannot be too small for any $0\le i,j\le B$.  We require that $|\pi_i-\pi_j|\ge 5\delta_t$ in our numerical illustrations. This ensures that each time interval is at least five bins wide.

\subsection{Model selection using MDL}
\label{sec:finalMDL}

In order to fit model~(\ref{eqn:piece}), we must estimate the number of time intervals, $B$, the change points $\bm{\pi}=(\pi_1,\dots,\pi_{B-1})^\intercal$, and the parameters \changes{for each interval,} $\bm{\Theta}=\{\bm{\theta}_1,\dots,\bm{\theta}_B\}$, where $\bm{\theta}_b = (\bm{\beta}_b, \bm{\eta}_b)$.  Once $B$ and $\bm{\pi}$ are specified, however, each $\lambda_b(\wl)$ along with its tuning parameters, $\gamma=\gamma_b$ and $\rho=\rho_b$, can be estimated individually using the method of Section~\ref{sec:hom}.  Thus, we aim to first estimate $B$, $\bm{\pi}$ and then compute $\hat\gamma_b(B, \bm{\pi})$, $\hat\rho_b(B, \bm{\pi})$, and $\hat\bm{\theta}_b(\hat\gamma_b, \hat\rho_b)$ as in Section~\ref{sec:hom} for $b=1,\ldots, B$.
This  notation signifies the dependence of
$\hat\gamma_b$ and $\hat\rho_b$ on $(B, \bm{\pi})$, and the dependence of $\hat\bm{\theta}_b$ on $(\hat\gamma_b, \hat\rho_b)$.
Similarly,
we also write ${\bm{\Theta}}$ and $\lambda(t_j,w_i)$ as
${\bm{\Theta}}(\bm{\gamma}, \bm{\rho})$ and
$\lambda\left\{t_j,w_i; B, \bm{\pi}, {\bm{\Theta}}(\bm{\gamma}, \bm{\rho})\right\}$,
respectively.

The profile penalized (pp) loglikelihood function for $B$, $\bm{\pi}$  under model~(\ref{eqn:piece}) depends on the tuning parameters $\bm{\gamma}=(\gamma_1,\dots,\gamma_B)$ and $\bm{\rho}=(\rho_1,\dots, \rho_B)$ and can be written
\begin{equation}
\log{\cal L}_{\mathsf{pp}}(B, \bm{\pi}; \bm{\gamma}, \bm{\rho}) =
\sum^N_{i=1}\sum^J_{j=1} q\left(Y_{ij}; s(w_i)\lambda\left\{t_j,w_i; B, \bm{\pi}, \widehat{\bm{\Theta}}(\bm{\gamma}, \bm{\rho})\right\}\right),
\label{eq:pplike}
\end{equation}
\changes{with $q$ the Poisson loglikelihood given in (\ref{eq:pois-loglike}).}
Because $\bm{\gamma}$ and $\bm{\rho}$ can be computed as functions of $B$ and $\bm\pi$, we can substitute
$\hat{\bm\gamma}(B,\bm\pi)$ and $\hat{\bm\rho}(B,\bm\pi)$ into
$\log{\cal L}_{\mathsf{pp}}(B, \bm{\pi}; \bm{\gamma}, \bm{\rho})$.  This is akin to
profiling over $\bm{\gamma}$ and $\bm{\rho}$, except that
$\hat{\bm\gamma}(B,\bm\pi)$ and $\hat{\bm\rho}(B,\bm\pi)$ are computed
according to the MDL criterion rather than by maximizing $\log{\cal L}_{\mathsf{pp}}(B,
\bm{\pi}; \bm{\gamma}, \bm{\rho})$. The result is the pseudo profile penalized
(ppp) loglikelihood given by
\begin{align}
  \log {\cal L}_{\mathsf{ppp}}(B, \bm{\pi}) &= \log{\cal L}_{\mathsf{pp}}(B, \bm{\pi}; \hat{\bm\gamma}(B,\bm\pi), \hat{\bm\rho}(B,\bm\pi))
  \nonumber \\
  &= \sum^N_{i=1}\sum^J_{j=1} \sum^B_{b=1} I_b(t_j) l_{\mathrm{one}}\left\{ Y_{ij};\hat\bm{\theta}_b(B,\bm\pi)  \right\},
\label{eq:llike}
\end{align}
where $\hat\bm{\theta}_b(B,\bm\pi) = \hat\bm{\theta}_b\big(\hat{\bm\gamma}(B,\bm\pi), \hat{\bm\rho}(B,\bm\pi)\big)$.
Because different values of $B$ lead to different numbers of parameters in the model,
we cannot estimate the parameters by maximizing (\ref{eq:llike}).  Instead we view this as a model selection problem and again use the MDL principle. In this way we can choose the ``best'' combination of $B$ and $\bm{\pi}$, and then compute $\hat\gamma_b(B, \bm\pi)$, $\hat\rho_b(B, \bm\pi)$, and $\hat\bm{\theta}_b(\hat\gamma_b, \hat\rho_b)$ \changes{for each time interval}. We derive the MDL criterion  in Appendix~\ref{app:mdl2} and show that the best model can be found by minimizing
\begin{align}
  \mathsf{mdl}_{\rm full}(B, \bm{\pi})
&= -\log {\cal L}_{\mathsf{ppp}}(B, \bm{\pi}) + \log {B} + \sum^{{B}}_{b=1}\log{c_b(\bm\pi)} \nonumber\\
&\quad\quad+
\sum^{{B}}_{b=1}\left[\frac{1}{2}\| \hat\bm{\theta}_b(B,\bm\pi)\|_0\log\{Nc_b(\bm\pi)\} + \log {N \choose \|\hat{\bm{\eta}}_b(B,\bm\pi)\|_0}\right], \label{eqn:mdl_break}
\end{align}
where ${c}_b(\bm\pi)$ is the cardinality of $\{ i: {\pi}_{b-1}\le t_i< {\pi}_{b}\}$ for $b=1,\dots,{B}$, i.e.,
${c_b(\bm\pi)}$ is the number of time bins in time interval $b$.

\ignore{
\subsection{Statistical Approaches to Change Point Selection}\label{sec:mdlref}
In addition to the MDL principle, other statistical approaches have also been proposed for change point estimation.  For example, \citet{Yao88} studied the use of the Schwartz criterion and showed that it is consistent.  This Schwartz criterion can be interpreted as the classical Bayesian information criterion (BIC) if the location of a change point is treated as a free parameter (or a dimension of the model).  However, as for change point problems the likelihood functions violate the regularity conditions behind the classical BIC, this Schwartz criterion can likely be improved.  To circumvent this issue, \citet{Zhang-Siegmund07} provided an ingenious method to approximate the Bayes factor and developed a modified BIC (mBIC) for change point estimation.  This mBIC has been shown to be consistent for various settings.  More recently, \citet{Ninomiya14} derived a specialized Akaike information criterion (AIC) for change point estimation.  However, results from \citet[][Section~3.3]{Aue-Lee11} can be used to show that this AIC method is not consistent.

As mentioned before, MDL is a general principle for solving model selection problems.  One of its earliest applications to change point estimation is the image segmentation work of \citet{Leclerc89}, where the corresponding MDL criterion can be easily reduced to the classical one-dimensional change point problem setting.  See also \citet{Lee98:segcor,Lee00:segind} and \citet{Aue-Lee11}.  Other successful MDL applications to change point problems include discontinuity detection in nonparametric regression \citep{Lee02}, structural break detection in nonstationary time series models \citep{Davis-et-al06,Davis-et-al08,Yau-et-al14}, and change point detection in quantile modeling \citep{Aue-et-al14a, Aue-et-al14b}.  Under some mild regularity assumptions, these MDL solutions are shown to be consistent.

Lastly, for many change point problems simpler than the one that we consider in this paper, we found that the major penalty terms of mBIC and MDL are remarkably similar; e.g., Equation~(8) of \citet{Zhang-Siegmund07} and Equation~(6.2) of \citet{Lee97c}.  This provides reassuring evidence that both mBIC and MDL are reliable methods for solving change point detection problems.
}

\subsection{Practical minimization}
\label{sec:fullfit}
We now consider  minimization of the MDL criterion given in~(\ref{eqn:mdl_break}).
Once the number of time intervals, $B$, and the locations of change points, $\bm{\pi}$, are specified, unique estimates of $\bm{\gamma}$ and $\bm{\rho}$ can be obtained using the method described in Section~\ref{sec:hom}.
Since the time intervals are independent the required computation is highly parallelizable.  Nonetheless, minimization of~(\ref{eqn:mdl_break}) is not trivial. It involves an expensive combinatoric optimization: for each $B$, there are ${J-1}\choose{B-1}$ possibilities for $\bm{\pi}$. (A small number of these possibilities is not considered by our optimiziation procedure because we impose a minimum width on each time interval.)
To simplify computation, we propose a fast stepwise forward algorithm to approximately minimize~(\ref{eqn:mdl_break}), rather than attempting to globally minimize it.  Briefly, at each step we add a single change point to the current collection by selecting the new change point that decreases $\mathsf{mdl}_{\rm full}(B,\bm{\pi})$ the most.

\begin{algorithm}[t]
  \caption{A stepwise search algorithm for minimizing~(\ref{eqn:mdl_break})}\label{alg:treegrow}
  \begin{flushleft}
  \textbf{Input:} Wavelength-cross-time bin counts $\{(Y_{ij})\}_{i=1,\dots,N, j=1,\dots, J}$\\
  \textbf{Description:} A fast algorithm for approximately minimizing $\mathsf{mdl}_{\rm full} (B, \bm{\pi}, \bm{\gamma}, \bm{\rho})$
\end{flushleft}
  \begin{algorithmic}[1]
    \State $B \leftarrow 1$ and $\bm{\pi} \leftarrow ()$ (null vector)
        \State $z \leftarrow {\mathsf{mdl}}_{\rm full} (1, \bm{\pi})$.
          \State $z' \leftarrow z + 1$
        \While  {$z < z'$}
        \State $B \leftarrow B+1$
        \State $z' \leftarrow z$
        \State $\bm{\pi}' \leftarrow \bm{\pi}$
        \State Use a grid search to find the change point that when added to the current collection, $\bm{\pi}'$, gives the largest reduction in ${\mathsf{mdl}}_{\rm full}$. Denote this new set of change points as $\bm{\pi}$.
        If there is no possible change point to add, break the while-loop.
        \State $z \leftarrow {\mathsf{mdl}}_{\rm full}(B, \bm{\pi})$
        \EndWhile
        \State {\bf end while}
    \State {Return $\{B-1, \bm{\pi}'\}$}
\end{algorithmic}
\end{algorithm}

More precisely, the algorithm first calculates (\ref{eqn:mdl_break}) for the model with no change points (i.e., only one time interval), that is, it computes $\mathsf{mdl}_1 =  \mathsf{mdl}_{\rm full}(1,\bm{\pi})$. (This does not depend on $\bm\pi$ because there are no change points.)  In the second step, the algorithm considers all of the possible models with a single change point and selects the one that minimizes $\mathsf{mdl}_{\rm full}(2,\bm{\pi})$.
Calling the obtained  minimum $\mathsf{mdl}_2$, the change point is added to the model if $\mathsf{mdl}_2 < \mathsf{mdl}_1$; otherwise the algorithm terminates.  The algorithm continues in this way, attempting to add an additional change point to the model at each step by minimizing (\ref{eqn:mdl_break}). It stops when there is no possible change point remaining that decreases $\mathsf{mdl}_{\rm full}(B,\bm{\pi})$ and this current model is taken as the (approximate) minimizer of~(\ref{eqn:mdl_break}).  Although this greedy algorithm may result in local minimization,  it is fast and provides a good compromise between computational efficiency and accuracy \citep{Lee02:genetic}.  The exact steps of the algorithm appear  in Algorithm~\ref{alg:treegrow}.

\subsection{A Monte Carlo Procedure for Testing for the Existence of Change Points}
\label{sec:MCtest}

\begin{algorithm}[t]
  \caption{Monte Carlo significance test}\label{alg:mctest}
  \begin{flushleft}
  {\textbf{Input:} Working data $\{(w_i, t_j, Y_{ij})\}_{i=1,\dots,N, j=1,\dots, J}$, number of simulation $N_{\mathrm{sim}}$, significance level $\alpha$}\\
  {\textbf{Description:} Monte Carlo test for change point}
\end{flushleft}
  \begin{algorithmic}[1]
  \State Fit the time homogenous Poisson moded (with no change point, Section \ref{sec:hom}) and full Poisson model (including change points, Section~\ref{sec:fullfit}) to $\{(w_i, t_j, Y_{ij})\}_{i=1,\dots,N, j=1,\dots,
    J}$ and denote the test statistic
  $m^\star = \mathsf{mdl}_{\rm full}(\hat B, \bm{\hat\pi})-\mathsf{mdl}_{\rm
  full}(1, \bm{\pi}_{\rm null})$. ($\bm{\pi}_{\rm null} = ()$, null vector)
    \For {$k=1$ to $N_{\mathrm{sim}}$}
    \State Generate the data set $\mathcal D_k = \{(w_i, t_j, Y_{i p_j})\}$ where $\{p_1,\dots,p_J\}$ is a uniformly random permutation of $\{1,2,\dots, J\}$.
      \State Fit the full Poisson model (including change points) using the procedure in Section~\ref{sec:fullfit} to  $\mathcal D_k$.
      \State Denote the resulting test statistic by $m^\star_k$.
    \EndFor
  \State Estimate the $p$-value as
    \[
      \hat{p} = \frac{\#\{m_k^\star \ge m^\star\} + 1}{N_{\mathrm{sim}}+1}
    \]
  \State Conclude that there is at least one change point if $\hat{p} < \alpha$ for some significance level, $\alpha$.
\end{algorithmic}
\end{algorithm}

In practice it is useful to quantify evidence for existence of the change points.
For a test statistic, we use the change in the
MDL function (\ref{eqn:mdl_break})
due to introduction of change points,
$m^\star = \mathsf{mdl}_{\rm full}(\hat B, \bm{\hat\pi})-\mathsf{mdl}_{\rm
  full}(1, \bm{\pi}_{\rm null})$,
where
$\hat B$ and $\bm{\hat\pi}$ are computed by minimizing (\ref{eqn:mdl_break}),
and $\bm{\pi}_{\rm null}$ is the null vector containing no change point.
Here $\mathsf{mdl}_{\rm full}(1,\bm{\pi}_{\rm null})$ is the value of the MDL
function (\ref{eqn:mdl_break}) with no change points.
We approximate the null distribution of $m^\star$ via Monte Carlo using permutation.
Specifically, we generate $N_{\rm sim}$ uniformly random permutations of $\{1,2,\dots, J\}$
independently and the corresponding replicate null
datasets, $(\mathcal D_1, \ldots, \mathcal D_{N_{\rm sim}})$,
by shuffling the time indices of the original dataset according to these permutations (see Step 3 of Algorithm~\ref{alg:mctest}).
Each replicate dataset is then fit with the full model (including change
points) according to the MDL criterion specified in (\ref{eqn:mdl_break}) and
the resulting test statistic, $m^\star_i$ is computed (for $i=1,\ldots, N_{\rm
  sim})$.
Note that the MDL criterion (\ref{eqn:mdl_break}) with no change points,
$\mathsf{mdl}_{\rm full}(1,\bm{\pi}_{\rm null})$,
is invariant to the permutation and thus need not to be re-computed for the permuted datasets
$(\mathcal D_1, \ldots, \mathcal D_{N_{\rm sim}})$.
The empirical distribution of $(m^\star_1, \ldots, m^\star_{N_{\rm
    sim}})$ is a Monte Carlo approximation to the null distribution of
$m^\star$ and a p-value can be computed by comparing $m^\star$ to $(m^\star_1,
\ldots, m^\star_{N_{\rm sim}})$. Details appear in Algorithm~\ref{alg:mctest}.

\section{Numerical Experiments}
\label{sec:sim}
We conducted a series of simulation experiments  to evaluate the empirical properties of the proposed methodology.
Binned photon counts were simulated under \changes{eight} different spectral-temporal test functions, $\lambda(\tm, \wl)$.  In order to make the experiments as realistic as possible, \changes{seven of the test functions} were constructed by fitting a spectral-temporal model, using the method described in Section \ref{sec:break}, to the seven observations of FK\,Com\  described in Section~\ref{sec:real}.
\changes{Because these test functions have at most two change points and because we would like to investigate the statistical properties of our method with more change points, an eighth test function was created by concatenating two of the original test functions. (I particular, those created from
observation identification (ObsID) numbers 13251 and 12298, see Table~\ref{tab:tfdetails}.)}
\changes{The analyses carried out here, including the binning of the data, is identical to that described in Section~\ref{sec:real}.}
\changes{Further details are summarized in Table~\ref{tab:tfdetails}.}
Together \changes{these eight test functions} cover a wide range of possible scenarios, ranging from models with no change point to models with \changes{four} change points.  For each of the \changes{eight} $\lambda(\tm, \wl)$, 200 data sets were simulated according to~(\ref{eq:pois-mdl}), using the real data values for $\delta_t$, $\delta_w$, and $A(w_i)=\sum^K_{k=1}A_k(w_i)$.  The proposed methodology was applied to all $\changes{8}\times 200 =1400$ simulated data sets  to estimate $\lambda(\tm, \wl)$ for each.

\begin{table}[t]
\begin{center}
  \changes{\caption{Details of simulation settings.
The columns report (1) the
chronologically ordered (except the last row)
ObsID numbers of the
datasets used to generate each test function,
(2) the number of wavelength bins,
(3) the number of time bins,
(4) the number of basis functions $P$, as in Equation (\ref{eqn:gensemi}),
used to fit the corresponding FK Com observation,
(5) the number of time intervals, i.e., the number of change points plus one,
estimated from the corresponding FK Com observation,
(6) vector of the numbers of emission lines with element $j$ corresponding
to time interval $j$,
estimated from the corresponding FK Com observation.
}
\label{tab:tfdetails}
\begin{tabular}{cccccc}
  \hline
ObsID & $N$ & $J$ & $P$ & $B$ & $M$ \\
  \hline
 12297 & 142 & 21 & 34 &  2 & (0, 1)  \\
 12356 & 142 & 32 & 34  & 1  & (1)\\
 13251 & 142 & 49 & 34  & 3 & (0, 0, 0)\\
 12298 & 142 & 20 & 34  & 2 & (2, 1)\\
 13259 & 142 & 23 & 34  & 2 & (0, 0)\\
 12357 & 142 & 18 & 34 & 1  & (1)\\
 12299 & 142 & 20 & 34  & 3 & (0, 0, 0)\\
 13251, 12298 & 142 & 69 & 34  & 5 & (0, 0, 0, 2, 1)\\
   \hline
\end{tabular}}
\end{center}
\end{table}

The simulation results are summarized in Table~\ref{tab:sim}.  The first column lists \changes{the ObsID number} for the observations of FK\,Com\ (in chronological rather than ascending order) used to generate the $\lambda(\tm, \wl)$, see Section~\ref{sec:real}.  The second column shows the (true) numbers of homogeneous time intervals, $B$, in each of the $\lambda(\tm, \wl)$; recall the number of change points is $B-1$. The third column reports the percentage of the simulated datasets for which the fitted $\hat B$ exactly equals the true value. The proposed methodology achieved a correct recovery rate of at least 94\% in 6 out of the \changes{8} simulations settings.

\begin{table}[t]
\begin{center}
\caption{Simulation results. The columns report (1) the chronologically ordered \changes{(except the last row)} ObsID numbers
 of the datasets
used to generate each test function,
(2) the true number of time intervals, i.e., the true number of change points plus one,
(3) the percentage of the fitted $\hat{B}$ that exactly equal the true $B$, and
(4)  the root mean square error of $\hat{\bm{\pi}}$ among those simulations with $\hat{B}=B$.}
\label{tab:sim}
\begin{tabular}{cccccc}
  \hline
ObsID & $B$ & \% correct $\hat{B}$ & RMSE($\hat{\bm{\pi}}$)\\% & RMSE($\hat{f}$) \\
  \hline
 12297 & 2 & 94 & 0.70\\% & 1.6452e-05 \\
 12356 & 1 & 100 & 0\\% & 4.8220e-06 \\
 13251 & 3 & 98 & 0.66\\% & 8.3488e-06 \\
 12298 & 2 & 98 & 0.29\\% & 9.2506e-06 \\
 13259 & 2 & 80 & 0.69\\% & 1.4185e-05 \\
 12357 & 1 & 100 & 0\\% & 5.5745e-06 \\
 12299 & 3 & 100 & 0\\% & 1.2936e-05 \\
 \changes{13251, 12298} & \changes{5} & \changes{85} & \changes{0.66}\\
   \hline
\end{tabular}
\end{center}
\end{table}

Our method also accurately estimates the change points, $\bm{\pi}$, \changes{as shown in Column~4 of Table~\ref{tab:sim}. For those simulations where $B$ is correctly estimated,
Column~4 gives the root mean squared error (RMSE) of $\hat{\bm{\pi}}$, in terms of the time bin width ($\delta_t =2000$ seconds in all cases).}   The RMSE is less than one bin width in all \changes{eight} cases.  Figure~\ref{fig:sim} displays histograms of the estimated change point location  among those simulations for which $\hat B=B$ under the \changes{five} test functions with less than a 100\% recovery rate  for $\bm{\pi}$. In  all \changes{five} cases, the estimated change point locations are narrowly centered around the true change points.
Finally, for those simulations for which $\hat{{B}}={B}$ and $\hat{\bm{\pi}}=\bm{\pi}$,
the uncertainty of $\hat{f}$ is summarized
in Figure~\ref{fig:cont}.
Here we only present results for one test function (ObsID~12299); plots for the other test functions appear in
Appendix~\ref{app:sim}.  For each of the $B=3$ homogeneous time intervals associated with ObsID~12299,
we overlay the fitted  $\hat{f}$ of the simulations  for which
$\hat{{B}}={B}$ and $\hat{\bm{\pi}}=\bm{\pi}$. The results are plotted in panels (a)-(c) of Figure \ref{fig:cont}. The relatively high uncertainty at high wavelengths is partly due to the small effective area at these wavelengths; see panel (d) of Figure \ref{fig:cont}.

\begin{figure}[t]
  \centering
\includegraphics[width=0.8\linewidth]{barplots2-1.pdf}
\includegraphics[width=0.8\linewidth]{barplots2-3.pdf}
\caption{Histograms of detected change point locations in the simulation study.
  The title of each panel refers to the ObsIDs  \changes{given in Tables~\ref{tab:tfdetails} and \ref{tab:sim}. The widths of the histogram bins is the same as the time binning of the data;
  bins containing} the true change-points locations are represented in black.
  The horizontal axis \changes{is elapsed spacecraft time and the panel corresponding to the concatenation of ObsID~13251 and ObsID~12298 has a broken x-axis because there is a gap in spacecraft time between the two observations.}}\label{fig:sim}
\end{figure}

\begin{figure}[p]
  \centering
\includegraphics[width=0.49\linewidth]{12299-cont-a.pdf}
\includegraphics[width=0.49\linewidth]{12299-cont-b.pdf}\\
\includegraphics[width=0.49\linewidth]{12299-cont-c.pdf}
\includegraphics[width=0.49\linewidth]{effect-d.pdf}
\caption{
  \changes{Panels (a) - (c):  Overlay plots of the $\hat{f}$ fitted to the 200 simulations generated with the test function corresponding to ObsID~12299.
  Notice that all 200 simulations resulted in $\hat{B}=B$ and $\hat{\bm{\pi}}=\bm{\pi}$.
  The three panels represent the $B=3$ homogeneous time intervals of ObsID~12299.
  The areas between the $\hat{f}$ computed with each
  simulated dataset and the underlying test function are colored and overlaid.
  Darker color represents more overlap of these areas  and thus summarizes the
  concentration of the
  the fitted spectra  around the underlying test function.
  The solid black line represents the
  smooth component, $f$, of test function 12299 and the dashed lines represent
  the 5\% and 95\% (pointwise) percentiles of $\hat{f}$.

  Similar plots for the other test functions appear in Appendix~\ref{app:sim}.
  Panel (d): Plot of the effective area $\sum^K_{k=1} A_k(w)$ curve for test function 12299.}
}\label{fig:cont}
\end{figure}

\section{Application to X-ray Observations of FK\,Com}
\label{sec:real}

FK\,Com\ is an evolved yellow giant star (spectral type G4\,{\tt III}) with a
mass of about  $8.4$ times that of the sun and at a distance of
$550\pm60$~light-years. Typically, late-spectral-type giants are slow rotators,
and thus do not maintain a strong magnetic dynamo that can sustain a
high-temperature corona. However, FK\,Com\ is an extremely fast rotator, with
equatorial velocities of 179~km/s, and a rotational period of 2.4~days
\citep{korh:etal:99,stras:09}.  (For comparison, the Sun rotates at about 2 km/s at its equator.)  Giant stars are formed when a sun-like star
depletes the hydrogen at its core and thus transitions from producing energy
though the fusion of hydrogen into helium to the fusion of heavier elements.
Typical sun-like progenitor stars would require rotational speeds that would
tear them apart in order to form a giant star with such a rapid rate of spin.
Thus, the favored hypothesis for the formation of FK\,Com-like stars is that they are the result of a binary star system that has coalesced through a process of mass transfer, dramatically speeding up the stellar rotation  in the process \citep{bopp:sten:81}.
If this is the case, the high rotation should induce a strong dynamo and consequently a dynamic corona, filled with high-temperature plasma and organized into loop-like structures via magnetic fields.

FK\,Com\ exhibits long duration spots---akin to sunspots---and flares associated with these spots. Optical Doppler imaging observations of FK\,Com\ demonstrate the existence of persistent spot-like features \citep{Elstner-Korhonen05} that remain localized to specific longitudes for long timescales ($>$months).
Analysis of an X-ray flare observed with the XMM-{\sl Newton} observatory showed that the flare is tied to these spots \citep{Drake-Chung-Kashyap08}.
Numerical MHD modeling of the stellar wind and the magnetic structure \citep{Cohen-Drake-Kashyap10} shows a highly complex coronal environment, with the most notable feature being even putatively stable coronal loops ``wrapping around'' due to the rotation.
Thus, even though the corona appears to be solar-like in its origin, it is an extreme case of such coronae, and is therefore of considerable interest to both theoretical and observational studies.

A long-duration, high-resolution X-ray observation of FK\,Com\ was carried out with the {\sl Chandra X-ray Observatory}  \changes{over seven discontinuous time periods} in April 2011.
The data show that the source brightness is constantly fluctuating with occasional periods of singular flaring activity, see Figures~\ref{fig:fkcom_lc} and \ref{fig:data}.
As a test case for our method,  \changes{we apply it to these seven {\sl Chandra} observations of FK\,Com.} Each was preprocessed by discarding photons with
 wavelength less than 1.65  or greater than 31~\AA ngstr\"{o}ms.
\changes{The natural resolution of {\sl Chandra}  data of this sort is extremely high, resulting in a very large, but very sparse table of photon counts; typically there is at most one photon recorded in each wavelength-cross-time bin. We can obtain substantial computational savings with only a modest sacrifice of scientific information by reducing the resolution of the data. Massive events in the corona of FK\,Com, for example, require appreciable time and are associated with dramatic spectral shifts. Change points typically appear with frequencies on the order of tens of thousands of seconds, whereas the inherent bin width is about 3 seconds. In our analyses, we reduce this resolution and use temporal bins of width $\delta_t=2000$~seconds. Similarly, the inherent wavelength bin width is 0.005~\AA ngstr\"{o}ms\ and we set $\delta_w=0.2$~\AA ngstr\"{o}ms. Though a substantial reduction in spectral resolution, our analysis is much higher resolution than the standard hardness ratio analysis, which can be viewed as using two or three wavelength bins of width 10-15~\AA ngstr\"{o}ms.}
Our fitted $\lambda(t_j,w_i)$ are plotted in Figure~\ref{fkcoms2}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.8\linewidth]{fkcom_phase_time_lc.pdf}
\caption{The phase-time light curve of FK\,Com.
The source brightness (i.e., light curve) over the course of the {\sl Chandra} observation is shown during periods when the source was observed.
The horizontal axis represents spacecraft elapsed time, and the vertical axis represents the rotational phase of the star.
The tilted red lines show the phase-time correspondence for each observation, showing the time and rotational phase coverage of the full dataset.
Overlaid on the phase-time line, displayed as a deviation from it, are sparklines depicting the   smoothed light curves (green curves).
The segments are labeled by the ObsID numbers.
This representation of the light curve shows that there is no obvious rotational signature in the light curve, and the variability is stochastic.
Notice the large flare that occurs in ObsID~12299: our analysis shows that the spectrum varies significantly during this flare, see Figure~\ref{fkcoms2}.
}
\label{fig:fkcom_lc}
\end{center}
\end{figure}

\begin{table}[H]
\begin{center}
\caption{Monte Carlo p-values for the test for change points in each of the seven observations of FK\,Com}
\label{tab:pvalue}
\begin{tabular}{cccccccc}
  \hline
ObsID & 12297 & 12356 & 13251 & 12298 & 13259 & 12357 & 12299\\
p-value & 0.12 & 1.00 & $<$0.01 & 0.01 & 0.14 & 1.00 & $<$0.01\\
   \hline
\end{tabular}
\end{center}

\end{table}

We also applied the Monte Carlo test for change points as described in Section~\ref{sec:MCtest}; the resulting p-values for the seven datasets appear in Table~\ref{tab:pvalue}.
There is no evidence for change points in two of the observations (ObsID~12356
and 12357), very weak evidence for two (ObsID~12297 and 13259), and strong
evidence for the other three. This agrees with the fitted models illustrated in
Figure~\ref{fkcoms2}, but quantifies the degree of evidence for the fitted
change points.  Finally, the two test functions that resulted in relatively low
recovery rates of $B$ in the simulations described in Section \ref{sec:sim}
correspond to the two data sets with the weakest evidence for the a change
point, compare ObsID 12297 and 13259 in Table~\ref{tab:sim} (\% correct $\hat{B}$)
with Table~\ref{tab:pvalue} (p-value).

\begin{figure}[t]
\begin{center}
\includegraphics[width=\linewidth]{fkcom_timewvl_scatterplot_wgap.png}
\caption{\changes{Observed data. Each point represents a photon detected with HETGS+ACIS-S configuration of {\it Chandra} (including both low- and high-energy
gratings data).  The wavelengths are given in \AA ngstr\"{o}ms, and the time
is in spacecraft clock seconds.  The ObsIDs of each segment
are marked to the right of the red vertical line representing the
starting time of the observation segment.  The collected data are sparse, but some
changes in the density of points can be discerned even by eye,
indicating temporal and spectral variations.}}
\label{fig:data}
\end{center}
\end{figure}

\changes{Although we find two or fewer change points in each observation, the number
was not limited in our analyses.
Operational constraints limit the continuous observation time available with
the space-based {\it Chandra X-ray Observatory}.
The change points that we observe stem from massive energetic shifts in the cornea of FK\,Com\ and take time to evolve. Keeping in mind that a single massive shift of this sort on the Sun would likely destroy human civilization, observing a total of seven in one (discontinuous) observation period spanning less than a fortnight (see Figure~\ref{fig:fkcom_lc}) is not a small number.
}

\begin{figure}[p!]
  \centering
\includegraphics[width=0.9\linewidth]{real2.pdf}
\caption{Fitted $\lambda(\tm, \wl)$ for seven {\sl Chandra} observations of FK\,Com. The ObsID of each dataset appears in the title of each panel. The left plot in each row shows  the fitted spectra for each identified time interval; the solid black, dashed red, and dotted green curves correspond to the first, second, and third time intervals, respectively.    The right plot in each row is a heat maps that represents the best fit values of $\lambda(\tm, \wl)$. Time is given in elapsed spacecraft time.
\label{fkcoms2}
}
\end{figure}

This analysis indicates that large changes in intensities are accompanied by significant changes in the spectrum of the source.  This is readily apparent in the results for ObsID 12299, which shows the observation separated into three segments with different intensities and spectral shapes; see ObsID~12299 in Figure~\ref{fig:fkcom_lc} and row 7 of Figure~\ref{fkcoms2}. The spectrum changes to being dominated by high-energy photons due to an increase in high-temperature plasma when the flare is set off.  As the flare decays, even when the observed intensity is relatively high, the spectrum can be seen returning to the pre-flare state, though with enhanced emission over the line-dominated region around $15$~\AA ngstr\"{o}ms.

From a theoretical point of view, this is not an obvious result, since for active stars like FK\,Com, it is believed that even the apparently quiescent corona emits X-rays as a superposition of a large number of weak flares \citep[see, e.g.,][]{Kashyap-Drake-Gudel02, Gudel04}.
However, our analysis favors the idea that hot plasma could become trapped in stretched field lines wrapping around the star, and thus cool radiatively, and flares are only triggered sporadically when the magnetic strain becomes too large to be sustained.
While a detailed modeling and interpretation of this picture is beyond the scope of this paper,  we point out that even a simple application of a statistically principled method is enough to obtain an important result: The quiescent corona of FK\,Com\ is {\sl unlike} the flaring corona, and in this sense is much more solar-like than previously anticipated.

\section{Discussion}

\label{sec:conc}
This article develops a novel approach for modeling change points in a class of marked Poisson processes, specifically the time-varying spectra of high-energy astrophysical sources.  The approach includes a family of change point models, a practical algorithm for estimation, and an MDL criterion for selecting the best fitting model. Despite numerous reported applications where the MDL principle leads to reliable model selection criteria, its use in the statistical literature remains rather sparse.  It is our hope that the present work, in which the MDL principle is applied to solve a complicated model selection problem that involves the ``large $p$ small $n$'' scenario, will encourage and stimulate renewed research interests in this useful principle among statisticians.

Although the methods developed in this article are motivated by the search for change points in high-energy time-varying spectra, it is possible to extend them to include spatial variations as well. Typically, astrophysical images include sources at many different scales, ranging from unresolved point sources to complex structures that may span the field of view. Owing to instrumental point-spread functions, point sources in close proximity may overlap and in many cases are overlaid on larger scale structures.  An extension
that can account for change points in the spatial plane or in the spatial cross spectral hyperplane would be quite useful in astronomical data analysis.  Statistically, this involves extending the dimension of the marks in the Poisson process from the one dimensional photon wavelengths to two dimension spatial locations or three dimensional spatial cross spectral measurements. This would, of course, involve extensions of the non-parametric models developed here for spectra to higher dimensional models for images or joint models for images and spectra.

\subsubsection*{Acknowledgments}

\changes{The authors are grateful to the Referees, the Associate Editor and the Editor, Professor Nicoleta Serban, for their many useful and constructive comments. They substantially improved the paper. This work was conducted under the auspices of the CHASC International Astrostatistics Center.  CHASC is supported by NSF grants DMS-12-08791, DMS-12-09232 and DMS-1513484.  In addition, Vinay Kashyap acknowledges support from Chandra grant GO1-12021X, and from NASA's contract to the Chandra X-ray Center, NAS8-03060, Thomas Lee from National Science Foundation grants DMS 12-09226 and DMS 15-12945, David van Dyk from a Wolfson Research Merit Award provided by the British Royal Society and from a Marie-Curie Career Integration Grant provided by the European Commission.
We also thank Xiao-Li Meng, Tom Ayres, Ofer Cohen, Jeremy Drake,
David Garcia-Alvarez, Dave Huenemoerder, Heidi Korhonen, Steve Saar,
Paola Testa, and Aad van Ballegooijen for many useful discussions.}

\bigskip