---
abstract: |
    We develop a nonlinear semi-parametric Gaussian process model to estimate periods of Miras with sparsely sampled light curves. The model uses a sinusoidal basis for the periodic variation and a Gaussian process for the stochastic changes. We use maximum likelihood to estimate the period and the parameters of the Gaussian process, while integrating out the effects of other nuisance parameters in the model with respect to a suitable prior distribution obtained from earlier studies. Since the likelihood is highly multimodal for period, we implement a hybrid method that applies the quasi-Newton algorithm for Gaussian process parameters and search the period/frequency parameter space over a dense grid.

    A large-scale, high-fidelity simulation is conducted to mimic the sampling quality of Mira light curves obtained by the M33 Synoptic Stellar Survey. The simulated data set is publicly available and can serve as a testbed for future evaluation of different period estimation methods. The semi-parametric model outperforms an existing algorithm on this simulated test data set as measured by period recovery rate and quality of the resulting Period-Luminosity relations.
author:
- 'Shiyuan He, Wenlong Yuan, Jianhua Z. Huang, James Long & Lucas M. Macri'
bibliography:
- 'm33gp.bib'
title: |
    Period estimation for sparsely sampled quasi-periodic\
    light curves applied to Miras
---

Introduction
============

The determination of reliable periods for variable stars has been an area of interest in astronomy for at least four centuries, since the discovery of the variability of Mira ($o$ Ceti) by Fabricius in 1596 and the first attempts to determine its period by Holwarda & Bouillaud in the mid-1600s. The availability of electronic computers for astronomical research half a century ago enabled the development of many algorithms to estimate periods quickly and reliably, such as @Lafler1965 [@Lomb1976; @Scargle1982].

The aforementioned algorithms work best in the case of periodic variations with constant amplitude and Mira variables present several challenges in this regard. While their periods of pulsation are stable except for a few intriguing cases [@Templeton2005], Mira light curves can exhibit widely varying amplitudes from cycle to cycle [see, for example, the historical light curve of Mira compiled by @Templeton2009]. In the case of C-rich Miras, the stochastic changes in mean magnitude across cycles [e.g., @Marsakova1999] only complicate the problem further. The wide variety of light curves for long-period variables, already recognized by @Campbell1925 and @Ludendorff1928, may complicate the identification of Miras among other stars. Lastly, from a purely practical standpoint, it is simpler to obtain light curves spanning several cycles for RR Lyraes or Cepheids (with periods ranging from $\sim 0.5$ to $\sim 100$ d) than for Miras (with periods ranging from $\sim 100$ to $\sim 1500$ d).

Despite these challenges, the identification and determination of robust periods for Miras — especially in the regime of sparsely sampled, low signal-to-noise light curves — would be very beneficial for the determination of distances to galaxies of any type. Thanks to the unprecedented temporal coverage of the Large Magellanic Cloud (LMC) by microlensing surveys, the availability of large samples of extremely well-observed Miras has led to a thorough characterization of their period-luminosity relations at various wavelengths [@Wood1999; @Ita2004; @Soszynski2007]. The dispersion of the $K$-band period-luminosity relation [@Glass2003 $\sigma=0.13$ mag], is quite comparable to that of Cepheids at the same wavelength [@Macri2015 $\sigma=0.09$ mag] and makes them competitive distance indicators.

The third phase of the OGLE survey [@Udalski2008] imaged most of the LMC with little interruption over 7.5 years and resulted in the discovery of 1663 Miras [@Soszynski2009] with a median of 466 photometric measurements per object. The temporal sampling of these light curves and their photometric precision are exceptional relative to typical astronomical surveys and make period estimation relatively easy. In comparison, a similar span of observations of M33 by the DIRECT [@Macri2001] and M33SSS projects [@Pellerin2011] in the $I$-band consists of a median number of 44 somewhat noisy measurements, heavily concentrated in a few observing seasons. Representative Mira light curves from the OGLE & DIRECT/M33SSS surveys are shown in Fig. \[fig:example.mira.lc\]. There are several reasons for the striking difference in quality between these two data sets. The LMC Miras are among the brightest objects in the OGLE fields, whereas their M33 counterparts are among the faintest in the aforementioned surveys of this galaxy. While the effective exposure times of all these surveys are quite comparable, after taking into account differences in collecting area of their respective telescopes, M33 lies approximately 6.2 mag farther in terms of its $I$-band apparent distance modulus. Furthermore, the main goal of the OGLE project (detection of microlensing events) requires a very dense temporal sampling of the survey fields; this is achieved by using a dedicated telescope and is helped by the fact that the LMC is observable nearly all year long from the site. In contrast, the observations of M33 were carried out using shared facilities (available only a few nights per month) with the primary purpose of studying Cepheids and eclipsing binaries (which do not require exceptionally dense temporal sampling), and the galaxy is only observable all night long for $\sim 1/3$ of the year. Standard period estimation algorithms, which work well for high signal-to-noise, well sampled light curves such as those obtained by OGLE, will fail on more typical data sets represented by the M33 observations. The purpose of this work is to develop and test a methodology for estimating periods for sparsely sampled, noisy, quasi-periodic light curves such as those of Miras observed in M33 by the aforementioned projects.

\[fig:example.mira.lc\] ![Representative Mira light curves observed by OGLE-III in the Large Magellanic Cloud (top) and DIRECT/M33SSS in M33 (bottom).](fig01a.eps "fig:"){width="49.00000%"} ![Representative Mira light curves observed by OGLE-III in the Large Magellanic Cloud (top) and DIRECT/M33SSS in M33 (bottom).](fig01b.eps "fig:"){width="49.00000%"}

The rest of the paper is organized as follows. In §\[sec.background\] we review several existing period estimation methods. In §\[sec.model\] we introduce a new semi-parametric (SP) model for Mira variables which uses a Gaussian process to account for deviations from strict periodicity. We use maximum likelihood to estimate the period and the parameters of the Gaussian process, while other nuisance parameters in the model are integrated out with respect to some prior distributions using earlier studies. Since the likelihood is highly multimodal for the period/frequency parameter, we implement a hybrid method that applies the quasi-Newton algorithm for Gaussian process parameters and a grid search for the period/frequency parameter. In order to assess the effectiveness of the SP model, in §\[sec.construct.test\] we carefully construct a simulated data set by fitting smooth functions to the light curves of well-observed OGLE LMC Miras and resampling them at the cadence, noise level, and completeness limits of the aforementioned M33 observations. Using the simulated data, in §\[sec.evaluation\] we compare the performance of existing period estimation methods to our SP model. We find that our proposed model shows an improvement over the generalized Lomb-Scargle (GLS) model under various metrics. In §\[sec.discussion\], we conclude and discuss some future applications. Simulated light curves for reproducing the results in the paper and performance benchmarking are made publicly available as supplementary material.

Period estimation techniques {#sec.background}
============================

Let $y_i$ be the magnitude of a variable star observed at time $t_i$ (in units of days) with uncertainty $\sigma_i$. The data set for this object, obtained as part of a time-series survey with $n$ epochs is $\{(t_i,y_i,\sigma_i)\}_{i=1}^n$. One common approach to estimate the primary frequency of such an object is to assume some parametric model for brightness variation and then use maximum likelihood to estimate parameters. @Zechmeister2009 define the GLS model as $$\label{eq:gls}
y_i = m + a\sin(2\pi f t_i + \phi) + \sigma_i\epsilon_i,$$ where $\epsilon_i \sim \mathcal{N}(0,1)$, $m$ is the mean magnitude, $a$ is the amplitude, $\phi \in [-\pi,\pi]$ is the phase, and $f$ is the frequency [see @Reimann1994 for early work in this model]. Using the sine angle addition formula and letting $\beta_1 = a\cos(\phi)$ and $\beta_2 = a\sin(\phi)$ one obtains $$\label{eq:gls2}
y_i = m + \beta_1\sin(2\pi f t_i) + \beta_2\cos(2\pi f t_i) + \sigma_i\epsilon_i.$$ The likelihood function of this model is highly multimodal in $f$. However at a fixed $f$ the model is linear in the parameters $(m,\beta_1,\beta_2)$. These two facts motivate the computation strategy of performing a grid search across frequency and minimizing a weighted least squares $$\begin{split}
& (\widehat{m}(f),\widehat{\beta}_1(f),\widehat{\beta}_2(f) ) \\
& \qquad = \operatorname*{arg\,min}_{m,\beta_1,\beta_2} \sum_{i=1}^n \frac{1}{\sigma_i^2} \left\{y_i - m \right.\\
& \qquad \qquad \left. -\beta_1\sin(2\pi f t_i) - \beta_2\cos(2\pi f t_i)\right\}^2,
\end{split}$$ at every frequency $f$ on the grid. Under the normality assumption, the weighted least squares minimization is equivalent to maximizing the likelihood. Since the model is linear, computation of $\widehat{m}(f),\widehat{\beta}_1(f),\widehat{\beta}_2(f)$ is straightforward. The residual sums of squares at $f$ is $$\begin{split}
{\rm RSS}(f) & = \sum_{i=1}^n \frac{1}{\sigma_i^2} \{y_i - \widehat{m}(f) \\
& \quad - \widehat{\beta}_1(f)\sin(2\pi f t_i) - \widehat{\beta}_2(f)\cos(2\pi f t_i)\}^2,
\end{split}$$ and the maximum likelihood estimator for $f$ is $$\widehat{f} = \operatorname*{arg\,min}_{f} {\rm RSS}(f).$$ Define ${\rm RSS}_0$ as the (weighted) sum of squared residuals when fitting a model with only an intercept term $m$. The periodogram is defined as $$\label{eq:glsP}
S_{\rm LS}(f) = \frac{(n-3)({\rm RSS}_0 - {\rm RSS}(f))}{2{\rm RSS}(f)}.$$ The periodogram has the property that if the light curve of the star is white noise (i.e., $y_i = m + \epsilon_i$), $S_{\rm LS}(f)$ has an $F_{2,n-3}$ distribution. Thus the periodogram may be used for controlling the “false alarm probability,” the potential that a peak in the periodogram is due to noise [@Schwarzenberg1996].

A large number of period estimation algorithms in astronomy are closely related to GLS. The LS method is identical to GLS but first normalizes magnitudes to mean $0$ and does not fit the $m$ term [@Lomb1976; @Scargle1982]. The “harmonic analysis of variance” includes an arbitrary number of harmonics in Equation [@Quinn1991; @Schwarzenberg1996]. @Bretthorst2013 incorporates Bayesian priors on the parameters $\beta_1$ and $\beta_2$. The method is similar to performing a discrete Fourier transform and selecting the frequency which maximizes the @Deeming1975 periodogram. However, @Reimann1994 showed that GLS has better consistency properties than the Deeming periodogram.

![Light curve of a Mira in the LMC observed by OGLE (black points), decomposed following Eqn. \[eqn:basicDecomposition\]. Top panel: fitted light curve; middle panel: periodic signal, $m+q(t)$; bottom panel: stochastic variations, $m+h(t)$.[]{data-label="fig:onemira"}](fig02.eps){width="49.00000%"}

It is also possible to use non-sinusoidal models but compute and minimize the residual sum of squares as above. For example, @Hall2000 consider the Nadaraya-Watson estimator and @Reimann1994 uses the Supersmoother algorithm. @Wang2012 used Gaussian processes with a periodic kernel and found the period with maximum likelihood or minimum leave-one-out cross-validation error.

None of the above methods account for the non-periodic variation present in Miras. While these methods are adequate for densely sampled Mira light curves (where the quantity of data overwhelms model inadequacy), their performance deteriorates in the sparsely sampled regime. In Section \[sec.evaluation\], we compare our proposed model with the LS method.

The SP model {#sec.model}
============

Suppose the data $\{(t_i,y_i,\sigma_i)\}_{i=1}^n$ are modeled by $$y_i = g(t_i) + \sigma_i\epsilon_i\, ,$$ where $g(t_i)$ is the light curve signal and the $\epsilon_i\sim N(0,1)$ is independent of other $\epsilon_j$s. The signal of the light curve is further decomposed into three parts, $$\begin{split} 
g(t) & = m + q(t) + h(t)\\
     & = m + \beta_1\cos(2\pi ft) + \beta_2\sin(2\pi ft) + h(t)\, ,
\label{eqn:basicDecomposition}
\end{split}$$ where $m$ is the long-run average magnitude, $q(t) = \beta_1\cos(2\pi ft) + \beta_2\sin(2\pi ft)$ with frequency $f$ is the exactly periodic signal, and $h(t)$ is the stochastic deviation from a constant mean magnitude, caused by the formation and destruction of dust in the cool atmospheres of Miras. Fig. \[fig:onemira\] provides an example of the decomposition for a Mira light curve. The first two terms $m + q(t)$ in Eqn. \[eqn:basicDecomposition\] are exactly the same as the GLS model of Eqn. \[eq:gls2\]. To simplify notation, we define ${\mathbf{b}}_f(t) =(\cos(2\pi ft),\sin(2\pi ft))^T$, so that $q(t) = {\mathbf{b}}_f(t)^T{\boldsymbol{\beta}}$. The subscript in ${\mathbf{b}}_f(t)$ emphasizes that the basis is parameterized by the frequency $f$.

An SP statistical model is constructed in Eqn. \[eqn:basicDecomposition\] if we assume $h(t)$ is a smooth function that belongs to a reproducing kernel Hilbert space $\mathcal{H}$ with norm $\Vert\cdot\Vert_{\mathcal{H}}$ and a reproducing kernel $K(\cdot,\cdot)$. For this model, if the frequency $f$ is known, we obtain a least squares kernel machine considered in @Liu2007. Because the frequency is unknown, the response function is nonlinear in $f$. This nonlinearity and the multimodality in $f$ of the residual sum of squares provide additional challenges that require a novel solution.

Besides the additive formulation in Eqn. \[eqn:basicDecomposition\], another possible solution to account for the quasi-periodicity is a multiplicative model such as $g(t)  = m + h(t) q(t)$, where the amplitude of the strictly periodic term $q(t)$ is modified by a smooth function $h(t)$. However, the multiplicative model is more computationally intensive in nature and requires imposing a positive constraint on $h(t)$. As we will show in the following subsections, the $h(t)$ term in the additive model can be easily absorbed into the likelihood function. Nevertheless, the multiplicative approach is an interesting alternative approach to model formulation and is open to future study.

Equivalent formulations
-----------------------

Following §5.2 of @Rasmussen2005, for fixed $f$, the parameters $m,\beta_1,\beta_2$ and $h(t)$ in Eqn. \[eqn:basicDecomposition\] are jointly estimated by minimizing $$\begin{split}
& \sum_{i=1}^n \frac{1}{\sigma_i^2} [
y_i - m - \beta_1\cos(2\pi ft_i) \\
& \qquad - \beta_2\sin(2\pi ft_i) - h(t_i) ]^2 + \lambda \Vert h(\cdot)\Vert_{\mathcal{H}}^2,
\label{eqn:generalFormulation}
\end{split}$$ where $\lambda$ is a regularization parameter. A smoothing/penalized spline model for $h(t)$ is a special case of the general formulation of Eqn. \[eqn:generalFormulation\] with a specifically defined kernel; see §6.3 of @Rasmussen2005. For fixed $\lambda$, the solution of $h(t)$ is a linear combination of $n$ basis functions $K(t_i,t)$, $i=1,2,\cdots,n$, by the representer theorem [@Kimeldorf1971; @OSullivan1986]. It is still left for us to choose the regularization parameter $\lambda$ to balance data fitting and the smoothness of the function $h(t)$.

An equivalent point of view to the above regularization approach is to impose a Gaussian process prior on the function $h(t)$; see §5.2.3 of @Rasmussen2005. The benefit of this view is that it provides an automatic method for selecting the regularization parameter $\lambda$. In particular, we can absorb $\lambda$ into the definition of the norm $\Vert \cdot\Vert_{\mathcal{H}}$ and assume the term $h(t)$ in Eqn. \[eqn:basicDecomposition\] follows a Gaussian process, $h(t) \sim {\mathcal{GP}}(0,k_{{\boldsymbol{\theta}}}(t,t'))$, with the squared exponential kernel $k_{{\boldsymbol{\theta}}} (t,t') = \theta_1^2\exp \left(-\frac{(t-t')^2}{2\theta_2^2}\right),$ and parameters ${\boldsymbol{\theta}}= (\theta_1,\theta_2)$. The Gaussian process assumption implies that at any finite number of time points $t_1,t_2,\cdots,t_s$, the vector $(h(t_1),\cdots,h(t_s))$ is multivariate normally distributed, with zero mean and covariance matrix $\mathbf{K} = (k(t_i,t_j))$. This imposes a prior on the function space of $h(t)$. We also impose priors on $m$ and ${\boldsymbol{\beta}}$ in Eqn. \[eqn:basicDecomposition\]. In particular, we assume $m \sim {\mathcal{N}}(m_0, \sigma_m^2)$ and $ {\boldsymbol{\beta}}\sim {\mathcal{N}}({\mathbf{0}},\sigma_b^2{\mathbf{I}})$. The prior mean $m_0$ can be interpreted as the average magnitude of Miras in a certain galaxy, and $\sigma_m^2$ is the variance of Miras in that galaxy; the prior variance $\sigma_b^2$ is the variance of the light curve amplitude. These prior parameters can be determined using previous studies. For example, in §\[sec.evaluation\], we use well-sampled light curves of LMC Miras [@Soszynski2009] to obtain values of these parameters. It is advisable to check the sensitivity of these prior specifications.

The benefit of using priors on $m$ and ${\boldsymbol{\beta}}$ is three-fold: first, they introduce regularization by using information from early studies; second, they provide a natural device for separating the estimation of frequency and the light curve signal component using Bayesian integration when the parameter of interest is the frequency; lastly, the regularization parameter ${\boldsymbol{\theta}}$ of the non-parametric function is allowed to be chosen by the maximum likelihood, without resorting to the computationally expensive cross-validation method.

In summary, we have built the following hierarchical model for a Mira light curve:

$$\label{eqn:hier}
\begin{split}
& y_i | m,{\boldsymbol{\beta}}, g(t_i)  \sim {\mathcal{N}}(g(t_i), \sigma_i^2), \\
& g(t)  = m + {\mathbf{b}}_f(t)^T{\boldsymbol{\beta}}+ h(t), \\
& m   \sim {\mathcal{N}}(m_0, \sigma_m^2), {\boldsymbol{\beta}}\sim {\mathcal{N}}({\mathbf{0}},\sigma_b^2{\mathbf{I}}),\\
& h(t)| {\boldsymbol{\theta}}\sim  {\mathcal{GP}}(0,k_{{\boldsymbol{\theta}}} (t,t')),
\end{split}$$

where ${\boldsymbol{\theta}}$ and $f$ are fixed parameters. In this model, the frequency parameter $f$ is of key interest to our study. We do not perform a fully Bayesian inference by imposing a prior distribution on $f$ because the likelihood function of $f$ is highly irregular, with numerous local maxima, and Monte Carlo computation of the posterior is expensive and intractable for large astronomical surveys.

Previously, @Baluev2013 applied a Gaussian process model to study the impact of red noise in radial velocity planet searches. While his maximum likelihood method is a classical frequentist approach in statistics, our approach can be considered as a hybrid of Bayesian and frequentist approaches. We treat the parameter of interest $f$, and the parameters for the kernel ${\boldsymbol{\theta}}$ of the Gaussian process as fixed, and impose a prior distribution on other parameters. This is similar to the type-II maximum likelihood estimation of parameters of a Gaussian process or regularization parameters in function estimation; see §5.2 of @Rasmussen2005. From the Bayesian point of view, ${\boldsymbol{\theta}}$ and $f$ are treated as hyper-parameters that in turn are estimated by the empirical Bayes method. Because the Gaussian process plays a critical role in modeling departure of light curves from periodicity, we may also refer to our model more precisely as the nonlinear SP Gaussian process model.

Estimation of the frequency and the periodogram
-----------------------------------------------

Let ${\mathbf{y}}= (y_1,y_2,\cdots, y_n)$ be the observation vector of the magnitudes of a light curve. By integrating out $m,{\boldsymbol{\beta}}$ and ${\mathbf{h}}$ from the joint distribution given by Eqn. \[eqn:hier\], we get the marginal distribution of ${\mathbf{y}}$, $p({\mathbf{y}}|{\boldsymbol{\theta}}, f)$, which is a multivariate normal with mean ${\boldsymbol{\mu}}= m_0{\mathbf{1}}$ and covariance matrix $${\mathbf{K}}_y = \left( \sigma_m^2 + \sigma_b^2 {\mathbf{b}}_f(t_i)^T {\mathbf{b}}_f(t_j) +
 k_{{\boldsymbol{\theta}}}(t_i, t_j) +\sigma^2_i\delta_{ij} \right)_{n\times n},$$ where $\delta_{ij} = 1$ if $i=j$ and $\delta_{ij} = 0$ if $i\neq j$. Therefore, the log likelihood of ${\boldsymbol{\theta}}$ and $f$ is $$\begin{split}
Q({\boldsymbol{\theta}},f) = &\log(p({\mathbf{y}}| {\boldsymbol{\theta}}, f))\\
             = &-\frac{1}{2} ({\mathbf{y}}-m_0{\mathbf{1}})^T{\mathbf{K}}_y^{-1}({\mathbf{y}}-m_0{\mathbf{1}}) \\
               &\qquad -\frac{1}{2}\log\det {\mathbf{K}}_y -\frac{n}{2}\log(2\pi) \label{eqn:mainObj}.
\end{split}$$

The maximum likelihood estimator of ${\boldsymbol{\theta}}$ and $f$ is obtained by maximizing $Q({\boldsymbol{\theta}},f)$. Since the likelihood function is differentiable with respect to ${\boldsymbol{\theta}}$ but highly multimodal in the parameter $f$, standard optimization methods cannot be directly used to jointly maximize over ${\boldsymbol{\theta}}$ and $f$.

We adopt a profile likelihood method as follows. For each frequency $f$ over a dense grid, we compute the maximum likelihood estimator $\widehat{{\boldsymbol{\theta}}}_f = \operatorname*{arg\,max}_{{\boldsymbol{\theta}}} Q({\boldsymbol{\theta}},f)$. This can be done using the quasi-Newton method. Then we perform a grid search to find the maximum profile likelihood estimator of $f$, i.e., $$\label{eqn:fhat}
\hat{f} = \operatorname*{arg\,max}_f Q(\widehat{{\boldsymbol{\theta}}}_f, f)\, ,$$ the estimated period is $\hat{P}=1/\hat{f}$. The details of the algorithm are given in §\[sec.quasi\]. The profile log-likelihood as a function of the frequency $f$ is adopted as the *periodogram* for our model, $$S_{SP}(f) = Q(\widehat{{\boldsymbol{\theta}}}_f, f)\, . \label{eqn:periodogram}$$ It contains the spectral information of the signal. The frequency of the dominant harmonic component is expected to be the location of the peak of this profile likelihood.

Computation of the periodogram {#sec.quasi}
------------------------------

Now we present the details of computing the profile likelihood. Because $Q({\boldsymbol{\theta}},f)$ is highly multimodal in the frequency parameter $f$, we follow the commonly used strategy of optimization through grid search. On the other hand, since $Q({\boldsymbol{\theta}},f)$ is differentiable in parameter ${\boldsymbol{\theta}}$, the quasi-Newton method can be employed to optimize over ${\boldsymbol{\theta}}$ for fixed $f$, and obtain the profile likelihood (Eqn. \[eqn:periodogram\]). The gradient of the log likelihood (Eqn. \[eqn:mainObj\]) with respect to $\theta_j (j=1,2)$ is $$\frac{\partial}{\partial \theta_j} Q({\boldsymbol{\theta}},f)\!=\! \frac{1}{2}\mathrm{tr}\left(({\boldsymbol{\alpha}}{\boldsymbol{\alpha}}^T\!-\!{\mathbf{K}}_y^{-1})\frac{\partial {\mathbf{K}}_y}{\partial \theta_j}\right)$$ where ${\boldsymbol{\alpha}}\!=\!{\mathbf{K}}^{-1}_y({\mathbf{y}}-m_0{\mathbf{1}})$. In general, the objective function for the Gaussian process model is not convex in its kernel parameters ${\boldsymbol{\theta}}$ and global optimization cannot be guaranteed. Fig. \[fig:lc.surface\] shows a surface plot of $Q({\boldsymbol{\theta}},f)$ as a function of ${\boldsymbol{\theta}}$ for one simulated light curve, with $f$ fixed at the true frequency. The surface exhibits unimodality in this case, although it is not convex.

The computation involved in calculating the profile likelihood through the quasi-Newton method can be intensive. Since the objective function (Eqn. \[eqn:mainObj\]) is non-convex in ${\boldsymbol{\theta}}$, generally multiple starting points should be attempted to find the global optimizer when applying the quasi-Newton method. In addition, evaluating the objective function and the gradient function requires inversion of the covariance matrix whose computation cost is of the order $O(n^3)$. During each quasi-Newton iteration, these evaluations could be repeated several times because multiple step size might be attempted. To make the computation more challenging, all of the above needs to be repeated at hundreds or even thousands of densely gridded $f$s per light curve. Furthermore, the method may need to be applied to hundreds of thousands or millions of light curves from large astronomical surveys.

In order to speed up computation over the dense grid of frequency values, we use the result of applying the quasi-Newton method at one frequency value as a warm start for the subsequent frequency value. Specifically, the optimizer $\widehat{{\boldsymbol{\theta}}}_f$ and its approximate inverse Hessian matrix are provided as quantities to start the quasi-Newton iterations for the next frequency value on the dense grid. When the initial point is near the local minimizer and the inverse Hessian matrix is a good approximation to the true Hessian matrix, the quasi-Newton algorithm will converge at superlinear rate; the step size of $\alpha=1$ will be accepted by the Wolfe descent condition, avoiding evaluation of the objective function multiple times to determine the appropriate step size during each iteration [see Ch. 6 of @Nocedal2006 for a more rigorous mathematical discussion]. We find that a warm start can speed up the computation significantly but sometimes we need to restart with random initial values to ensure convergence to the global optimum. The pseudocode provided in the Appendix describes our algorithm.

Estimation of the signal and its components
-------------------------------------------

After the parameters $f$ and ${\boldsymbol{\theta}}$ are fixed at their maximum likelihood estimates $\widehat{f}$ and $\widehat{{\boldsymbol{\theta}}}_{\hat{f}}$, we can perform the inference of the light curve signal $g(t)$ and its components in the standard Bayesian framework. Interested readers may consult Ch. 2 of @Rasmussen2005 for a detailed discussion of this topic.

Firstly, we could obtain the posterior distribution of ${\boldsymbol{\gamma}}= (m,{\boldsymbol{\beta}}^T)$, the parameters for the long run average magnitude and the exactly periodic term. The prior of ${\boldsymbol{\gamma}}$ is ${\mathcal{N}}({\boldsymbol{\gamma}}_0,{\boldsymbol{\Sigma}}_\gamma)$ with ${\boldsymbol{\gamma}}_0=(m_0,0,0)^T$ and ${\boldsymbol{\Sigma}}_\gamma = \mathrm{diag}(\sigma_m^2,\sigma_b^2,\sigma_b^2)$. Its posterior distribution is ${\boldsymbol{\gamma}}| {\mathbf{y}}\sim  {\mathcal{N}}(\bar{{\boldsymbol{\gamma}}}, \bar{{\boldsymbol{\Sigma}}}_\gamma)$ with $$\label{equ:gammapost}
\begin{split}
\bar{{\boldsymbol{\gamma}}}  = & \left({\mathbf{H}}^T{\mathbf{K}}_c^{-1}{\mathbf{H}}+{\boldsymbol{\Sigma}}_\gamma^{-1}\right)^{-1} \\
                 & \left({\boldsymbol{\Sigma}}_\gamma^{-1}{\boldsymbol{\gamma}}_0+{\mathbf{H}}^T{\mathbf{K}}_c^{-1}{\mathbf{y}}\right)\, , \\
\bar{{\boldsymbol{\Sigma}}}_\gamma  = &  \left({\mathbf{H}}^T{\mathbf{K}}_c^{-1}{\mathbf{H}}+ {\boldsymbol{\Sigma}}_\gamma^{-1}\right)^{-1}, 
\end{split}$$ [where]{} $${\mathbf{h}}(t)\!=\!(1,{\mathbf{b}}_{\hat{f}}(t)^T)^T, {\mathbf{H}}\!=\!({\mathbf{h}}(t_1),{\mathbf{h}}(t_2),\cdots,{\mathbf{h}}(t_n))^T,$$ and$\,{\mathbf{K}}_c\!=\!\big(k_{\widehat{{\boldsymbol{\theta}}}_{\hat{f}}}(t_i,\!t_j)\!+\!\sigma_i^2\delta_{ij}\!\big)\!_{{\tiny\textit{n}}\times\!{\tiny\textit{n}}}\,$with$\,\widehat{f}\,$and$\,\widehat{{\boldsymbol{\theta}}}_{\hat{f}}\,$plugged in.

Consider the prediction of light curve magnitude at a specific time point $t^*$. Define the vector ${\mathbf{k}}^* = (k_{\widehat{{\boldsymbol{\theta}}}}(t^*,t_1),\ \cdots,\ k_{\widehat{{\boldsymbol{\theta}}}}(t^*,t_n))^T$. Conditional on $({\mathbf{y}}, {\boldsymbol{\gamma}})$, the distribution of $g(t^*)| {\mathbf{y}}, {\boldsymbol{\gamma}}$ is a multivariate normal with mean $ {\mathbf{h}}(t^*)^T \gamma + {\mathbf{k}}_{{\boldsymbol{\theta}}}(t^*,{\mathbf{t}}) {\mathbf{K}}_c^{-1} ({\mathbf{y}}- {\mathbf{H}}\gamma)$ and variance $k_{\widehat{{\boldsymbol{\theta}}}}(t^*,t^*) - ({\mathbf{k}}^*)^T {\mathbf{K}}_c^{-1}  {\mathbf{k}}^*$. With the posterior distribution of ${\boldsymbol{\gamma}}$ given in Eqn. \[equ:gammapost\], we are able the remove ${\boldsymbol{\gamma}}$ from the above conditional distribution of $g(t^*)$. Finally, we get the posterior distribution of the signal at $t^*$ as $g(t^*)| {\mathbf{y}}\sim {\mathcal{N}}(\bar{g}^*, \bar{\sigma}_{g^*}^2)$ with $$\begin{split}
\bar{g}^* = & {\mathbf{h}}(t^*)^T \bar{\gamma} + {\mathbf{k}}(t^*,{\mathbf{t}}) {\mathbf{K}}_c^{-1}
            ({\mathbf{y}}- {\mathbf{H}}\bar{\gamma})\, ,\\
\bar{\sigma}^2_{g^*} = & k_{\widehat{{\boldsymbol{\theta}}}}(t^*,t^*) - ({\mathbf{k}}^*)^T {\mathbf{K}}_c^{-1}
                         {\mathbf{k}}^*+
{\mathbf{r}}^T\bar{{\boldsymbol{\Sigma}}}_\gamma{\mathbf{r}}\, , 
\end{split}
\label{equ:gpredict}$$ where ${\mathbf{r}}={\mathbf{h}}(t^*)-{\mathbf{H}}^T{\mathbf{K}}_c^{-1}{\mathbf{k}}^*$.

Simulation of M33 light curves {#sec.construct.test}
==============================

It is not possible to evaluate the period estimation accuracy of our method directly on the M33 data because the “ground truth” is unknown. Instead, we construct a test data set by smoothing the well-sampled OGLE light curves to infer continuous functions, then resample these functions to match the observational patterns of the M33 data, and at last add noise to the light curves. This data set can serve as a testbed for future studies of comparing different period estimation methods. We will now describe the M33 observations and the construction of the test data set. As the whole simulation procedure is a complicated process, we will discuss its components in detail from §4.1 to §4.4. The whole simulation procedure will be summarized in §4.5.

Characteristics of the M33 observations
---------------------------------------

Most of the disk of M33 was observed by the DIRECT [@Macri2001] and M33SSS [@Pellerin2011] projects in the $BVI$ bands, with a combined baseline of $7-9$ years and a sampling pattern that depends on the exact location within the disk (see Fig. \[fig:obs.gaps\]). The large area of coverage and long baseline of these observations make them suitable for Mira searches. We use the $I$-band observations to carry out the simulations, as this is the wavelength range where Miras are brightest (out of the three bands used by these projects). Detailed descriptions of the M33 observations can be found in the above referenced papers. We use the data products from a new reduction that will be presented in a companion paper (W. Yuan et al. 2016, in prep.). $I$-band light curves are available for $\sim 2.5\times 10^5$ stars, with a median of 44 measurements and a maximum of 170.

.

We model the relation between a magnitude measurement $m$ and its uncertainty $\sigma$ as $$\sigma = a(t_i',F)^{[m-b(t_i',F)]} + c(t_i',F)\, , \label{equ.sigma.mag}$$ for each observation field $F$ and each observation night $t_i'$, where $a(t_i',F)$, $b(t_i',F)$ and $c(t_i',F)$ are field- and night-specific constants. There are 31 different fields in total, $F=0,1,\cdots, 9, a,b,\cdots, u$. The parameters are determined via least-squares fitting using all the measurements for the specific field $F$ and night $t_i'$. Fig. \[fig:simu.sigmag\] shows the $m-\sigma$ relation for a typical field.

![$m-\sigma$ relation for a given night and field within M33. The solid red line is the best-fit relation using the empirical function $\sigma = a^{(m-b)} + c$, with $a=2.666,\,b=23.117,\,c=0.008$.[]{data-label="fig:simu.sigmag"}](fig05.eps){width="49.00000%"}

In order to test the SP periodogram we need sparsely sampled, moderately noisy Mira light curves with known periods. Thus, we characterize the sampling patterns and noise levels of the M33 observations and simulated Mira light curves of known periods using the OGLE observations of these objects in the LMC.

Matching the M33 observation pattern
------------------------------------

The first step in simulating a Mira light curve is to randomly select a sampling pattern based on the light curve of an actual star in some field $F$, $\{t_i'\}_{i=1}^n$ with $n\in [10,170]$. A random time shift $s$ is added, $t_i = t_i' + s$ for $i=1,2,\cdots, n$. The random shift $s$ follows a uniform distribution over the interval $[0, P_0]$, where $P_0$ is the true period of the LMC Mira selected during the artificial light curve generation process. This helps to simulate a large number of unique light curves sampled at random phases using the limited number of template light curves.

The Mira template light curves
------------------------------

The template Mira light curves are obtained by using our SP model to fit the Mira light curves in the LMC, collected by the OGLE project [@Soszynski2009]. A total number of 1663 Miras have been observed in $I$ with very high accuracy, excellent phase coverage, and a long baseline (the median and mean number of observations are 466 and 602, respectively, with a baseline of $\sim 7.5$ years for most fields). Because the LMC light curves are densely sampled with high quality, we can adopt a more complicated model to provide a higher fidelity fit. Following §5.4.3 of @Rasmussen2005, instead of Eqn. \[eqn:basicDecomposition\], the signal light curve $g(t)$ is decomposed into $$\label{eqn:fullgpmodel}
g(t)  = m + l(t) + q(t) + h(t),$$ where $m$ is the long run average magnitude, $l(t)$ is the long-term (low-frequency) trend across different cycles, $q(t)$ is the periodic term, and $h(t)$ is small-scale (high-frequency) variability within each cycle. The latter three terms are modeled by the Gaussian process with different kernels. In particular, we use the squared exponential kernel $k_l(t_1,t_2) = \theta_1^2 \exp(-\frac{1}{2}\frac{(t_1-t_2)^2}{\theta_2^2})$ for $l(t)$, another squared exponential kernel $k_h(t_1,t_2) = \theta_6^2 \exp(-\frac{1}{2}\frac{(t_1-t_2)^2}{\theta_7^2})$ for $h(t)$, and lastly a periodic kernel $$\begin{split}
k_q(t_1,t_2) = \theta_3^2  \exp \bigg ( -\frac{1}{2}&\frac{(t_1-t_2)^2}{\theta_4^2}  \\
 &- \frac{2\sin^2(2\pi f(t_1-t_2))}{\theta_5^2} \bigg )
\end{split}$$ for $q(t)$. Note the periodic kernel allows the light curve amplitude to change across cycles. The maximum likelihood method is applied to fit each LMC light curve, fixing $f$ to the OGLE value and solving for the unknown parameters $(\theta_1,\theta_2, \cdots, \theta_7)$. Fig. \[fig:complexDecomposition\] is an illustration of the model fitting result using Eqn \[eqn:fullgpmodel\] based on the same light curve as in Fig. \[fig:onemira\]. Notice that the more complex model in Fig. \[fig:complexDecomposition\] is only suitable for a densely sampled light curve.

Once the sampling pattern is chosen, one of the template light curves will be selected according to the luminosity function described in the next subsection. With the selected template, the magnitude of the simulated light curve signal at $t_i'$ with shift $s$ is $g(t_i'+s)$, which is computed with Eqn. \[eqn:fullgpmodel\] in a similar way as Eqn. \[equ:gpredict\].

![Light curve of a Mira in the LMC observed by OGLE (black points), decomposed following Eqn \[eqn:fullgpmodel\]. Top panel: the fitted light curve; second panel: long-term signal, $m+l(t)$; third panel: periodic term, $m + q(t)$; bottom panel: stochastic variations, $m+h(t)$.[]{data-label="fig:complexDecomposition"}](fig06.eps){width="49.00000%"}

Matching the luminosity function\
to the M33 observations
---------------------------------

While the OGLE observations of LMC Miras are deep enough to detect these objects over their entire range of luminosities, the M33 observations become progressively more incomplete for fainter and redder objects. We derived an empirical completeness function for the M33 observations as follows. We fitted the observed luminosity function $\mathcal{F}_0(I)$ using an exponential for $I \in [18.5,20]$ mag and extrapolated to fainter magnitudes, obtaining $\mathcal{F}_1(I)$. The empirical completeness function is then $\mathcal{C}(I) = \mathcal{F}_1(I) / \mathcal{F}_0(I)$.

We randomly picked $\{t_i'\}_{i=1}^n$ from the M33 light curves. For each $\{t_i'\}_{i=1}^n$, we selected a (LMC-based) template using $\mathcal{C}(I+6.2)$ as the probability distribution. The value of $+6.2$ mag accounts for the approximate difference in distance modulus between the LMC and M33. In this way the resulting luminosity function of the simulated light curves is statistically the same as that of the real M33 observations.

The simulation procedure
------------------------

With all the components discussed above, we are able to present the whole simulation procedure here. In order to generate one simulated Mira light curve matching the sampling characteristics of the M33 observations, the first step is to randomly select a sampling pattern $\{t_i'\}_{i=1}^n$, and then add a random shift $s$, $t_i = t_i' + s$, $i=1,2,\cdots, n$. The second step is to randomly select a template light curve according to the luminosity function, then compute the light curve signal $g(t_i'+s)$ for the selected sampling pattern $\{t_i'\}_{i=1}^n$. The third step is to use the best-fit relations (Eqn. \[equ.sigma.mag\]) to add photometric noise via $$y_i = g(t_i' + s)+6.2+ \sigma_i\epsilon_i\, ,$$ where $+6.2$ mag is the approximate relative distance modulus, $\epsilon_i$ is drawn from $\mathcal{N}(0,1)$, and $\sigma_i$ is computed from $$\sigma_i = a(t_i',F)^{[g(t_i)+6.2-b(t_i',F)]} + c(t_i',F).$$ for the selected observation pattern $t_i'$ and field $F$. Following this procedure, we generate one simulated light curve $\{t_i',y_i,\sigma_i\}_{i=1}^n$. The procedure is repeated until $10^5$ suitable light curves are generated, excluding any with $<10$ data points or sampling on $<7$ nights.

Performance evaluation {#sec.evaluation}
======================

Having generated the test data set, we evaluate the performance of the SP model and compare it with the GLS model. We choose prior parameters for the SP model of $m_0 = 15.62 + 6.2$, $\sigma_m = 10$ and $\sigma_b =1$. The adopted value of $m_0$ is the average $I$ magnitude of Miras in the LMC and once again $+6.2$ is the approximate relative distance modulus between M33 and the LMC. The values of $\sigma_m$ and $\sigma_b$ are larger than those derived from the LMC samples in order to make those priors non-informative. Although fitting the SP model is computationally slower than the LS model, we find that our model gives an overall improvement in various metrics. For both methods, the periodograms are computed on a dense frequency grid from $1/2000$ to $1/100$ with a spacing of the order of $10^{-5}$. For the GLS method, we chose a spacing of (0.05/time span) or $\sim2.5\times 10^{-5}$, which results in optimal performance for this simulation. For our SP method, we chose a slightly smaller value of $10^{-5}$ to facilitate the warm start mechanism in our algorithm (see Appendix) given that small changes in frequency result in tiny changes of the objective function.

The aliasing effect
-------------------

We fit the entire simulated data set using the SP model. Fig. \[fig:lc.and.spec\] gives an example of a simulated light curve and its SP periodogram (Eqn. \[eqn:periodogram\]). In this example, the true frequency (labeled by the blue dotted line) is successfully recovered.

Aliasing frequencies at $f\pm 1/365$ d affect most periodograms when dealing with sparsely observed astronomical data. The red dashed line in Fig. \[fig:lc.and.spec\] indicates the aliasing frequency at $f+1/365$ where a strong peak exists. This is not a rare case, and for some light curves the one-year beat aliasing frequencies have higher log likelihoods than the true frequencies. Fig. \[fig:fvsf\] compares the recovered and true frequencies for all simulated light curves. Two secondary strips parallel to the main one and offset by $\pm0.00274$ represent $\hat{f} = f \pm 1/365$, respectively. Other aliasing frequencies, such as $2f$, $3f$, $0.5f$, etc., are also noticeable. Lastly, due to the sampling pattern of some light curves, the side lobes of the main peak can be higher than the central value. These manifest as close parallel strips to the aforementioned features.

Accuracy assessment
-------------------

The estimated frequency is considered as correct if $\Delta f=|\hat{f}-f_0|<C_f$ for each light curve. The estimation accuracies for the two methods are summarized in Table \[tbl:accuracy\] for several different values of $C_f$. We choose $C_f = 2.7\times 10^{-4}$ to stringently bind the one-to-one strip in Fig. \[fig:fvsf\]. Overall, the SP correctly estimates the period for 69.4% of the light curves, while the LS model has an accuracy of 63.6%. The improvement of SP over LS is more evident for C-rich Miras, with about 10% higher accuracy, while the improvement for O-rich Miras is smaller, with about 3% higher accuracy. The difference in performance arises because C-rich Miras often exhibit larger stochastic deviations that can be better captured by the SP model, while O-rich Miras have more stable light curves that can be modeled reasonably well with the LS method. We also compute the estimation accuracy of each method by grouping the light curves according to the number of observations, as shown in the left panels of Fig. \[fig:accuracy\]. The top and bottom rows show results for C- and O-rich Miras, respectively. The performance difference is once again more evident in the C-rich category.

Note that accuracy is not a monotonic function of the number of observations, implying this is not a good indicator [*per se*]{} of the information content of the light curves for frequency (period) estimation. Thus, we define another metric, called [*phase coverage*]{}. Recall that the times of observation for a given light curve are $t_1,t_2,\cdots,t_n$. Given a period of $P$, these are converted into corresponding phases by $s_i = (t_i\ \mathrm{ mod }\ P)/P,\ i=1,2,\cdots, n\,$ in the closed interval $[0,1]$. Now, define $$J = \Big(\bigcup_i (s_i-l,s_i+l)\Big)\cap [0,1]\, ,$$ for a specific $l>0$, the phase coverage can be measured by $\lambda(J)$ where $\lambda(\cdot)$ is the Lebesgue measure (we choose $l=0.02$). $\lambda(J)$ describes the “length” of the union of the intervals $J$. For example, $\lambda(J) = 0.1$ for $J = (0.1,0.2)$, and $\lambda(J) = 0.2$ for $J = (0.1,0.2)\cup (0.5, 0.6)$.

We divide the light curves into 100 groups such that their $\lambda(J)$ is in one of the intervals $(k/100,(k+1)/100]$ for $k=0,1,\cdots,99$ and compute the estimation accuracy for each subset. The results for the two models are plotted in the middle column of Fig. \[fig:accuracy\]. Now the estimation accuracy is monotonically increasing as a function of phase coverage. The accuracy improvement of our method is highest when the phase coverage is around 0.5 for C-rich Miras. As the phase coverage approaches the extremes (0 or 1), the performance difference between the two methods diminishes. At $\lambda(J) \approx 0$, both methods will fail because this is a hopeless situation. At the other extreme, when $\lambda(J)\approx 1$ and abundant information is available for frequency estimation, both methods have an accuracy close to 1.

The periodogram $S_{\rm SP}(f)$ of our model defined in Eqn. \[eqn:periodogram\] provides more information than just the optimal frequency. Suppose $f_1$ is the largest local maximal (global maximum) of $S_{\rm SP}(f)$, and $f_2$ is the second largest local maximal of $S_{\rm SP}(f)$. Now define [conf]{} $=S_{\rm SP}(f_1)-S_{\rm SP}(f_2)\ge 0$. The value of [conf]{} serves as a confidence measurement of the global optimal estimate in Eqn. \[eqn:fhat\]. Larger values of [conf]{} indicate smaller uncertainty in our estimate, and thereby the estimate is more reliable. Now, let $c_0$ be the smallest value, and let $c_1, \dots, c_{100}$ be the $1$st–$100$th percentiles of all the [conf]{} values computed for all the light curves. Each light curve can be assigned to a percentile group if its [conf]{} is in $(c_{k-1},c_{k}]$ for some $k\in\{1,2,\cdots,100\}$. After assigning all light curves by [conf]{} to their corresponding percentile groups, the estimation accuracy in each group can be computed. The same procedure is applied to the GLS model, with the $p$-value of the F-statistics given in Eqn. \[eq:glsP\] for the top peak being used as its [conf]{} measurement. The result is plotted in the right column of Fig. \[fig:accuracy\]. The accuracy of our SP method is much higher than the LS model in the top 40 groups. In particular, the accuracy of our method is higher than 90% in the top 20 groups for both C- and O-rich Miras.

Light curves with high values of [conf]{} are particularly reliable for constructing Period-Luminosity relations (hereafter, PLRs) based on the “Wesenheit” function [@Madore1982]. This function enables a simultaneous correction for the effects of dust attenuation and finite width of the instability strip by defining a new magnitude $W_I\!=\!I\!-\!1.55(V\!-\!I)$, where $V$ and $I$ are the mean magnitudes in those filters. Figure \[fig:pl.compare\] compares PLRs based on $W_I$ magnitudes and periods determined by OGLE and estimated with each of the two models. The top and bottom rows display the PLRs for C- and

[ccrrr]{}

& SP & 58.1 & 55.3 & 56.5\
& LS & 49.4 & 51.6 & 50.6\
& SP & 69.6 & 63.8 & 66.3\
& LS & 60.1 & 60.6 & 60.4\
& SP & 73.5 & 66.2 & 69.4\
& LS & 63.7 & 63.5 & 63.6\
& 43,116 & 56,884 &

[O-rich Miras, respectively. The leftmost column shows the PLRs based on the actual OGLE periods, while the next two sets of columns show the corresponding relations based on SP or LS periods for the simulated light curves with the top 10% and 40% values of [conf]{}.]{}

In order to provide a quantitative comparison of the improvement obtained with our SP method, we calculated the dispersion of the actual $W_I$ PLRs and their recovered counterparts as a function of [conf]{} value as follows, separately for C- and O-rich Miras. First, we selected all objects of a given class with $2\!<\!\log P\!<\!3$. If the @Soszynski2009 catalog did not provide a $V$ measurement for a given variable, the missing value was estimated through linear interpolation of the $(I,V\!-\!I)$ relation for objects of the same class within $|\Delta\log P|<0.05$ dex. We fitted a quadratic PLR $$m=a+b(\log P-2.3)+c(\log P-2.3)^2$$ with iterative $3\sigma$ clipping (removing $\sim 5$% of the data). We then computed the dispersion of the initially selected OGLE sample about the best-fit relation, including outliers. This yielded “benchmark” dispersions of 0.45 & 0.54 mag for C- & O-rich variables, respectively. Keeping the best-fit relation fixed, we computed the dispersion of recovered PLRs using all artificial light curves within a certain range of [conf]{} (top 10%, top 20%, $\dots$), using the periods and [conf]{} values derived by the SP or the LS method. As in the case of the OGLE samples, we only considered objects with $2\!<\!\log P\!<\!3$. The results are plotted in Fig. \[fig:pl.sig\]. The SP subsamples exhibit lower

[llrrrr]{}

00082 & O & 14.241 & 16.509 & 164.84 &\
00094 & C & 15.120 & 18.885 & 332.30 &\
00098 & C & 15.159 & 17.921 & 323.10 &\
00115 & C & 14.932 & 16.947 & 176.13 &\
00355 & O & 14.199 & 16.219 & &\

[(or at worst, equal) dispersions than their LS counterparts for all percentiles and for both subtypes. As discussed previously, the improvement provided by our method is strongest for C-rich Miras and diminishes in significance as one includes light curves with progressively lower confidence values.]{}

Summary {#sec.discussion}
=======

In this paper, we developed a nonlinear SP Gaussian process model for estimating the periods of sparsely sampled quasi-periodic light curves, motivated by the desire to detect Miras in an existing set of observations of M33. We conducted a large-scale high-fidelity simulation of Mira light curves as observed by the DIRECT/M33SSS surveys to compare our model with the GLS method. Our model shows improved accuracy under various metrics. The simulation data set is provided as a testbed for future comparison with other methods. The SP model will be used in a companion paper to search for Miras in M33, estimate their periods, and study the resulting PLRs.

 

SH was partially supported by Texas A&M University-NSFC Joint Research Program. WY & LMM acknowledge financial support from the NSF through AST grant \#1211603 and from the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University. JZH was partially supported by NSF grant DMS-1208952. The authors acknowledge the Texas A&M University Brazos HPC cluster that contributed to the research reported here.

[[*Simulated light curves:*]{} A tarfile, containing $10^5$ simulated light curves are generated following the procedure of §\[sec.construct.test\]. Each light curve is stored in one file with three columns: MJD, $I$ magnitude, and uncertainty. The file name, e.g., lc006788.dat is generated sequentially and is only meant for bookkeeping purposes. A mapping between simulated light curve ID and the original OGLE object is given in the file “lc.dat”, which can also be found in the tarfile.]{}

[[*Mira variables:*]{} Table\[tbl:oglemira\] summarizes the relevant properties of OGLE LMC Miras from @Soszynski2009 that were used to simulate the light curves: OGLE ID, main period and mean $I$ & $V$ magnitudes. It includes some extrapolated values of $V$ for objects with missing data (suitably identified with a “\*”). This table can be used to compare true versus derived periods and to generate Period-Luminosity relations.]{}

[[*Software:*]{} The related software package, `varStar`, has been released under a GPL3 license [@He2016]. The active software development repository can be found at [github.com/shiyuanhe/varStar](github.com/shiyuanhe/varStar).]{}

* *Quasi-Newton’s Method with Grid Search**\
**Input:** Maximal and minimal trial frequencies $f_M\!>\!f_m\!>\!0$; frequency step $\Delta f$; $n$ observations $\{t_i,y_i,\sigma_i\}$.\
**Output:** Periodogram $S(f)$ evaluated at the trial frequencies.\

Initialize ${\boldsymbol{\theta}}^{(0)}$ and $\mathbf{H}^{(0)}$, and $f\gets f_m$;

$p \gets 0$;

$\mathbf{t}_p \gets -\mathbf{H}^{(p)} \frac{\partial}{\partial {\boldsymbol{\theta}}}
Q({\boldsymbol{\theta}}^{(p)}, f)$;

${\boldsymbol{\theta}}^{(p+1)}\gets {\boldsymbol{\theta}}^{(p)} +\alpha_p \mathbf{t}_p$ and the step size $\alpha_p$ satisfying the Wolfe condition;

$\mathbf{d}_p \gets {\boldsymbol{\theta}}^{(p+1)}-{\boldsymbol{\theta}}^{(p)}$, $\mathbf{e}_p \gets \frac{\partial}{\partial {\boldsymbol{\theta}}} Q({\boldsymbol{\theta}}^{(p+1)}, f)- 
\frac{\partial}{\partial {\boldsymbol{\theta}}} Q({\boldsymbol{\theta}}^{(p)}, f)$;

$\rho_p= 1/\mathbf{d}_p^T\mathbf{e}_p$;

$\mathbf{H}^{(p+1)}\gets (\mathbf{I}-\rho_p
\mathbf{d}_p\mathbf{e}_p^T)
\mathbf{H}^{(p)} (\mathbf{I}-\rho_p \mathbf{e}_p\mathbf{d}_p^T) 
+\rho_p\mathbf{d}_p\mathbf{d}_p^T$;

$p \gets p + 1$;

$\widehat{{\boldsymbol{\theta}}}_{f} \gets {\boldsymbol{\theta}}^{(p)}$, and $S(f)\gets Q(\widehat{{\boldsymbol{\theta}}}_{f},f)$;

$\mathbf{H}^{(0)} \gets \mathbf{H}^{(p)}$ and ${\boldsymbol{\theta}}^{(0)}\gets {\boldsymbol{\theta}}^{(p)}$;
