# Simultaneous inference of periods and period-luminosity relations for Mira
variable stars

Shiyuan<EMAIL_ADDRESS>[ Zhenfeng
<EMAIL_ADDRESS>[ Wenlong<EMAIL_ADDRESS>[ Lucas M<EMAIL_ADDRESS>[ Jianhua Z.
<EMAIL_ADDRESS>[  Institute of Statistics and Big Data,
Renmin University of China<EMAIL_ADDRESS>Microsoft,
<EMAIL_ADDRESS>Department of Physics and Astronomy, Johns Hopkins
University<EMAIL_ADDRESS>Department of Physics and Astronomy, Texas A&M
University<EMAIL_ADDRESS>Department of Statistics, Texas A&M University
<EMAIL_ADDRESS>
###### Abstract

The Period–Luminosity relation (PLR) of Mira variable stars is an important
tool to determine astronomical distances. The common approach of estimating
the PLR is a two-step procedure that first estimates the Mira periods and then
runs a linear regression of magnitude on log period. When the light curves are
sparse and noisy, the accuracy of period estimation decreases and can suffer
from aliasing effects. Some methods improve accuracy by incorporating complex
model structures at the expense of significant computational costs. Another
drawback of existing methods is that they only provide point estimation
without proper estimation of uncertainty. To overcome these challenges, we
develop a hierarchical Bayesian model that simultaneously models the quasi-
periodic variations for a collection of Mira light curves while estimating
their common PLR. By borrowing strengths through the PLR, our method
automatically reduces the aliasing effect, improves the accuracy of period
estimation, and is capable of characterizing the estimation uncertainty. We
develop a scalable stochastic variational inference algorithm for computation
that can effectively deal with the multimodal posterior of period. The
effectiveness of the proposed method is demonstrated through simulations, and
an application to observations of Miras in the Local Group galaxy M33. Without
using ad-hoc period correction tricks, our method achieves a distance estimate
of M33 that is consistent with published work. Our method also shows superior
robustness to downsampling of the light curves.

variational inference,

hierarchical Bayesian modeling,

astrostatistics,

Mira variables,

###### keywords:

, , ,
and

##  1 Introduction

One essential step to determine the age and composition of our Universe is to
calibrate its current expansion rate (commonly referred to as the “Hubble
constant” or $H_{0}$ in astronomy) using the cosmological distance ladder.
Observational astronomers have extensively relied on the Period–Luminosity
relations (PLRs) of classical Cepheid variable stars to measure distances to
galaxies within 150 million light years of the Milky Way, enabling the
calibration of farther-reaching techniques (such as type Ia supernovae) to
obtain increasingly more robust estimates of $H_{0}$ (Freedman et al., 2001;
Riess et al., 2011; Riess et al., 2016; Riess et al., 2019). Recent studies
(Whitelock et al., 2014; Whitelock and Feast, 2014; Yuan et al., 2017; Huang
et al., 2018) have shown that Mira variable stars (hereafter, Miras) also
exhibit tight PLRs at near-infrared (NIR) wavelengths and thus can serve as
independent distance indicators in lieu of Cepheids. Miras are highly-evolved
stars with two subclasses based on their atmospheric chemistry: oxygen- and
carbon-rich (O- and C-rich, respectively). They constitute a class of long-
period variable stars whose luminosity changes dramatically (in a quasi-
periodic manner) as a function of time. One pulsation cycle of Miras usually
lasts more than one hundred days, with extreme cases beyond several years.

In a given system, such as the Triangulum Galaxy (also known as Messier 33;
hereafter, M33), all its member stars essentially are at the same distance
($D$) from us. Compared to their intrinsic luminosity, their apparent
brightness is scaled by the same factor of $1/D^{2}$. As a result, the
observed PLR of Miras in different galaxies located at various distances are
expected to follow the same functional form, up to a zero-point offset.
Astronomical time-series imaging surveys, using modern CCD cameras on meter-
class or larger telescopes, routinely perform accurate brightness measurements
(called photometry) of millions of stars in other galaxies. These measurements
are usually reported using units of magnitude, which are inversely
proportional to the logarithm of the light flux. Thus, brighter objects have
smaller magnitudes. Repeated imaging of a given field containing a large
number of stars yields time series measurements of stellar brightness called
light curves. To determine the colors of stars, images are usually taken using
filters (also called bands) which only allow certain wavelengths to pass. As a
result, several light curves in various colors are collected for each star;
these constitute the so-called multi-band light curve data.

The third phase of the OGLE survey (Udalski et al., 2008) obtained densely-
sampled light curves of 1,663 Miras (Soszyński et al., 2009) in the Large
Magellanic Cloud (LMC), a satellite galaxy of the Milky Way. The survey was
conducted using two filters ($V$ and $I$). The OGLE-III $I$\- and $V$-band
light curves have a median of about 470 and 45 photometric measurements per
object, respectively. Figure 1(a) shows the light curve obtained by this
survey for one O-rich Mira with a period of $164.84$ days. The horizontal axis
represents time, while the vertical axis shows the magnitude. Note that as
brighter objects have smaller magnitude values, the vertical axis is flipped.
The $I$\- and $V$-band light curves are plotted using black circles and orange
triangles, respectively. The dashed curves are the best sinusoidal fit to each
band separately. Using additional measurements in the NIR $J$, $H$ and $K_{s}$
bands, Yuan et al. (2017) obtained tight PLRs shown in Figure 1(b). In this
figure, each point represents one O-rich Mira. The vertical axis is the
average magnitude of each star, while the horizontal axis is the period of
each object, using a logarithmic scale (with base $10$). The solid lines in
Figure 1(b) are least-squares fits to each band, using a quadratic polynomial
in the logarithm of the period. The zero-order terms of the polynomials are
used for astronomical distance determination.

(a) One LMC O-rich Mira.

(b) PLRs for LMC.

Figure 1: The left panel illustrates observations of the O-rich Mira by the
OGLE survey (ID: OGLE-LMC-LPV-00082). The right panel shows the $J$-, $H$\-
and $K_{s}$-band PLRs for O-rich Miras in the LMC, along with the best-fit
quadratic polynomials.

The OGLE survey has obtained light curves of unusual high quality thanks to
decades-long observations with a dedicated telescope. More typical
astronomical time-series surveys have lower data quality which hinders
accurate period estimation and PLR construction. Yuan et al. (2018) analyzed
Miras in M33, a galaxy in our Local Group located $\sim 20\times$ farther than
the LMC. Each object was measured in four bands ($I$, $J$, $H$, and $K_{s}$)
with a median number of observations of 68, 5, 6, and 11, respectively. The
left panel of Figure 2 shows the light curves of one of these Miras. It is
evident that both the number of data points and their signal-to-noise ratio
are limited. This dataset is not alone in terms of the difficulty to
accurately recover periods. For example, the second data release of the Gaia
survey (Mowlavi et al., 2018) consists of 550,737 variable stars observed in
three filters each with an median of about 25 data points, while the second
data release of the Pan-STARRS survey (Flewelling, 2018), contained an average
of less than 10 measurements in each of the five bands it used.

Figure 2: In the left panel is a typical light curve for a Mira in M33. In the
right panel is resulting periodogram applying the SP method (He et al., 2016)
to its $I$-band light curve. The vertical dashed line indicates the estimated
frequency.

When light curves are sparsely sampled and subject to a high level of noise,
the task of period estimation is no doubt challenging. The Lomb-Scargle method
(Lomb, 1976; Scargle, 1982) is the most commonly used for this purpose. For a
general periodic star, the method fits a simple sinusoidal curve to a single-
band light curve. The fitting is performed over a sequence of trial
frequencies (the inverse of periods), and the frequency with the minimal
residual sum of squares is selected as the estimate. However, a light curve of
one Mira variable is not exactly periodic. From Figure 1(a), it is evident
that either the $I$\- or the $V$-band light curves have additional
fluctuations that are not accounted for by the fitted sinusoids (dashed
curve). The semi-parametric (SP) model of He et al. (2016) addresses the issue
by fitting a Gaussian Process to the residuals. The Lomb-Scargle method has
also been extended by some authors (Vanderplas and Ivezic, 2015; Long et al.,
2016) via exploiting multi-band light curves of a periodic star. As the number
of parameters usually grows with the number of bands, Long et al. (2016)
employed explicit regularization to constrain the model degree of freedom. For
multi-band Mira observations, Yuan et al. (2018) enhanced the SP model and
have achieved the best performance for Mira variables to date.

After all these efforts and explorations, the next bottleneck in raising
estimation accuracy even further has been encountered: existing methods are
susceptible to the so-called “aliasing effect” due to the discrete sampling of
a continuous function. For illustration, we apply the SP model (He et al.,
2016) to the $I$-band light curve shown in the left panel of Figure 2. The
model output is in the right panel of Figure 2, with profile log-likelihood on
the vertical axis and frequency on the horizontal axis. Notice the evaluated
log-likelihood is highly multimodal, and the frequency with maximum log-
likelihood (the primary peak marked by the vertical dashed line) is identified
as the estimate. For a light curve with a weak signal, a false frequency can
show stronger evidence (larger likelihood) than the true one. When this
happens, an incorrect period estimation is acquired. The work of Yuan et al.
(2018) devised an ad-hoc trick for remedy: instead of the primary peak (with
the largest log-likelihood value), they selected the secondary peak (with the
second largest log-likelihood value) as the final estimate if it better fits
the PLR.

In addition to the estimation accuracy, computational cost is another curb on
the application of some existing methods. Due to the multimodal nature of the
likelihood, a grid search over the frequency domain seems inevitable to find
the optimal frequency. This is especially troublesome for the methods of He et
al. (2016) and Yuan et al. (2018). When these methods estimate the period of a
Mira, the related Gaussian process likelihood has to be computed repeatedly
over the grid of trial frequencies. Moreover, this computationally intensive
process has to be carried out for each Mira in the dataset.

In this work, we propose a novel period estimation framework to address these
challenges. Based on the state-of-the-art method of He et al. (2016) and Yuan
et al. (2018), a Bayesian hierarchical model is constructed for Mira
variables. The model embeds the PLR as one of its model components, and
simultaneously models all Mira light curves in a dataset. The resulting model
is capable of automatically reducing the aliasing effect, thereby increasing
estimation accuracy. Under this framework, each Mira supplies information to
update their common PLR. In return, the PLR component guides parameter
estimation of individual Miras. In our numerical studies, this mechanism is
demonstrated to be key for improving estimation accuracy.

To reduce the computational cost, we develop an efficient algorithm for model
inference. Given a Bayesian model, Monte Carlo Markov Chain (MCMC) is the
standard tool for inference. However, computing MCMC will be prohibitively
expensive for our model. As we are fitting a collection of Mira light curves,
tens of thousands of parameters will be involved even for a dataset of
moderate size. Moreover, the multimodal posteriors will also obstruct the
convergence of the Markov chain. Fortunately, advances in variational
inference (VI; Jordan et al., 1999; Blei et al., 2017; Zhang et al., 2019)
bring us the opportunity to exploit the problem in the new framework. VI
converts the Bayesian inference into an optimization problem, and finds the
best density in a pre-specified family to approximate the actual posterior
distribution. To further improve computational efficiency for large scale
problems, we adopt the stochastic variational inference framework by Hoffman
et al. (2013). Using an unbiased estimate of the natural gradient during each
iteration, our algorithm is able to avoid costly numerical integration for the
frequency parameter, and avoid scanning over all samples per iteration.

Under the same simulation as Yuan et al. (2018), we show our approach only
fails to recover the period for $\sim 2$% of the Miras, as opposed to their
$\sim 9$%. The absolute deviation error is remarkably reduced by 77%.
Furthermore, the method of Yuan et al. (2018) costs hundreds of CPU hours,
which is not possible to run without a computer cluster. In contrast, our
method finishes within a few hours on a personal computer.

The benefit of the proposed model goes beyond fast and accurate period
estimation. Firstly, existing methods only deliver a point estimate for the
frequency (or period) without an uncertainty measurement. Our work is the
first attempt to quantify and approximate the multimodal posterior
uncertainty. This is achieved by allowing the variational distribution for
frequency to be any continuous function over a compact interval. An analytical
solution has been derived for its update in our work. Secondly, the PLR is
traditionally determined by a two-step procedure: average magnitudes and
periods are estimated for individual Miras separately; then the PLR is
obtained via a least-squares fit. The proposed method is a more systematic way
to estimate the PLR, which is the fundamental tool for distance measurements.
By adopting Bayesian modeling, we directly harvest its estimate and
uncertainty quantification as direct model products.

The rest of the paper is organized as follows. Our proposed method is
developed in Section 2. Its computational algorithm via the stochastic
variational inference is developed in Section 3. We then conduct simulation
studies to compare our method with existing methods in Section 4. Finally, in
Section 5, the proposed method is applied to a real data set of M33 Miras
(Yuan et al., 2018). Section 6 studies the robustness of the proposed method
to downsampling. The Supplementary Material (He et al., 2020) contains proofs
of theoretical results and additional numerical results.

##  2 The Hierarchical Model

$\displaystyle y_{ibj}|s_{ibj}$
$\displaystyle\sim\mathcal{N}(s_{ibj},\sigma_{ibj}^{2})\text{ for
}i\in\mathcal{I},\ b\in\mathcal{B},\ j\in\mathcal{J}_{ib},$ (1) $\displaystyle
s_{ibj}$
$\displaystyle=m_{ib}+\mathbf{b}_{ibj}^{T}\boldsymbol{\beta}_{ib}+h_{ibj};$
(2) $\displaystyle m_{ib}|f_{i},$
$\displaystyle\boldsymbol{\alpha}_{b},\gamma_{b}\sim\mathcal{N}(\mathbf{d}_{i}^{T}\boldsymbol{\alpha}_{b},1/\gamma_{b})\text{
for }i\in\mathcal{I},\ b\in\mathcal{B};$ (3)
$\displaystyle\boldsymbol{\beta}_{i}|\boldsymbol{\Omega}$
$\displaystyle\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Omega}^{-1})\text{ for
}i\in\mathcal{I};$ (4) $\displaystyle h_{ib}(\cdot)|f_{i}$
$\displaystyle\sim\mathcal{GP}(0,k_{b}(\cdot,\cdot))\text{ for
}i\in\mathcal{I},\ b\in\mathcal{B};$ (5) $\displaystyle f_{i}$
$\displaystyle\sim\mathrm{Unif}(f_{\min},f_{\max})\text{ for
}i\in\mathcal{I};$ (6) $\displaystyle\boldsymbol{\alpha}_{b}$
$\displaystyle\sim\mathcal{N}\left(\bar{\boldsymbol{\alpha}}_{b},(1/\bar{\delta})\mathbf{I}\right)\text{
for }b\in\mathcal{B};$ (7) $\displaystyle\gamma_{b}$
$\displaystyle\sim\text{Gamma}(\bar{\gamma}_{b}\bar{r},\bar{r})\text{ for
}b\in\mathcal{B};$ (8) $\displaystyle\boldsymbol{\Omega}$
$\displaystyle\sim\text{Wishart}(\bar{\boldsymbol{\Omega}}/\bar{n},\bar{n}).$
(9) Figure 3: The hierarchical model for period estimation and PL relation
construction.

Now we elucidate our Bayesian hierarchical model for period estimation with
multi-band Mira data. Our dataset is a collection of $N$ Miras in a galaxy and
the observations of each Mira are collected with $B$ filters (bands) indexed
by $b\in\\{1,2,\cdots,B\\}$. For the $i$-th star $(i=1,\cdots,N)$, the data
collected with the $b$-th filter is denoted as
$\mathcal{D}_{ib}=\\{(t_{ibj},y_{ibj},\sigma_{ibj})\\}_{j=1}^{n_{ib}}$, where
$y_{ibj}$ is the magnitude measured through the $b$-th filter at time
$t_{ibj}$, and the corresponding measurement uncertainty is $\sigma_{ibj}$. As
each Mira is observed through the $b$-th filter, all measurements for the
$i$-th Mira are denoted as $\mathcal{D}_{i}=\cup_{b}\mathcal{D}_{ib}$. The
period of the $i$-th Mira will be denoted as $p_{i}$, and its frequency is
$f_{i}=1/p_{i}$. Note $p_{i}$ (or $f_{i}$) is one of our key quantities of
interest.

In the rest of the paper, our whole collection of Mira data will be denoted as
$\mathcal{D}=\cup_{i}\mathcal{D}_{i}$. The following index sets will be also
used: $\mathcal{I}=\\{1,\cdots,N\\}$, $\mathcal{B}=\\{1,\cdots,B\\}$, and
$\mathcal{J}_{ib}=\\{1,\cdots,n_{ib}\\}$. Our model is based on the work of He
et al. (2016) and Yuan et al. (2018), where each Mira is modeled separately.
We jointly model the whole collection $\mathcal{D}$ of Miras, and
simultaneously estimate the PLR. The overall model is listed in Figure 3 and
its diagram plotted in Figure 4. This is a hierarchical Bayesian model, where
Eqn. (1)–(2) constitute the level-1, Eqn. (3)–(6) constitute the level-2, and
Eqn. (7)–(9) constitute the level-3 of the model. The following subsections
explain the model components in details.

$f_{i}$ $\boldsymbol{\alpha}_{b}$$\bar{\boldsymbol{\alpha}}_{b}$
$\bar{\delta}$$m_{ib}$$\gamma_{b}$ $\bar{\gamma}_{b}$ $\bar{r}$
$\boldsymbol{\beta}_{ib}$ $\boldsymbol{\Omega}$ $\bar{\boldsymbol{\Omega}}$
$\bar{n}$ $y_{ibj}$ $h_{ibj}$
$b\in\mathcal{B}$$j\in\mathcal{J}_{ib}$$i\in\mathcal{I}$$b\in\mathcal{B}$

Figure 4: The hierarchical model diagram.

###  2.1 The light curve decomposition

Denote $\mathcal{D}_{ib}=\\{(t_{ibj},y_{ibj},\sigma_{ibj})\\}_{j=1}^{n_{ib}}$
as the observed light curve for the $i$-th Mira with the $b$-th filter. The
magnitude $y_{ibj}$ is decomposed as signal plus noise terms
$y_{ibj}=s_{ib}(t_{ibj})+\sigma_{ibj}\epsilon_{ibj},$ where $s_{ib}(\cdot)$ is
the signal as a function of time, and $\epsilon_{ibj}$’s are standard Gaussian
noise. In other words, the observed light curve is noisy observation of
$s_{ib}(t)$ at sparse time points. Following He et al. (2016), this signal
function is further decomposed as

$s_{ib}(t)=m_{ib}+\mathbf{b}(f_{i}\cdot
t)^{T}\boldsymbol{\beta}_{ib}+h_{ib}(f_{i}\cdot t),$ (10)

where $m_{ib}$ is the average magnitude, $\mathbf{b}(u)=(\cos(2\pi
u),\sin(2\pi u))^{T}$ is the sinusoidal basis vector and
$\boldsymbol{\beta}_{ib}=(\beta_{ib1},\beta_{ib2})^{T}$ is the corresponding
coefficient, and $h_{ib}(f_{i}\cdot t)$ is a smooth function that captures the
non-periodic variations in Mira magnitudes. Note the sinusoidal function
$\mathbf{b}(f_{i}\cdot t)^{T}\boldsymbol{\beta}_{ib}$ has a period of
$p_{i}=1/f_{i}$. Figure 5 illustrates the decomposition (10) for the $I$-band
light curve in Figure 1(a). In the upper panel, the solid curve represents the
fitted signal $s_{ib}(t)$. The middle panel shows the fitted periodic
component $m_{ib}+\mathbf{b}(f_{i}\cdot t)^{T}\boldsymbol{\beta}_{ib}$, and
the lower panel is for the non-periodic component $m_{ib}+h_{ib}(f_{i}\cdot
t)$.

With simplified notations $s_{ibj}=s_{ib}(t_{ibj})$,
$\mathbf{b}_{ibj}=\mathbf{b}(f_{i}\cdot t_{ibj})$ and
$h_{ibj}=h_{ib}(f_{i}\cdot t_{ibj})$, Eqn. (10) constitutes the level-1 of the
hierarchical model in Figure 3. This top level models the individual observed
light curve $\mathcal{D}_{ib}$ separately. Our joint modeling comes in by
introducing the PLR model component in the next subsection.

Figure 5: An example of Mira light curve decomposition. The solid curve in the
upper panel is the fitted light curve. The solid curve in the middle panel is
the exact periodic component. The lower panel is the non-periodic function.

######  Remark 1.

The Generalized Lomb-Scargle method (Zechmeister and Kürster, 2009) and its
multi-band version only have the average magnitude and the sinusoidal
component in the decomposition (10). The additional term $h_{ib}(\cdot)$ is
introduced by He et al. (2016) and Yuan et al. (2018) to account for the
additional variability of Mira light curves. All of these methods model each
Mira separately, without the hierarchical information introduced in the
following subsections.

###  2.2 The Period-Luminosity relation

For a specific band $b\in\mathcal{B}$, the average magnitude $m_{ib}$ and the
period $p_{i}$ follow the PLR for all $i\in\mathcal{I}$. As in Figure 1(b),
for each band, the PLR can be fitted by a quadratic polynomial in
$\log_{10}(p)$. Since our model is parameterized by frequency, the following
equivalent PLR (noting $\log_{10}p_{i}=-\log_{10}f_{i}$) is employed,

$m_{ib}=\alpha_{b3}\log^{2}_{10}f_{i}+\alpha_{b2}\log_{10}f_{i}+\alpha_{b1}+e_{ib}=\boldsymbol{\alpha}_{b}^{T}\mathbf{d}_{i}+e_{ib}\,.$
(11)

In the above, $\mathbf{d}_{i}=(1,\log_{10}(f_{i}),\log^{2}_{10}(f_{i}))^{T}$
is the PLR basis vector for the $i$-th Mira,
$\boldsymbol{\alpha}_{b}=(\alpha_{b1},\alpha_{b2},\alpha_{b3})^{T}$ is the PLR
coefficient for the $b$-th band, and the PLR residual $e_{ib}$ is assumed to
follow $\mathcal{N}(0,1/\gamma_{b})$ with precision $\gamma_{b}$. In summary,
given $f_{i},\boldsymbol{\alpha}_{b},\gamma_{b}$, the average magnitude is
dictated by
$m_{ib}|f_{i},\boldsymbol{\alpha}_{b},\gamma_{b}\sim\mathcal{N}(\boldsymbol{\alpha}_{b}^{T}\mathbf{d}_{i},1/\gamma_{b}).$
In addition, we place on $f_{i}$ a uniform prior over a compact interval
$[f_{\min},f_{\max}]$, i.e., $f_{i}\sim\mathrm{Unif}(f_{\min},f_{\max})$. The
compact interval $[f_{\min},f_{\max}]$ can be selected as the search range of
interest. For most Miras, their frequency lies inside the interval
$[1/1000,1/100]$. These constitutes parts of the level-2 of our model in
Figure 3.We remark that while the level-1 of the model treats each light curve
separately, the level-2 dictates the common a-priori information shared by
individual Mira parameters.

In our proposed model, the coefficient $\boldsymbol{\alpha}_{b}$ and the
precision $\gamma_{b}$ are also unknown parameters with conjugate prior
$\boldsymbol{\alpha}_{b}\sim\mathcal{N}(\bar{\boldsymbol{\alpha}}_{b},(1/\bar{\delta})\mathbf{I})$,
$\gamma_{b}\sim\text{Gamma}(\bar{\gamma}_{b}\bar{r},\bar{r})$, where
$\bar{\boldsymbol{\alpha}}_{b}$, $\bar{\gamma}_{b}$, $\bar{\delta}$, $\bar{r}$
are fixed hyper-parameters. These become Eqn. (7)–(8), and part of level-3 of
our model in Figure 3. Notice that the level-3 parameters
$\boldsymbol{\alpha}_{b},\gamma_{b}$ and $\boldsymbol{\Omega}$ (to be
introduced in the next subsection) govern the probability distribution of the
Mira samples.

######  Remark 2.

Eqn. (3) plays a key role in our hierarchical model. It governs the relation
between $m_{ib}$ and $f_{i}$, and discourages implausible parameter values.
This model component is the key to reduce the aliasing effect and improves the
period estimation accuracy. Moreover, a direct product of our proposed method
is the PLR coefficients $\boldsymbol{\alpha}_{b}$’s, which is one of the
important pursuits for modern astronomical surveys.

###  2.3 The periodic component

For the $i$-th star, let
$\boldsymbol{\beta}_{i}=(\boldsymbol{\beta}_{i1}^{T},\boldsymbol{\beta}_{i2}^{T},\cdots,\boldsymbol{\beta}_{iB}^{T})^{T}$
be the collection of periodic coefficients across all bands. This vector has
multivariate Gaussian prior
$\boldsymbol{\beta}_{i}|\boldsymbol{\Omega}\sim\mathcal{N}(\mathbf{0},\boldsymbol{\Omega}^{-1})$
in Eqn. (4) of Figure 3. The precision matrix $\boldsymbol{\Omega}$ is also
unknown with a Wishart prior
$\boldsymbol{\Omega}\sim\text{Wishart}(\bar{\boldsymbol{\Omega}}/\bar{n},\bar{n})$
in Eqn. (9), for fixed hyper-parameters $\bar{\boldsymbol{\Omega}}$ and
$\bar{n}$.

The multivariate Gaussian prior (4) on $\boldsymbol{\beta}_{i}$ empowers us
the ability to modeling the correlation between the periodic components across
all bands. Long et al. (2016) and Yuan et al. (2018) have noticed both
periodic amplitude and phase are linearly correlated across bands. The prior
distribution as given in (4) on $\boldsymbol{\beta}_{i}$ allows us to take
this correlation into account. Nevertheless, the assumption imposed by Eqn.
(4) is weaker than the explicit regularization of Long et al. (2016).

###  2.4 The non-periodic component

Recall the decomposition (10) has an additional term $h_{ib}(\cdot)$
accounting for the non-periodic variation. For each $j\in\mathcal{J}_{ib}$, we
denote $u_{ibj}=f_{i}\cdot t_{ibj}$ as the phase, which is the frequency
$f_{i}$ multiplied by the time points $t_{ibj}$. Note $h_{ib}(\cdot)$ is a
function of $u_{ibj}$’s. For the a specific band $b\in\mathcal{B}$, the random
function $h_{ib}(\cdot)$ is modeled by a Gaussian process (GP, Rasmussen and
Williams, 2005) with zero mean and a kernel function $k_{b}(\cdot,\cdot)$.
Under the GP assumption, the vector
$\mathbf{h}_{ib}=\left(h_{ib}(u_{ib1}),\cdots,h_{ib}(u_{ibn_{ib}})\right)^{T}$
follows a multivariate Gaussian distribution with zero mean and covariance
matrix
$\mathbf{H}_{ib}=\big{[}k_{b}(u_{ibp},\,u_{ibq})\big{]}_{p,q=1}^{n_{ib}}$. The
Eqn. (5) in Figure 3 is exactly this GP component. Note the $b$-th band data
of all Miras share a common kernel $k_{b}(\cdot,\cdot)$. In this work, we
employed the following kernel function,

$k_{b}(u,u^{\prime})=\tau_{b1}\cdot\exp\big{\\{}-(u-u^{\prime})^{2}/\tau_{b2}\big{\\}}+\tau_{b3}\cdot\mathrm{I}(u=u^{\prime}),$

where $\mathrm{I}(\cdot)$ is the indicator function, and
$\boldsymbol{\tau}_{b}=(\tau_{b1},\tau_{b2},\tau_{b3})^{T}$ is the kernel
parameter for the $b$-th band. The first term of the above kernel is the
squared exponential kernel with variance parameter $\tau_{b1}$ and length-
scale parameter $\tau_{b2}$. The second term is a nugget term accounting for
additional variability. Our adoption of the squared exponential kernel follows
a common practice in the kernel machine field, though other kernel functions
can also be used. The Gaussian process curve fitting is generally not
sensitive to the choice of kernel.

In this formulation, for each band, the non-period components of light curves
of all Miras are assumed to be realizations of the same Gaussian process. This
allows Miras with sparse observations or relatively low data quality to
“borrow strength” from one another when fitting individual light curves.
Letting all Miras share the same set of kernel parameters has the additional
advantage of avoiding the intensive optimization of kernel parameters for
individual light curves, and thereby greatly speeding up computation.

######  Remark 3.

In the work of He et al. (2016) and Yuan et al. (2018), the term
$h_{ib}(\cdot)$ is chosen as a function of time $t$, and the kernel parameter
is distinct for different bands and different Miras. In contrast, our
$h_{ib}(\cdot)$ is a function of phase $u=f_{i}\cdot t$, and the kernel
parameter is shared by all Miras in a given band. Our choice of using
$h_{ib}(f_{i}\cdot t)$ makes the function argument of the non-periodic
component consistent with the periodic component $\mathbf{b}(f_{i}\cdot
t)^{T}\boldsymbol{\beta}_{ib}$ as shown in (10). This choice is supported by
two numerical results presented in Section S.3 of the Supplementary Material
(He et al., 2020).

###  2.5 The joint likelihood

All of the above leads to our model listed in Figure 3, which has a few
unknown parameters. In the rest of the paper, to simplify the notation, the
PLR coefficients for all bands are collectively denoted as
$\boldsymbol{\alpha}=\\{\boldsymbol{\alpha}_{1},\cdots,\boldsymbol{\alpha}_{B}\\}$,
and all precision parameters as
$\boldsymbol{\gamma}=(\gamma_{1},\cdots,\gamma_{B})^{T}$. For the $i$-th Mira,
its average magnitudes across all bands are denoted by
$\mathbf{m}_{i}=(m_{i1},\cdots,m_{iB})^{T}$. For the whole collection of
Miras, all their average magnitude parameters are denoted by
$\mathbf{m}=\\{\mathbf{m}_{1},\cdots,\mathbf{m}_{N}\\}$. In addition, we set
$\mathbf{f}=\\{f_{1},\cdots,f_{N}\\}$ for the frequency parameters of all
Miras and
$\boldsymbol{\beta}=\\{\boldsymbol{\beta}_{1},\cdots,\boldsymbol{\beta}_{N}\\}$
for all periodic coefficients.

Given the hierarchical model in Figure 3, Bayesian inference requires
computing the posterior of parameters, which is proportional to the joint
distribution of all parameters. For the data $\mathcal{D}_{i}$ of the $i$-th
Mira, its related joint likelihood (of Eqn. (1) through (6)) is

$\displaystyle\prod_{b=1}^{B}\prod_{j=1}^{n_{ib}}\exp\Big{\\{}-\frac{1}{2\sigma_{ibj}^{2}}(y_{ibj}-m_{ib}-\mathbf{b}_{ibj}^{T}\boldsymbol{\beta}_{ib}-h_{ibj})^{2}\Big{\\}}\times\prod_{b=1}^{B}|\mathbf{H}_{ib}|^{-1/2}\exp\Big{\\{}-\frac{1}{2}\mathbf{h}_{ib}^{T}\mathbf{H}_{ib}^{-1}\mathbf{h}_{ib}\Big{\\}}$
$\displaystyle\qquad\qquad\times|\boldsymbol{\Omega}|^{1/2}\exp\Big{\\{}-\frac{1}{2}\boldsymbol{\beta}_{i}^{T}\boldsymbol{\Omega}\boldsymbol{\beta}_{i}\Big{\\}}\times\prod_{b=1}^{B}\gamma_{b}^{1/2}\exp\Big{\\{}-\frac{\gamma_{b}}{2}(m_{ib}-\mathbf{d}_{i}^{T}\boldsymbol{\alpha}_{b})^{2}\Big{\\}}.$

The Gaussian process component $\mathbf{h}_{ib}$ can be integrated out. In
fact, given $m_{ib},\boldsymbol{\beta}_{ib},f_{i}$, the vector
$\mathbf{y}_{ib}=(y_{ib1},\cdots,y_{ibn_{ib}})^{T}$ follows a Gaussian
distribution with mean vector
$m_{ib}\mathbf{1}+\mathbf{B}_{ib}\boldsymbol{\beta}_{ib}$ and covariance
matrix
$\boldsymbol{\Sigma}_{ib}=\mathbf{H}_{ib}+\text{diag}(\sigma_{ib1}^{2},\cdots,\sigma_{ibn_{ib}}^{2})$,
where $\mathbf{1}\in\mathbb{R}^{n_{ib}}$ is a vector of ones and
$\mathbf{B}_{ib}=\big{(}\mathbf{b}_{ib1},\cdots,\mathbf{b}_{ibn_{ib}}\big{)}^{T}\in\mathbb{R}^{n_{ib}\times
2}$. Consequently, the likelihood for the observation of the $i$-th Mira
simplifies to

$\displaystyle L_{i}=$
$\displaystyle\prod_{b=1}^{B}|\boldsymbol{\Sigma}_{ib}|^{-1/2}\exp\Big{\\{}-\frac{1}{2}(\mathbf{y}_{ib}-m_{ib}\mathbf{1}-\mathbf{B}_{ib}\boldsymbol{\beta}_{ib})^{T}\boldsymbol{\Sigma}_{ib}^{-1}(\mathbf{y}_{ib}-m_{ib}\mathbf{1}-\mathbf{B}_{ib}\boldsymbol{\beta}_{ib})\Big{\\}}$
$\displaystyle\qquad\times|\boldsymbol{\Omega}|^{1/2}\exp\Big{\\{}-\frac{1}{2}\boldsymbol{\beta}_{i}^{T}\boldsymbol{\Omega}\boldsymbol{\beta}_{i}\Big{\\}}\times\prod_{b=1}^{B}\gamma_{b}^{1/2}\exp\Big{\\{}-\frac{\gamma_{b}}{2}(m_{ib}-\mathbf{d}_{i}^{T}\boldsymbol{\alpha}_{b})^{2}\Big{\\}}.$
(12)

Given the whole dataset $\mathcal{D}$, the overall joint distribution (of Eqn.
(1) through (9)) is

$\displaystyle
p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D})=$
$\displaystyle\Big{(}\prod_{i=1}^{N}L_{i}\Big{)}\times\exp\Big{\\{}-\frac{1}{2}\langle\bar{n}\bar{\boldsymbol{\Omega}}^{-1},\boldsymbol{\Omega}\rangle+\frac{\bar{n}-2B-1}{2}\log\det\boldsymbol{\Omega}\Big{\\}}$
$\displaystyle\quad\times\prod_{b=1}^{B}\exp\Big{\\{}-\bar{r}\gamma_{b}+(\bar{r}\bar{\gamma}_{b}-1)\log\gamma_{b}-\frac{\bar{\delta}}{2}\|\boldsymbol{\alpha}_{b}-\bar{\boldsymbol{\alpha}}_{b}\|^{2}\Big{\\}}.$
(13)

The right hand side is the product of the likelihoods for all observations and
the priors of all related parameters. We know the posterior
$p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha}|\,\mathcal{D})$
is proportional to the above distribution up to an unknown constant.
Traditionally, inference for this hierarchical Bayesian model is carried out
by the Monte Carlo Markov Chain (MCMC). However, it is well known that the
traditional MCMC lacks scalability to a large data set. To make it
computationally feasible to apply our model to real Mira light-curve datasets,
we will develop an efficient algorithm based on the stochastic variational
inference in the next section.

######  Remark 4.

Before fitting the model, we need to specify the hyper-parameters
$\bar{\boldsymbol{\alpha}}_{b}$, $\bar{\delta}$, $\bar{\gamma}_{b}$,
$\bar{r}$, $\bar{\boldsymbol{\Omega}}$, $\bar{n}$. We can set relative small
values for $\bar{r}$, $\bar{n}$, and large values for $\bar{\delta}$. This
makes the priors less informative. The other hyper-parameters
($\bar{\boldsymbol{\Omega}}$ and $\bar{\delta}$) can be easily specified based
on information from existing high-quality astronomical surveys. This is
applicable to the first and second order parameters $\bar{\alpha}_{b2}$,
$\bar{\alpha}_{b3}$ of PLRs, as the shape of PLRs is the same for Miras in
different galaxies. The only exception is the intercept terms
$\bar{\alpha}_{b1}$ for the PLR prior, as the distance for a new dataset is
unknown. Fortunately, for a new survey, we only need to set a roughly reliable
value of $\bar{\alpha}_{b1}$. For this purpose, we can apply multi-band
generalized Lomb-Scargle (MGLS) method to a small subset of the data, and then
use a robust method (e.g. minimal absolute deviation) to get a fitted value of
$\bar{\alpha}_{b1}$.

It is known in the literature that Bayesian hierarchical model may perform
poorly when vague priors are used. Gelman (2006) observed that choice of “non-
informative” prior distributions on hierarchical variance parameters may make
a big impact on the inference, and improper choice may lead to inaccurate
posterior distributions. In our numerical studies, by setting non-extreme
values $\bar{r}=1$, $\bar{n}=1$ for the variance parameters, our prior
distributions do not have the adverse shrinkage effect pointed out in Gelman
(2006) and our model shows accurate estimation results. Moreover, we have
sufficient data so that the posterior distribution is dominated by the
information provided by the data.

######  Remark 5.

We also need to specify a value for the kernel parameter $\boldsymbol{\tau}$.
This is achieved by applying the MGLS method to a small subset of the dataset.
With the estimated frequency from MGLS, we fit sinusoidal model to the light
curves and obtain the sinusoidal residuals. After that, we fit a Gaussian
process to the residual and obtain maximum likelihood estimates of the kernel
parameters.

##  3 Main Algorithm

Variational inference (Jordan et al., 1999; Wainwright et al., 2008) aims at
approximating the true posterior by finding an optimal density function
$q(\cdot)$ among a specific distribution family $\mathcal{Q}$. This $q(\cdot)$
is then employed for subsequent inference. To find the optimal approximation,
the commonly used strategy is to minimize the Kullback-Leibler divergence (KL
divergence) between $q(\cdot)$ and the actual posterior distribution
$p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha}|\,\mathcal{D})$.
In our context, we are trying to solve

$\min_{q(\cdot)\in\mathcal{Q}}\int
q(\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\Omega},\mathbf{f},\mathbf{m},\boldsymbol{\beta})\log\frac{q(\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\Omega},\mathbf{f},\mathbf{m},\boldsymbol{\beta})}{p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha}|\,\mathcal{D})}\,\mathrm{d}\mathbf{m}\,\mathrm{d}\boldsymbol{\beta}\,\mathrm{d}\mathbf{f}\,\mathrm{d}\boldsymbol{\Omega}\,\mathrm{d}\boldsymbol{\gamma}\,\mathrm{d}\boldsymbol{\alpha}$
(14)

for $q(\cdot)$ among a specific distribution family $\mathcal{Q}$.
Equivalently, the variational inference aims at maximizing the evidence lower
bound (ELBO)
$\underline{L}(\mathcal{D}):=\mathbb{E}_{q}\big{\\{}\log\frac{p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D})}{q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})}\big{\\}}$,
where the expectation is taken with respect to
$q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})$.
By Jensen’s inequality, it holds that

$\displaystyle\log L(\mathcal{D})$ $\displaystyle=\log\Big{[}\int
p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D})\,\,\mathrm{d}\mathbf{m}\,\mathrm{d}\boldsymbol{\beta}\,\mathrm{d}\mathbf{f}\,\mathrm{d}\boldsymbol{\Omega}\,\mathrm{d}\boldsymbol{\gamma}\,\mathrm{d}\boldsymbol{\alpha}\Big{]}$
$\displaystyle\geq\int
q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})\log\frac{p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D})}{q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})}\,\,\mathrm{d}\mathbf{m}\,\mathrm{d}\boldsymbol{\beta}\,\mathrm{d}\mathbf{f}\,\mathrm{d}\boldsymbol{\Omega}\,\mathrm{d}\boldsymbol{\gamma}\,\mathrm{d}\boldsymbol{\alpha}=\underline{L}(\mathcal{D})\,.$

The term $\log L(\mathcal{D})$ is the marginal likelihood for the dataset, and
it is interpreted as the model evidence. The above implies the ELBO
$\underline{L}(\mathcal{D})$ acts as a lower bound for $\log L(\mathcal{D})$.

The variational inference problem (14) requires specifying the distribution
family $\mathcal{Q}$. This family should be large enough to approximate the
true posterior accurately, while not too complicated for ease of computation.
The mean field family (Blei et al., 2017) is a popular choice, for which the
distribution $q(\cdot)$ factors as a product of distributions for each
parameter. In this work, we build a structured mean-field variational family
$\mathcal{Q}$. The densities inside $\mathcal{Q}$ have the following form,

$\displaystyle
q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})=$
$\displaystyle\prod_{i=1}^{N}q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})\times
q(\boldsymbol{\alpha})q(\boldsymbol{\gamma})\times q(\boldsymbol{\Omega})$
$\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\left[q(\mathbf{m}_{i},\boldsymbol{\beta}_{i}|f_{i})q(f_{i})\right]\times\prod_{b=1}^{B}\left[q(\boldsymbol{\alpha}_{b})q(\gamma_{b})\right]\times
q(\boldsymbol{\Omega})\,.$ (15)

The distribution specification for each factor of Eqn. (15) is listed in Table
1. Most of them belong to the exponential family. Section S.1 in the
Supplementary Material (He et al., 2020) provides a review on these
exponential family distributions and our employed notions for their canonical
parameters.

For most parameters (e.g.
$\boldsymbol{\alpha}_{b},\gamma_{b},\boldsymbol{\Omega}$) of our model, their
full conditional distributions belong to the exponential distribution family.
Their variational distributions are assumed to be their conjugate
counterparts. For example, from the joint distribution (13), we can find the
full conditional for $\gamma_{b}$ is a Gamma distribution. Correspondingly,
the variational distribution $q(\gamma_{b})$ is assumed to be a Gamma
distribution.

The variational distribution
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})=q(\mathbf{m}_{i},\boldsymbol{\beta}_{i}|f_{i})q(f_{i})$
of each Mira is constructed with special consideration. Recall the actual
posterior is multimodal in $f_{i}$. To facilitate uncertainty quantification
of the estimated frequency, the variational distribution $q(f_{i})$ is allowed
to be any continuous density on the compact interval
$[f_{\mathrm{min}},f_{\mathrm{max}}]$. Besides, for convenience, the
distribution $q(\mathbf{m}_{i},\boldsymbol{\beta}_{i}|f_{i})$ is assummed to
be Gaussian and depends on the frequency $f_{i}$. This conditional dependence
is necessary, because the amplitude and phase of the fitted sinusoidal curve
vary with $f_{i}$.

Table 1: Notation for the mean field distributions with their canonical
parameters
Local | Global
---|---
Factor | Distribution | Parameter | Factor | Distribution | Parameter
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i}|f_{i})$ | Gaussian | $\eta_{i}^{\theta}=(\eta_{i}^{\theta 1},\eta_{i}^{\theta 2})$ | $q(\boldsymbol{\alpha}_{b})$ | Gaussian | $\eta_{b}^{\alpha}=(\eta_{b}^{\alpha 1},\eta_{b}^{\alpha 2})$
$q(f_{i})$ | Free Form | None | $q(\gamma_{b})$ | Gamma | $\eta_{b}^{\gamma}=(\eta_{b}^{\gamma 1},\eta_{b}^{\gamma 2})$
|  |  | $q(\boldsymbol{\Omega})$ | Wishart | $\eta^{\boldsymbol{\Omega}}=(\eta^{\boldsymbol{\Omega}1},\eta^{\boldsymbol{\Omega}2})$

After specifying $\mathcal{Q}$, the model inference is completely framed as an
optimization problem. Our target is maximizing the evidence lower bound, i.e.
$\max_{\ q(\cdot)\in\mathcal{Q}}\,\underline{L}(\mathcal{D})$, with
$\mathcal{Q}$ determined in (15). The optimal solution
$q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})$
will then be used for model inference, in place of the actual posterior
$p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha}|\,\mathcal{D})$.
This optimization is typically carried out by the coordinate ascent algorithm
(Blei et al., 2017), where one factor in (15) gets updated with the others
fixed during each round of iteration. Though this updating strategy improves
computational speed over MCMC in lots of cases, it is still not fast enough
for our problem at hand. For efficient computation, we will exploit our
hierarchical model structure under the framework of stochastic variational
inference (SVI, Hoffman et al., 2013). According to the SVI scheme, the
unknown parameters can be divided into two groups: local parameters and global
parameters. The local parameters are specific to each Mira. For the $i$-th
Mira, its own local parameters include the average magnitude $\mathbf{m}_{i}$,
the periodic coefficient $\boldsymbol{\beta}_{i}$, and its frequency $f_{i}$.
On the other hand, the global parameters govern all Mira samples. These global
parameters include the PLR coefficients $\boldsymbol{\alpha}$, the PLR
residual precision vector $\boldsymbol{\gamma}$, and the precision matrix
$\boldsymbol{\Omega}$ for the sinusoidal coefficients. By exploiting the model
structure, SVI becomes much more efficient than the classical coordinate
ascent algorithm. In fact, conditional on the global parameters, the local
parameters of each Mira can be updated independently of each other. On the
other hand, the global parameters can be updated based on mini-batch samples.
The details to update local and global parameters will be provided in Section
3.1 and Section 3.2, respectively.

###  3.1 Update local parameters

In this subsection, we derive analytic formulae to update the local
variational distribution
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})=q(\mathbf{m}_{i},\boldsymbol{\beta}_{i}|f_{i})q(f_{i})$
for each $i\in\mathcal{I}$, while the global $q(\boldsymbol{\alpha})$,
$q(\boldsymbol{\Omega})$ and $q(\boldsymbol{\gamma})$ are fixed. For one
specific $i\in\mathcal{I}$, we need to maximize the ELBO with respect to
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})$. Note the ELBO can be
factored as

$\displaystyle\underline{L}(\mathcal{D})=$
$\displaystyle\mathbb{E}_{q}\Big{\\{}\log\frac{p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D})}{q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})}\Big{\\}}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{N}\mathbb{E}_{q}\Big{\\{}\log\frac{p(\mathbf{m}_{j},\boldsymbol{\beta}_{j},f_{j}|\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D}_{j})}{q(\mathbf{m}_{j},\boldsymbol{\beta}_{j},f_{j})}\Big{\\}}+\mathbb{E}_{q}\Big{\\{}\log\frac{\prod_{b=1}^{B}\left[p(\boldsymbol{\alpha}_{b})p(\gamma_{b})\right]\times
p(\boldsymbol{\Omega})}{\prod_{b=1}^{B}\left[q(\boldsymbol{\alpha}_{b})q(\gamma_{b})\right]\times
q(\boldsymbol{\Omega})}\Big{\\}}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{q}\Big{\\{}\log\frac{p(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i}|\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D}_{i})}{q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})}\Big{\\}}+C\,.$
(16)

In the above, the second equation uses the conditional independence between
observations, and the last term $C$ is a constant irrelevant to the update of
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})$. According to the last line
(16), we only need to consider the data $\mathcal{D}_{i}$ of the $i$-th Mira.
This means given the global $q(\boldsymbol{\alpha})$, $q(\boldsymbol{\Omega})$
and $q(\boldsymbol{\gamma})$, the local variational distribution
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})$ can be updated independently
of the other Miras.

The average magnitude $\mathbf{m}_{i}$ and the sinusoidal coefficient
$\boldsymbol{\beta}_{i}$ can be conveniently updated as a whole. For this
purpose, we introduce a few notations. Let
$\boldsymbol{\theta}_{ib}=(m_{ib},\boldsymbol{\beta}_{ib}^{T})^{T}$ be the
parameter for the $b$-th band. Then, we arrange $\boldsymbol{\theta}_{ib}$’s
of all bands into a vector as
$\boldsymbol{\theta}_{i}=(\boldsymbol{\theta}_{i1}^{T},\cdots,\boldsymbol{\theta}_{iB}^{T})^{T}\in\mathbb{R}^{3B}$.
Correspondingly, two index sets can be defined as
$\mathcal{K}=\\{1,4,7,\cdots,3B-2\\}$ and its complementary
$\mathcal{L}=\\{1,2,\cdots,3B\\}\backslash\mathcal{K}$. The index set
$\mathcal{K}$ records the locations of $m_{ib}$’s in the vector
$\boldsymbol{\theta}_{i}$, while $\mathcal{L}$ corresponds to the locations of
$\boldsymbol{\beta}_{ib}$’s inside $\boldsymbol{\theta}_{i}$.

Recall our target of updating the local
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})=q(\mathbf{m}_{i},\boldsymbol{\beta}_{i}|f_{i})q(f_{i})$
with fixed $q(\boldsymbol{\alpha})$, $q(\boldsymbol{\Omega})$ and
$q(\boldsymbol{\gamma})$. With the definition of $\boldsymbol{\theta}_{i}$, it
is equivalent to state the update in terms of
$q(\boldsymbol{\theta}_{i},f_{i})=q(\boldsymbol{\theta}_{i}|f_{i})q(f_{i})$.
According to our distribution specification in Table 1, $q(f_{i})$ is allowed
to be any continuous density on a compact interval
$[f_{\mathrm{min}},f_{\mathrm{max}}]$, while
$q(\boldsymbol{\theta}_{i}|f_{i})\propto\exp\Big{\\{}\big{\langle}\eta^{\theta
1}_{i},\boldsymbol{\theta}_{i}\boldsymbol{\theta}_{i}^{T}\big{\rangle}+\big{\langle}\eta^{\theta
2}_{i},\boldsymbol{\theta}_{i}\big{\rangle}\Big{\\}}$ is assumed to be
Gaussian with canonical parameters $\eta^{\theta 1}_{i}\in\mathbb{R}^{3B\times
3B}$ and $\eta^{\theta 2}_{i}\in\mathbb{R}^{3B}$. Note the value of
$\eta^{\theta 1}_{i}$ and $\eta^{\theta 2}_{i}$ depends on the frequency
$f_{i}$. With these notations, our goal boils down to finding the optimal
$q(f_{i})$, and for each $f_{i}$ computing the optimal $\eta^{\theta
1}_{i},\eta^{\theta 2}_{i}$. The optimal solutions should maximize the ELBO
(16), which is now equivalently expressed as

$\displaystyle\underline{L}(\mathcal{D})=$
$\displaystyle\mathbb{E}_{q}\Big{\\{}\log\frac{p(\boldsymbol{\theta}_{i},f_{i}|\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D}_{i})}{q(\boldsymbol{\theta}_{i},f_{i})}\Big{\\}}+C\,.$
(17)

This equivalent ELBO involves the conditional likelihood
$p(\boldsymbol{\theta}_{i},f_{i}|\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D}_{i})$,
which is proportional to the likelihood in (12). To convert the expression in
a compact form, we need to define a few more notations. According to (3) and
(4) of our model, the conditional distribution of $\boldsymbol{\theta}_{i}$
given $(f_{i},\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\Omega})$ is
$\mathcal{N}(\boldsymbol{\theta}_{0},\boldsymbol{\Theta}^{-1})$, with its mean
vector $\boldsymbol{\theta}_{0}$ and precision matrix $\boldsymbol{\Theta}$
defined as follows. The sub-vectors of $\boldsymbol{\theta}_{0}$ indexed by
$\mathcal{K}$ and $\mathcal{L}$ are
$\boldsymbol{\theta}_{0}[\mathcal{K}]=(\mathbf{d}_{i}^{T}\boldsymbol{\alpha}_{1},\cdots,\mathbf{d}_{i}^{T}\boldsymbol{\alpha}_{B})^{T}$
and $\boldsymbol{\theta}_{0}[\mathcal{L}]=\mathbf{0}$, respectively. In fact,
$\boldsymbol{\theta}_{0}[\mathcal{K}]$ contains the a-priori mean of
$\mathbf{m}_{i}$ in (3), while $\boldsymbol{\theta}_{0}[\mathcal{L}]$ is the
a-priori mean of $\boldsymbol{\beta}_{i}$ in (4). Similarly, the sub-matrices
of $\boldsymbol{\Theta}$ indexed by $\mathcal{K}$ and $\mathcal{L}$ are

$\boldsymbol{\Theta}[\mathcal{K},\mathcal{K}]=\mathrm{diag}(\gamma_{1},\cdots,\gamma_{B}),\quad\boldsymbol{\Theta}[\mathcal{L},\mathcal{L}]=\boldsymbol{\Omega},\quad\text{and}\quad\boldsymbol{\Theta}[\mathcal{K},\mathcal{L}]=\big{(}\boldsymbol{\Theta}[\mathcal{L},\mathcal{K}]\big{)}^{T}=\mathbf{0}\,,$

respectively. In the above, for example,
$\boldsymbol{\Theta}[\mathcal{K},\mathcal{L}]$ is the submatrix of
$\boldsymbol{\Theta}$ with rows indexed by $\mathcal{K}$, and columns indexed
by $\mathcal{L}$. All of the above is simply writing (3) and (4) of all bands
in a compact form, i.e.,
$\boldsymbol{\theta}_{i}|(f_{i},\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\Omega})\sim\mathcal{N}(\boldsymbol{\theta}_{0},\boldsymbol{\Theta}^{-1})$.
In addition, for the $b$-th band of the $i$-th Mira, define a basis matrix
$\mathbf{C}_{ib}=(\mathbf{1},\mathbf{B}_{ib})\in\mathbb{R}^{n_{ib}\times 3}$.
Then, by rewriting (12) with these new notations, it follows that

$\displaystyle\log
p(\boldsymbol{\theta}_{i},f_{i}|\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D}_{i})=$
$\displaystyle-\frac{1}{2}\sum_{b=1}^{B}(\mathbf{y}_{ib}-\mathbf{C}_{ib}\boldsymbol{\theta}_{ib})^{T}\boldsymbol{\Sigma}_{ib}^{-1}(\mathbf{y}_{ib}-\mathbf{C}_{ib}\boldsymbol{\theta}_{ib})-\frac{1}{2}\log\det\boldsymbol{\Sigma}_{ib}$
$\displaystyle\quad\qquad-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0})^{T}\boldsymbol{\Theta}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0})+C^{\prime},$

with some irrelevant constants $C^{\prime}$. The above expression will be
plugged into (17) to find the optimal solution. Our updating formulae for the
local parameter is summarized by the following theorem, whose proof is
provided in Section S.2 of the Supplementary Material (He et al., 2020).

######  Theorem 3.1.

Suppose the data $\mathcal{D}_{i}$ of the $i$-th Mira and
$q(\boldsymbol{\alpha})$, $q(\boldsymbol{\gamma})$, and
$q(\boldsymbol{\Omega})$ are fixed. Then, for each $f_{i}$, the optimal
$q(\boldsymbol{\theta}_{i}|f_{i})$ maximizing ELBO
$\underline{L}(\mathcal{D})$ has canonical parameters

$\eta^{\theta
1}_{i}=\mathbb{E}_{q(\boldsymbol{\gamma})}\mathbb{E}_{q(\boldsymbol{\Omega})}(\xi^{\theta
1}_{i}),\quad\text{and}\quad\eta^{\theta
2}_{i}=\mathbb{E}_{q(\boldsymbol{\alpha})}(\xi^{\theta 2}_{i}).$ (18)

In the above, $\xi^{\theta
1}_{i}=-\frac{1}{2}\times\mathrm{diag}\left(\mathbf{C}_{i1}^{T}\boldsymbol{\Sigma}_{i1}^{-1}\mathbf{C}_{i1},\cdots,\mathbf{C}_{iB}^{T}\boldsymbol{\Sigma}_{iB}^{-1}\mathbf{C}_{iB}\right)-\frac{1}{2}\boldsymbol{\Theta}$
and

$\xi^{\theta
2}_{i}=\left((\mathbf{C}_{i1}^{T}\boldsymbol{\Sigma}_{i1}^{-1}\mathbf{y}_{i1})^{T},\cdots,(\mathbf{C}_{iB}^{T}\boldsymbol{\Sigma}_{iB}^{-1}\mathbf{y}_{iB})^{T}\right)^{T}+\boldsymbol{\Theta}\boldsymbol{\theta}_{0}\,.$

Moreover, the optimal $q(f_{i})$ is given by

$q(f_{i})\propto\exp\Big{\\{}-\frac{1}{2}g(f_{i})-\frac{1}{2}\log\det(-\eta_{i}^{\theta
1})\Big{\\}},$ (19)

with

$\displaystyle g(f_{i})=\frac{1}{2}(\eta_{i}^{\theta 2})^{T}(\eta_{i}^{\theta
1})^{-1}\eta_{i}^{\theta
2}+\left\langle\mathbb{E}_{q(\boldsymbol{\Omega})}\mathbb{E}_{q(\boldsymbol{\gamma})}\boldsymbol{\Theta},\
\mathbb{E}_{q(\boldsymbol{\alpha})}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{T}\right\rangle+\sum_{b=1}^{B}\left[\mathbf{y}_{ib}^{T}\boldsymbol{\Sigma}_{ib}^{-1}\mathbf{y}_{ib}+\log\det\boldsymbol{\Sigma}_{ib}\right].$

This theorem develops analytical formulae to update
$q(\boldsymbol{\theta}_{i},f_{i})$ for the $i$-th Mira. At each fixed $f_{i}$,
the canonical parameter for $q(\boldsymbol{\theta}_{i}|f_{i})$ is directly
available via (18). Notice the basis matrix $\mathbf{C}_{ib}$, the covariance
matrix $\boldsymbol{\Sigma}_{ib}$ as well as the a-priori mean vector
$\boldsymbol{\theta}_{0}$ all implicitly depend on the frequency $f_{i}$.
Besides, the optimal density $q(f_{i})$ is computed by (19). In our
implementation, $q(f_{i})$ is discretized and evaluated on a dense grid of
$[f_{\min},f_{\max}]$. The grid is gaped by a small value of $\Delta f$. After
that $q(f_{i})$ is normalized to become a proper density. The final algorithm
for updating
$q(\boldsymbol{\theta}_{i},f_{i})=q(\boldsymbol{\theta}_{i}|f_{i})q(f_{i})$ is
summarized in Algorithm 1. In our numerical studies, the computation is
performed on $500$ equally spaced points on the interval $[1/1000,1/100]$.

Algorithm 1 Updating Local Parameters

Input: $f_{\text{min}},f_{\text{max}},\Delta f$. The dataset
$\mathcal{D}_{i}$.

Output: $q(\boldsymbol{\theta}_{i},f_{i})$.

1:$f_{i}\leftarrow f_{\text{min}}$

2:while $f_{i}\leq f_{\text{max}}$ do

3: Update $q(\boldsymbol{\theta}_{i}|f_{i})$ by (18).

4: Compute $q(f_{i})$ by (19).

5: $f_{i}\leftarrow f_{i}+\Delta f$

6:end while

7:Normalize $q(f_{i})$ to become a proper density.

###  3.2 Update global parameters

Now, we address updating the global
$q(\boldsymbol{\alpha}),q(\boldsymbol{\Omega})$ and $q(\boldsymbol{\gamma})$,
with fixed local variational distributions
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})$ for all $i\in\mathcal{I}$. We
firstly develop the update of $q(\boldsymbol{\Omega})$ in detail. The way to
update $q(\boldsymbol{\alpha})$ and $q(\boldsymbol{\gamma})$ is derived
similarly.

When $q(\boldsymbol{\alpha})$, $q(\boldsymbol{\gamma})$ and all
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})$’s are fixed, optimizing the
ELBO with respect to $q(\boldsymbol{\Omega})$ is equivalent to maximizing
$\mathbb{E}_{q}\log\left[\frac{p(\boldsymbol{\Omega}|\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})}{q(\boldsymbol{\Omega})}\right]$.
This is because

$\displaystyle\underline{L}(\mathcal{D})=$
$\displaystyle\mathbb{E}_{q}\Big{\\{}\log\frac{p(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha},\mathcal{D})}{q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})}\Big{\\}}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{q}\log\left[\frac{p(\boldsymbol{\Omega}|\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})}{q(\boldsymbol{\Omega})}\right]+\mathbb{E}_{q}\log\left[\frac{p(\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})}{q(\mathbf{m},\boldsymbol{\beta},\mathbf{f},\boldsymbol{\Omega},\boldsymbol{\gamma},\boldsymbol{\alpha})/q(\boldsymbol{\Omega})}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{q}\log\left[\frac{p(\boldsymbol{\Omega}|\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})}{q(\boldsymbol{\Omega})}\right]+C^{\prime\prime}\,,$
(20)

where $C^{\prime\prime}$ is some constant invariant to the change of
$q(\boldsymbol{\Omega})$. In the above, the expectation
$\mathbb{E}_{q}(\cdot)$ is taken with respect to $q\in\mathcal{Q}$ in (15) at
the current iteration. From the joint distribution (13), we can find the full
conditional for $\boldsymbol{\Omega}$ as

$\displaystyle
p(\boldsymbol{\Omega}|\boldsymbol{\alpha},\boldsymbol{\gamma},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})\propto\exp\Big{\\{}-\frac{1}{2}\Big{\langle}\bar{n}\bar{\boldsymbol{\Omega}}^{-1}+\sum_{i=1}^{N}\boldsymbol{\beta}_{i}\boldsymbol{\beta}_{i}^{T},\boldsymbol{\Omega}\Big{\rangle}+\frac{N+\bar{n}-2B-1}{2}\log|\boldsymbol{\Omega}|\Big{\\}}\,.$

This full conditional is a Wishart distribution with canonical parameter
$\xi^{\boldsymbol{\Omega}1}=-(\bar{n}\bar{\boldsymbol{\Omega}}^{-1}+\sum_{i=1}^{N}\boldsymbol{\beta}_{i}\boldsymbol{\beta}_{i}^{T})/2$
and $\xi^{\boldsymbol{\Omega}2}=(N+\bar{n}-2B-1)/2$. In addition, note
$q(\boldsymbol{\Omega})\propto\exp\big{\\{}\langle\eta^{\boldsymbol{\Omega}1},\boldsymbol{\Omega}\rangle+\eta^{\boldsymbol{\Omega}2}\log|\boldsymbol{\Omega}|\big{\\}}$
is also Wishart with canonical parameters $\eta^{\boldsymbol{\Omega}1}$ and
$\eta^{\boldsymbol{\Omega}2}$. It can be verified the optimal
$q(\boldsymbol{\Omega})$ maximizing (20) has canonical parameters

$\eta^{\boldsymbol{\Omega}1}=\mathbb{E}_{q}(\xi^{\boldsymbol{\Omega}1})=-\frac{1}{2}\Big{[}\bar{n}\bar{\boldsymbol{\Omega}}^{-1}+\sum_{i=1}^{N}\mathbb{E}_{q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})}\left(\boldsymbol{\beta}_{i}\boldsymbol{\beta}_{i}^{T}\right)\Big{]}\,,$
(21)

and $\eta^{\boldsymbol{\Omega}2}=\xi^{\boldsymbol{\Omega}2}$. Although it is
straightforward to compute
$\eta^{\boldsymbol{\Omega}2}=\xi^{\boldsymbol{\Omega}2}=(N+\bar{n}-2B-1)/2$,
the involved computation in (21) is expensive for a large dataset. In
particular, the expectation with respect to
$q(\mathbf{m}_{i},\boldsymbol{\beta}_{i},f_{i})$ does not have a close form
expression. It can only be carried out by costly numerical integration.
Moreover, this has to be done for all of the $N$ observations during each
algorithm iteration, which is especially troublesome when $N$ is large.

Fortunately, the intensive computation can be avoided by SVI (Hoffman et al.,
2013), which only requires an unbiased estimation of natural gradient. For the
canonical parameter $\eta^{\boldsymbol{\Omega}1}$ of $q(\boldsymbol{\Omega})$,
its natural gradient of the ELBO is
$\mathbb{E}_{q}(\xi^{\boldsymbol{\Omega}1})-\eta^{\boldsymbol{\Omega}1}$. It
is the difference between the expected parameters of the full conditional and
the parameters of the current $q(\boldsymbol{\Omega})$. The gradient ascent
algorithm would update $q(\boldsymbol{\Omega})$ by
$\eta^{\boldsymbol{\Omega}1}\leftarrow\eta^{\boldsymbol{\Omega}1}+\kappa_{t}\left[\mathbb{E}_{q}(\xi^{\boldsymbol{\Omega}1})-\eta^{\boldsymbol{\Omega}1}\right]$
with some step size $\kappa_{t}$ at the $t$-th iteration. In the context of
stochastic gradient algorithm (Ghadimi and Lan, 2013), it is plausible to
employ an unbiased estimation of the natural gradient via sampling. For this
purpose, at the $t$-th iteration, we sample a subset of indices
$\widetilde{\mathcal{I}}_{t}=\\{i_{1},i_{2},\cdots,i_{I}\\}$ from
$\mathcal{I}=\\{1,\cdots,N\\}$ without replacement and with equal probability.
The cardinality of the set is $|\widetilde{\mathcal{I}}_{t}|=I$ for some fixed
$I$. This corresponds to the mini-batch strategy, and the batch size $I$ is
usually set to be $1,4,8,16$, etc. Given the set $\widetilde{\mathcal{I}}_{t}$
and for each $j\in\widetilde{\mathcal{I}}_{t}$, we further sample a frequency
value $\tilde{f}_{j}$ from the current variational distribution $q(f_{j})$.
Then, an unbiased estimate of the natural gradient
$\mathbb{E}_{q}(\xi^{\boldsymbol{\Omega}1})-\eta^{\boldsymbol{\Omega}1}$ is

$-\frac{1}{2}\Big{[}\bar{n}\bar{\boldsymbol{\Omega}}^{-1}+\frac{N}{I}\sum_{j\in\widetilde{\mathcal{I}}_{t}}\mathbb{E}_{q(\mathbf{m}_{j},\boldsymbol{\beta}_{j}|\tilde{f}_{j})}\Big{(}\boldsymbol{\beta}_{j}\boldsymbol{\beta}_{j}^{T}\Big{)}\Big{]}-\eta^{\boldsymbol{\Omega}1}.$

Based on the unbiased natural gradient estimate, the update for
$q(\boldsymbol{\Omega})$ becomes

$\displaystyle\eta^{\boldsymbol{\Omega}1}$
$\displaystyle\leftarrow\eta^{\boldsymbol{\Omega}1}+\kappa_{t}\Big{(}-\frac{1}{2}\Big{[}\bar{n}\bar{\boldsymbol{\Omega}}^{-1}+\frac{N}{I}\sum_{j\in\widetilde{\mathcal{I}}_{t}}\mathbb{E}_{q(\mathbf{m}_{j},\boldsymbol{\beta}_{j}|\tilde{f}_{j})}\Big{(}\boldsymbol{\beta}_{j}\boldsymbol{\beta}_{j}^{T}\Big{)}\Big{]}-\eta^{\boldsymbol{\Omega}1}\Big{)},$
(22) $\displaystyle\eta^{\boldsymbol{\Omega}2}$
$\displaystyle\leftarrow(N+\bar{n}-2B-1)/2.$ (23)

Notice we only need to compute with $I$ samples in the udpate (22). In
addition, the expectation with respect to
$q(\mathbf{m}_{j},\boldsymbol{\beta}_{j}|\tilde{f}_{j})$ is straightforward
because it is a Gaussian distribution. The optimal value for
$\eta^{\boldsymbol{\Omega}2}$ is directly provided in (23), and
$\eta^{\boldsymbol{\Omega}2}$ keeps fixed at this value during the iteration.
These features make the algorithm scale easily to a large collection of
observations.

Under the same principle, we can find the full conditional for $\gamma_{b}$ as
$q(\gamma_{b}|\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})\propto\exp\big{\\{}\xi_{b}^{\gamma
1}\gamma_{b}+\xi_{b}^{\gamma 2}\log\gamma_{b}\big{\\}}$, which is gamma
distribution with canonical parameter $\xi_{b}^{\gamma
1}=-\bar{r}-\sum_{i=1}^{N}(m_{ib}-\mathbf{d}(f_{i})^{T}\boldsymbol{\alpha}_{b})^{2}/2$
and $\xi_{b}^{\gamma 2}=N/2+\bar{\gamma}_{b}\bar{r}-1$. Note
$\mathbf{d}(f)=(1,\log_{10}(f),\log^{2}_{10}(f))^{T}$ is the PLR basis. With
the sampled subset of data and the sampled frequencies, the SVI algorithm
updates $q(\gamma_{b})$ by

$\displaystyle\xi_{b}^{\gamma 1}$ $\displaystyle\leftarrow\xi_{b}^{\gamma
1}+\kappa_{t}\Big{(}-\bar{r}-\frac{N}{I}\sum_{j\in\widetilde{\mathcal{I}}_{t}}\mathbb{E}_{q(\mathbf{m}_{j},\boldsymbol{\beta}_{j}|\tilde{f}_{j})}\mathbb{E}_{q(\boldsymbol{\alpha}_{b})}\left[m_{jb}-\mathbf{d}(\tilde{f}_{j})^{T}\boldsymbol{\alpha}_{b}\right]^{2}/2\Big{)},$
(24) $\displaystyle\xi_{b}^{\gamma 2}$ $\displaystyle\leftarrow
N/2+\bar{\gamma}_{b}\bar{r}-1.$ (25)

Similarly, the full conditional for $\boldsymbol{\alpha}_{b}$ is
$q(\boldsymbol{\alpha}_{b}|\boldsymbol{\gamma},\boldsymbol{\Omega},\boldsymbol{\beta},\mathbf{m},\mathbf{f},\mathcal{D})\propto\exp\Big{\\{}\left\langle\eta^{\alpha
1}_{b},\
\boldsymbol{\alpha}_{b}\boldsymbol{\alpha}_{b}\right\rangle+\left\langle\eta^{\alpha
2}_{b},\ \boldsymbol{\alpha}_{b}\right\rangle\Big{\\}}$. This full conditional
is multivariate Gaussian with canonical parameter $\eta^{\alpha
1}_{b}=-\frac{\bar{\delta}}{2}\mathbf{I}-\frac{\gamma_{b}}{2}\sum_{i=1}^{N}\mathbf{d}(f_{i})\mathbf{d}^{T}(f_{i})$,
and $\eta^{\alpha
2}_{b}=\bar{\delta}\bar{\alpha}_{b}+\gamma_{b}\sum_{i=1}^{N}m_{ib}\mathbf{d}(f_{i})$.
The SVI algorithm updates $q(\boldsymbol{\alpha}_{b})$ by

$\displaystyle\eta^{\alpha 1}_{b}$ $\displaystyle\leftarrow\eta^{\alpha
1}_{b}+\kappa_{t}\Big{(}-\frac{\bar{\delta}}{2}\mathbf{I}-\frac{\gamma_{b}}{2}\times\frac{N}{I}\sum_{j\in\widetilde{\mathcal{I}}_{t}}\mathbf{d}(\tilde{f}_{j})\mathbf{d}(\tilde{f}_{j})^{T}\Big{)},$
(26) $\displaystyle\eta^{\alpha 2}_{b}$ $\displaystyle\leftarrow\eta^{\alpha
2}_{b}+\kappa_{t}\Big{(}\bar{\delta}\bar{\alpha}_{b}+\gamma_{b}\times\frac{N}{I}\sum_{j\in\widetilde{\mathcal{I}}_{t}}\mathbb{E}_{q(\mathbf{m}_{j},\boldsymbol{\beta}_{j}|\tilde{f}_{j})}(m_{jb})\times\mathbf{d}(\tilde{f}_{j})\Big{)}.$
(27)

Algorithm 2 Main Algorithm

1:for $t=1,2,\cdots$ do

2: Set step size $\kappa_{t}=(c_{1}+t)^{-c_{2}}$.

3: Sample $\widetilde{\mathcal{I}}_{t}=\\{i_{1},i_{2},\cdots,i_{I}\\}$ from
$\mathcal{I}=\\{1,2,\cdots,N\\}$ without replacement.

4: for $j\in\widetilde{\mathcal{I}}_{t}$ do

5: Update local $q(\boldsymbol{\theta}_{j},f_{j})$ by Algorithm 1.

6: Sample $\tilde{f}_{j}$ according to $q(f_{j})$.

7: end for

8: Update $q(\boldsymbol{\Omega})$ by (22) and (23).

9: for $b=1,2,\cdots,B$ do

10: Update $q(\gamma_{b})$ by (24) and (25).

11: Update $q(\boldsymbol{\alpha}_{b})$ by (26) and (27).

12: end for

13:end for

The final algorithm is summarized in Algorithm 2. The algorithm chooses the
step size $\kappa_{t}=(c_{1}+t)^{-c_{2}}$. At the start of each iteration, a
small subset of samples is selected and their local parameters get updated in
Line 3–7 of Algorithm 2. With the sampled and updated local variational
distribution, the global parameters get updated by (22)–(27). In our numerical
studies, we set $c_{1}\in[1000,2000]$ and $c_{2}\in[0.5,1]$. We also take
mini-batch size $I=8$, and find that 1000 iterations are enough for
convergence.

##  4 Simulation

We have conducted two simulation experiments to compare the performance of our
method (SVI) with some existing methods. These include the generalized Lomb-
Scargle method (GLS, Zechmeister and Kürster, 2009) and its multi-band version
(MGLS, Vanderplas and Ivezic, 2015), the single-band semi-parametric model
(SP, He et al., 2016) and its multi-band extension (MSP, Yuan et al., 2018).

Table 2: Performance for Simulation I | GLS | MGLS | SP | MSP | SVI
---|---|---|---|---|---
ADE ($\times 10^{-4}$) | 5.98 | 3.97 | 5.70 | 2.18 | 0.52
RR (%) | 72.00 | 82.44 | 78.58 | 91.26 | 98.04

Table 3: SVI Confidence Set | Coverage
---|---
Nomimal (%) | 90 | 95 | 99 | 99.5
Actual (%) | 78.5 | 86.2 | 95.3 | 97.8

The first simulation uses the same simulated dataset used in Yuan et al.
(2018), where the light curves for 5,000 O-rich Miras were simulated in the
$I$, $J$, $H$, and $K_{s}$ bands. The cadence and sampling quality of the
light curves closely matched their actual survey of M33. Given the observation
time points, the signal light curves were generated using the Mira template of
Yuan et al. (2017). This template is trained from the well-observed OGLE LMC
data set. After this, the magnitude of the generated light curves were shifted
by $+6.27$ mag to account for the relative distance between the LMC and M33.
Finally, the light curve magnitude values were contaminated by realistic
noise. Four simulated multi-band Mira light curves are shown in the left
column of Figure 6.

Figure 6: In the left column are four simulated Mira light curves from
Simulation I. In the right column are the estimated variational distributions
$q(f_{i})$ for the corresponding Mira in the same row.

The five methods GLS, MGSL, SP, MSP, and SVI are applied to this dataset. The
single-band models (GLS and SP) are applied to the $I$-band dataset as it has
the most data points. For our proposed method, we use maximum a posterior
(MAP) as a point estimate for the frequency. In other words, for the $i$-th
star, its estimated frequency is $\hat{f}_{i}=\arg\max_{f_{i}}q(f_{i})$ with
the final estimated variational distribution $q(f_{i})$.

The accuracy of period estimation is measured by two metrics: recovery rate
and absolute deviation error. The recovery rate (RR) is the percentage of
light curves with accurately estimated frequency. Suppose for the $i$-th
simulated Mira, $f_{i0}$ is its true frequency and $\hat{f}_{i}$ is the
estimated frequency. The estimation is considered to be accurate if their
absolute difference is less than a threshold value, i.e., if
$|f_{i0}-\hat{f}_{i}|\leq\lambda$. We choose $\lambda=2.7\times 10^{-4}$ in
this work, which is the average half distance to the sidelobes in the
frequency spectra (He et al., 2016). Then, the recovery rate (RR) is computed
as
$\mathrm{RR}=\frac{1}{N}\sum_{i=1}^{N}\mathrm{I}(|f_{i0}-\hat{f}_{i}|\leq\lambda)$,
where $\mathrm{I}(\cdot)$ is the indicator function. The other employed metric
is the absolute deviation error (ADE), computed by
$\mathrm{ADE}=\frac{1}{N}\sum_{i=1}^{N}|f_{i}-\hat{f}_{i}|$.

For all methods, their recovery rate and ADE on this dataset are reported in
Table 3. Notice the single-band GLS method has worst performance in both
metrics. By exploiting more observational information, the multi-band methods
(MGLS and MSP) increase the RR by 10% and reduces ADE by more than 40%,
compared with their single-band counterparts (GLS and SP). It is also notable
that SP (or MSP) outperforms GLS (or MGLS resp.). Recall the SP and MSP have
modeled the sinusoidal residuals by Gaussian process. The improvement
demonstrates the effectiveness of modeling the quasi-periodic signal of Mira
light curves. Lastly, we note our SVI method has the best performance. The
proposed SVI has remarkably estimated almost all samples accurately. It has a
recovery rate of $98.20\%$ and ADE accuracy of $0.52\times 10^{-4}$.

In the right column of Figure 6, we plot the fitted variational distribution
$q(f_{i})$ for each simulated Mira light curve shown in the left column of the
same row. The red reversed triangles mark the estimated frequencies
$\hat{f}_{i}$, and the blue triangles mark the true frequencies. The first two
rows demonstrate the cases when the true frequencies are correctly identified.
The last two rows show the cases when the MAP $\hat{f}_{i}$ fails to estimate
the true frequency. However, in these cases when the signal is weak, our
proposed method still recovers the multi-modal posterior of the frequency, and
the true frequency lies around one of the posterior modes.

The plots of $q(f_{i})$ illustrate the power of this method to quantify the
frequency estimation uncertainty. To further assess the uncertainty
quantification from $q(f_{i})$, we construct a $1-\alpha$ level confidence set
for each Mira in the simulated dataset. The confidence set
$\mathcal{C}_{i}=\\{f_{i}:q(f_{i})>c_{i}\\}$ is constructed with some $c_{i}$
such that $\int_{\mathcal{C}_{i}}q(f_{i})\,\mathrm{d}f_{i}=1-\alpha$. Note
each $\mathcal{C}_{i}$ may be composed of a few short intervals due to the
multimodal nature of $q(f_{i})$. For different value of $\alpha$, the
confidence set is computed, and the actual coverage rate is evaluated over the
simulated dataset. The empirical results are reported in Table 3. For example,
when the nominal coverage rate is $95\%$ (with $\alpha=5\%$), the actual
coverage rate is only $86.2\%$. It indicates the confidence sets are not wide
enough and have an under-coverage issue. This problem is inherent for mean
filed variational inference which underestimates estimation uncertainty (Wang
and Blei, 2019).

Taking a closer look of the coverage rate issue, we found that, in most cases
when the confidence set misses the true frequency, the true frequency resides
very close to the confidence set; it is rare that the confidence set
completely misses the true frequency. More precisely, let $\mathcal{C}_{i}$ be
the 99% confidence set for the $i$ sample. Let $\lambda=2.7\times 10^{-4}$ be
the threshold to determine the recovery rate (RR) as used in the paper. We say
that the set $\mathcal{C}_{i}$ “entirely misses” the true frequency, if
$\text{dist}(f_{0i},\mathcal{C}_{i})>\lambda$; and it “almost covers” the true
frequency, if otherwise. In other words, the $\mathcal{C}_{i}$ is considered
as completely failed if the distance of the true frequency $f_{i0}$ to the set
$\mathcal{C}_{i}$ is larger than $\lambda$. Out of the total 5,000 simulated
samples, for 234 samples (4.68%) the 99% confidence set fails to cover the
corresponding the true frequencies; among them, for only 19 (0.38%) samples,
the confidence set entirely misses the true frequency, and for 215 samples
(4.3%), the confidence set almost covers the true frequency.

To further assess the performance of the proposed method, we conducted
another, more comprehensive simulation study where 90 datasets with distinct
characteristics were generated by considering 90 different combinations of
sampling cadence patterns, different observation noise levels, and various
number of sampling points. Some of these datasets are more challenging than
others with higher level of data sparsity and noise. Our proposed method still
delivers the best performance in terms of the two metrics (RR and ADE) for
each dataset. Details of this simulation study are presented in Section S.4 of
the Supplementary Material (He et al., 2020).

(a)

(b)

(c)

Figure 7:  The horizontal axis is the logarithm of SVI period, and the
vertical axis is the difference of average magnitude between the results of
SVI and Yuan et al. (2018). Left, center and right panels show $J-$, $H$ and
$K_{s}$ magnitudes, respectively. The horizontal dashed lines indicate the
average difference between the two methods.

##  5 The M33 Dataset

In this section we analyze the dataset of M33 Miras from Yuan et al. (2018).
The authors collected images and photometric measurements from several sources
(Javadi et al., 2015; Macri et al., 2001; Pellerin and Macri, 2011) to
generate light curves in four bands ($I$, $J$, $H$, and $K_{s}$) with a median
number of observations of 68, 5, 6, and 11, respectively. Yuan et al. (2018)
identified 1,781 candidate Miras and classified 1,265 of them as O-rich, 88 as
C-rich, and 428 as unknown, respectively. The O-rich Miras were used to
estimate distance moduli (logarithms of the distance) for M33. From the $J$-,
$H$\- and $K_{s}$-band PLRs, they obtained values of $24.82\pm 0.06$ mag,
$24.82\pm 0.06$ mag, and $24.75\pm 0.06$ mag, respectively. Earlier works
based on Cepheid variables in the same galaxy yielded estimates of its
distance modulus as $24.65\pm 0.12$ mag (Macri, 2001), $24.53\pm 0.11$ mag
(Scowcroft et al., 2009), and $24.62\pm 0.07$ mag (Gieren et al., 2013).
Detached eclipsing binaries resulted in an estimate of $24.92\pm 0.12$ mag
(Bonanos et al., 2006), while RR Lyrae variables provided $24.67\pm 0.08$ mag
(Sarajedini et al., 2006).

Figure 8:  Comparison of the period estimated by SVI and by Yuan et al.
(2018). The red points are Miras with inconsistent periods from the two
methods.

Figure 9: Scatterplot for the sinusoidal coefficients for M33 Miras.

We carry out an analysis with our proposed method using the 1,265 candidate
Miras classified as O-rich by Yuan et al. (2018). Our method fits the light
curves reasonably well, as can be visually assessed by plotting the MAP
estimator; see Section S.5 of the Supplementary Material for some randomly
selected examples. We obtain the frequency and average magnitude for each
Mira. The frequency of the $i$-th Mira is directly estimated by
$\hat{f}_{i}=\arg\max_{f_{i}}q(f_{i})$. The uncertainty for its estimated
period $\hat{p}_{i}=1/\hat{f}_{i}$ is computed by
$\sigma^{2}(\hat{p}_{i})=\int(1/f_{i})^{2}q(f_{i})\mathrm{d}f_{i}-\big{[}\int(1/f_{i})q(f_{i})\mathrm{d}f_{i}\big{]}^{2}$.
The computation of $\hat{f}_{i}$ is straightforward and more reliable, as our
model has built-in PLR to suppress aliasing frequencies. In contrast, the work
of Yuan et al. (2018) involves an ad-hoc post-estimation correction. For each
Mira, they test both the primary and secondary peaks in their log-likelihood.
Their final estimated frequency is selected as the one that best fits the
overall PLR.

Table 4: $J$-, $H$\- & $K_{s}$-band PLR coefficients for candidate M33 O-rich
Miras

galaxy | band |  linear (period$<400$d) | quadratic
---|---|---|---
|  | $a_{0}$ | $a_{1}$ | $\sigma$ | $a_{0}$ | $a_{1}$ | $a_{2}$ | $\sigma$
LMC (Yuan+17) | $J$ | $12.67\pm 0.01$ | $-3.48\pm 0.09$ | $0.15$ | $12.70\pm 0.01$ | $-3.49\pm 0.09$ | $-1.54\pm 0.23$ | $0.15$
M33 (Yuan+18) | $J$ | $18.92\pm 0.01$ | $\dots$ | $0.25$ | $18.94\pm 0.01$ | $\dots$ | $\dots$ | $0.25$
M33 (SVI) | $J$ | $18.97\pm 0.01$ | $-3.32\pm 0.06$ | $0.22$ | $19.01\pm 0.01$ | $-3.36\pm 0.05$ | $-1.52\pm 0.02$ | $0.23$
LMC (Yuan+17) | $H$ | $11.91\pm 0.01$ | $-3.64\pm 0.09$ | $0.16$ | $11.96\pm 0.01$ | $-3.59\pm 0.10$ | $-3.40\pm 0.31$ | $0.16$
M33 (Yuan+18) | $H$ | $18.17\pm 0.01$ | $\dots$ | $0.24$ | $18.27\pm 0.01$ | $\dots$ | $\dots$ | $0.26$
M33 (SVI) | $H$ | $18.24\pm 0.01$ | $-3.59\pm 0.06$ | $0.21$ | $18.27\pm 0.01$ | $-3.11\pm 0.04$ | $-3.24\pm 0.02$ | $0.22$
LMC (Yuan+17) | $K_{s}$ | $11.56\pm 0.01$ | $-3.77\pm 0.07$ | $0.12$ | $11.59\pm 0.01$ | $-3.77\pm 0.08$ | $-2.23\pm 0.20$ | $0.12$
M33 (Yuan+18) | $K_{s}$ | $17.78\pm 0.01$ | $\dots$ | $0.21$ | $17.83\pm 0.01$ | $\dots$ | $\dots$ | $0.22$
M33 (SVI) | $K_{s}$ | $17.85\pm 0.01$ | $-3.93\pm 0.06$ | $0.20$ | $17.87\pm 0.01$ | $-3.68\pm 0.04$ | $-2.20\pm 0.02$ | $0.20$

A comparison of the periods estimated by our method and the MSP method is
shown in Figure 9. Most of the estimated periods from the two methods are
consistent, with only $5.8\%$ Miras (red points) deviating by more than 10%.
We also compare the average $JHK_{s}$ magnitude between our determinations and
that of Yuan et al. (2018) in Figure 7(a), 7(b) and 7(c). Our determination
has slighter higher average magnitude values. Nevertheless, the two results
are generally consistent, despite there are some outliers. The standard
deviations of the difference are 0.16, 0.13, and 0.10 mag for the $J$, $H$,
and $K_{s}$ bands, respectively.

Table 5: Derived Distance Moduli for M33

band | $\Delta a_{0}$ | $\Delta\bar{m}$ | $\Delta A_{\lambda}$ | $\Delta\mathrm{ct}$ | $\Delta\mu$ | $\mu_{\mathrm{LMC}}$ | $\mu_{\mathrm{M33}}$
---|---|---|---|---|---|---|---
$J$ | $6.311\pm 0.014$ | $-0.036\pm 0.035$ | $0.029\pm 0.008$ | $0.016\pm 0.036$ | $6.320\pm 0.053$ | $18.493\pm 0.048$ | $24.813\pm 0.071$
$H$ | $6.312\pm 0.014$ | $-0.033\pm 0.029$ | $0.018\pm 0.005$ | $0.010\pm 0.040$ | $6.307\pm 0.052$ | $18.493\pm 0.048$ | $24.800\pm 0.070$
$K_{s}$ | $6.288\pm 0.014$ | $-0.026\pm 0.024$ | $0.012\pm 0.003$ | $-0.007\pm 0.032$ | $6.267\pm 0.042$ | $18.493\pm 0.048$ | $24.760\pm 0.064$

In addition, we compute the posterior mean of the sinusoidal coefficient
$\boldsymbol{\beta}_{i}$ for each Mira. The scatter plot of the expected
$\boldsymbol{\beta}_{i}$’s (the expectation is taken with respect to
$q(\boldsymbol{\beta}_{i})$) is shown in Figure 9. Notice that the
coefficients for the $JHK_{s}$ bands have a tight correlation, and their
relation with the $I$-band coefficients $\beta_{I1},\beta_{I2}$ is also non-
negligible. The correlation structure in the real data justifies our model
Eqn. (4) on $\boldsymbol{\beta}_{i}$. By exploiting this structure, our method
is able to constrain the model degree of freedom and estimate parameters more
accurately.

We note the doughnut shape of the posterior joint distribution of
$\beta_{I1},\beta_{I2}$ shown in Figure 9 does not suggest model
misspecification. In fact, it is the Bayesian method in action. The posterior
distribution needs not be consistent with the prior. In Bayesian statistics,
it is expected that when there is enough data, the information provided in the
data will “wash out the prior.”

Next, we will compare our estimated PLRs with the results from Yuan et al.
(2017) and Yuan et al. (2018). For each band, Yuan et al. (2018) fit a
quadratic PLR

$M=a_{0}+a_{1}(\log_{10}(P)-2.3)+a_{2}(\log_{10}(P)-2.3)^{2}$ (28)

to all Miras; and a linear PLR $M=a_{0}+a_{1}(\log_{10}(P)-2.3)$ to Miras with
period less than 400 days. Notice the above quadratic PLR is simply an
equivalent reparameterization of our PLR in (11). In the work of Yuan et al.
(2018), the PLRs are firstly fitted to the LMC Miras (Yuan et al., 2017).
After that, they fix the coefficients $a_{1}$ and $a_{2}$, and only fit the
intercept $a_{0}$ to the M33 Miras. Table 4 presents their results. The
coefficients for the LMC Miras are in the rows of LMC (Yuan+17), and their
results on the M33 Miras are in the rows named M33 (Yuan+18).

The quadratic PLRs from our models are plotted as red curves in Figure 10. In
the figure, each point represents one Mira with our estimated average
magnitude and period. For Miras with period less than 400 days, we also fit
additional linear PLRs, which are shown as black lines. Our coefficients are
presented in the rows named M33 (SVI) in Table 4. For comparison, the reported
coefficient values have been converted to respect the equivalent
parameterization (28). Notice our estimated $a_{1}$ and $a_{2}$ are consistent
with the LMC results from Yuan et al. (2017) within the error budget.

In Table 4, the column $\sigma$ is the standard deviation of the PLR
residuals. Smaller $\sigma$ means tighter PLRs. It is remarkable that our
model achieves competitive performance as the M33 result of Yuan et al.
(2018), without their ad-hoc period correction tricks. This verifies the
quality the automatically estimated PLRs embedded in our model.

Figure 10: The PLRs of M33 Miras in the $J$ (left), $H$ (center) and $K_{s}$
(right) bands. The vertical axes shown the average magnitude in the respective
band while the horizontal axes are the logarithm of the period. Red curves are
the fitted PLRs from our SVI method. Black lines are fitted linear PLRs for
Miras with period less than 400 days.

Lastly, we compute the distance modulus of the M33 galaxy. Denote
$a_{0}^{\mathrm{LMC}}$ as the zero-order coefficient of the LMC quadratic PLR
from Yuan et al. (2017). Denote $a_{0}^{\mathrm{M33}}$ as the zero-order
coefficient of our M33 quadratic PLR . Then, up to a few corrections, $\Delta
a_{0}=a_{0}^{\mathrm{M33}}-a_{0}^{\mathrm{LMC}}$ is equal to the difference in
distance moduli of the two galaxies. These corrections include: conversion to
flux-averaged magnitude $\Delta\bar{m}$, correction for the different
interstellar extinction $\Delta A_{\lambda}$ towards LMC and M33; and the
color term bias $\Delta\mathrm{ct}$ of the photometric calibration. We note
the flux correction $\Delta\bar{m}$ is necessary because the average magnitude
was computed in units of flux for the LMC Miras. This correction is carried
out as follows: We fit the signal light curve $s_{ib}(t)$ in (10) for each
Mira over the whole observational time domain. The signal $s_{ib}(t)$ is in
the unit of magnitude and can be converted to flux via
$s^{\prime}_{ib}(t)=10^{-0.4\cdot s_{ib}(t)}$. The average value $l_{ib}$ of
$s^{\prime}_{ib}(t)$ is computed, and it is converted back to magnitude by
$m_{ib}^{\prime}=-2.5\log_{10}(l_{ib})$. Note some irrelevant constants have
been dropped during these conversions. Then, the correction value
$\Delta\bar{m}$ for each band is the average difference
$m_{ib}^{\prime}-m_{ib}$ from all Miras.

After making these corrections, we obtain the relative distance modulus
$\Delta\mu=\Delta a_{0}+\Delta\bar{m}+\Delta A_{\lambda}+\Delta\mathrm{ct}$
between the LMC and M33. We adopt a distance modulus
$\mu_{\mathrm{LMC}}=18.493\pm 0.048$ for LMC (Pietrzyński et al., 2013). Thus,
the distance modulus for M33 is
$\mu_{\mathrm{M33}}=\mu_{\mathrm{LMC}}+\Delta\mu$. This procedure is applied
separately to the $J$-, $H$\- and $K_{s}$-band data and summarized in Table 5.
From the respective PLRs, our analysis yields distance moduli of $24.813\pm
0.071$ mag, $24.800\pm 0.070$ mag, and $24.760\pm 0.064$ mag, respectively.
These are also consistent with earlier estimates.

##  6 Effect of downsampling

In Section 5, we showed that our proposed SVI method produces distance
estimates and uncertainties which are consistent with other methods without
resorting to ad hoc period estimate corrections. This section studies how
robust these estimates are to downsample the light curves, e.g., how do
distance estimates change as the number of photometric measurements and/or
light curve goes down.

We conducted an experiment of downsampling the actual M33 data set used in
Section 5. We compared our SVI method with two other multi-band methods, MGLS
and MSP, without using any ad hoc period estimate corrections. To assess the
reliability and quality of the estimated PLR for distance estimation in
downsampling, we measured the changes in the PLR intercept and in the PLR
scattering between the downsampled estimates and the full sample estimates.

In each replication of the experiment, $\tilde{N}=50$ or $\tilde{N}=100$ Miras
are randomly drawn from the original $N=1,265$ candidate O-rich Miras. For
each selected Mira, we also put a restriction $\bar{n}_{b}$ over the number of
photometric measurements for its $b$-th band light curve
($b\in\\{I,J,H,K_{s}\\}$). In particular, for the $b$-th band light curve of
the $i$-th selected Mira, we sample
$\tilde{n}_{ib}=\min\\{n_{ib},\bar{n}_{b}\\}$ photometric measurements without
replacement from the original $n_{ib}$ measurements. Six combinations of
$\tilde{N}$ and $\bar{n}_{b}$ are considered as in Table 6.

Table 6: Experiment settings for downsampling the actual M33 observations Setting | $\tilde{N}$ | $\bar{n}_{I}$ | $\bar{n}_{J}$ | $\bar{n}_{H}$ | $\bar{n}_{K_{s}}$
---|---|---|---|---|---
S1 | 50 | 10 | 5 | 5 | 5
S2 | 50 | 20 | 10 | 10 | 10
S3 | 50 | 30 | 15 | 15 | 15
S4 | 100 | 10 | 5 | 5 | 5
S5 | 100 | 20 | 10 | 10 | 10
S6 | 100 | 30 | 15 | 15 | 15

Each setting in Table 6 is repeated 20 times. In each replication of one
setting, three multi-band methods (MGLS, MSP and SVI) are applied to estimate
the period and the average magnitude for each down-sampled Mira. To fit a
quadratic PLR for MGLS and MSP, we fix the coefficient $a_{1}$ and $a_{2}$
(estimated from LMC data) and only estimate the intercept term $a_{0}$.
Suppose the intercept estimated from the subsample is denoted as
$\tilde{a}_{0}$ and the full sample estimate is denoted $\hat{a}_{0}$. The
first row of Figure 11 shows the boxplot of the difference
$\tilde{a}_{0}-\hat{a}_{0}$. The three panels in the first row correspond to
the $J$, $H$ and $K_{s}$ bands, respectively; and the results for MGLS, MSP
and SVI are colored as red, green, and blue, respectively. The intercept
difference $\tilde{a}_{0}-\hat{a}_{0}$ is closest to $0$ for our SVI method,
and the variability is also small. As comparison, MGLS and MSP have much
larger bias and variability than SVI. Compared with MGLS, the MSP method has
smaller bias in the $H$ band but larger bias in the $J$ and $K_{s}$ bands.

Figure 11: Method comparison via down-sampling the real M33 observation. The
first row shows the difference between the downsampled PLR intercept and the
full sample PLR intercept. The second row shows the standard deviation ratio
of the downsampled PLR residuals and the full sample PLR residuals. The three
columns correspond to the $J$ ,$H$ and $K_{s}$ bands, respectively. The result
of 20 replications for MGLS, MSP and SVI are shown in the red, green and blue
boxplots, respectively.

We also compare the scattering of the PLR residuals. Let $\tilde{\sigma}$ be
the standard deviation of the PLR residuals for each method in each experiment
replication. Also, let $\hat{\sigma}$ be the standard deviation of the PLR
residuals obtained by our SVI method for the full sample. The second row of
Figure 11 shows the boxplot of the ratio $\tilde{\sigma}/\hat{\sigma}$ for the
20 replications of the three methods. The residual scattering ratio is only
slightly above $1$ for our SVI method over the six scenarios of downsampling.
As comparison, MGLS and MSP have substantially larger residual scattering
ratio. The residual scattering of MSP is smaller than that of MGLS.

These results suggest that our SVI is much more robust to downsample Mira
light curves than MGLS and MSP methods. It would be an interesting research
topic to study the implication of these results for telescope cadence
selection and the necessary density of time domain sampling.

## Acknowledgements

He’s Project 11801561 was supported by NSFC. L. Macri acknowledges support
from the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University. Huang’s work was partially supported by NSF grants DMS-1208952,
IIS-1900990, and CCF-1956219. The authors would like to thank two anonymous
reviewers for helpful comments which led to inclusion of Section 6 and other
significant improvements of the paper.

Supplementary proof and results. The supplementary material contains a review
of the notations of the exponential family used in the paper. It also contains
the proof of Theorem 3.1 and additional numerical results.

Supplementary code and dataset. The code and datasets are released online at
https://github.com/shiyuanhe/sviperiodplr. They serve to reproduce the
simulation and real data analysis results in this work.

## References

  * Blei et al. (2017) Blei, D. M., A. Kucukelbir, and J. D. McAuliffe (2017).  Variational inference: A review for statisticians.  Journal of the American Statistical Association 112(518), 859–877.
  * Bonanos et al. (2006) Bonanos, A. Z., K. Z. Stanek, R. P. Kudritzki, L. M. Macri, D. D. Sasselov, J. Kaluzny, P. B. Stetson, D. Bersier, F. Bresolin, T. Matheson, et al. (2006).  The first direct distance determination to a detached eclipsing binary in m33.  The Astrophysical Journal 652(1), 313.
  * Flewelling (2018) Flewelling, H. (2018, January).  Pan-STARRS Data Release 2.  In American Astronomical Society Meeting Abstracts #231, Volume 231 of American Astronomical Society Meeting Abstracts, pp. 436.01.
  * Freedman et al. (2001) Freedman, W. L., B. F. Madore, B. K. Gibson, L. Ferrarese, D. D. Kelson, S. Sakai, J. R. Mould, R. C. Kennicutt Jr, H. C. Ford, J. A. Graham, et al. (2001).  Final results from the hubble space telescope key project to measure the hubble constant.  The Astrophysical Journal 553(1), 47.
  * Gelman (2006) Gelman, A. (2006).  Prior distributions for variance parameters in hierarchical models (comment on article by browne and draper).  Bayesian analysis 1(3), 515–534.
  * Ghadimi and Lan (2013) Ghadimi, S. and G. Lan (2013).  Stochastic first-and zeroth-order methods for nonconvex stochastic programming.  SIAM Journal on Optimization 23(4), 2341–2368.
  * Gieren et al. (2013) Gieren, W., M. Górski, G. Pietrzyński, P. Konorski, K. Suchomska, D. Graczyk, B. Pilecki, F. Bresolin, R.-P. Kudritzki, J. Storm, et al. (2013).  The araucaria project. a distance determination to the local group spiral m33 from near-infrared photometry of cepheid variables.  The Astrophysical Journal 773(1), 69.
  * He et al. (2020) He, S., Z. Lin, W. Yuan, L. M. Macri, and J. Z. Huang (2020).  Supplement to “simultaneous inference of periods and period-luminosity relations for mira variable stars”.
  * He et al. (2016) He, S., W. Yuan, J. Z. Huang, J. Long, and L. M. Macri (2016).  Period estimation for sparsely sampled quasi-periodic light curves applied to miras.  The Astronomical Journal 152(6), 164.
  * Hoffman et al. (2013) Hoffman, M. D., D. M. Blei, C. Wang, and J. Paisley (2013).  Stochastic variational inference.  The Journal of Machine Learning Research 14(1), 1303–1347.
  * Huang et al. (2018) Huang, C. D., A. G. Riess, S. L. Hoffmann, C. Klein, J. Bloom, W. Yuan, L. M. Macri, D. O. Jones, P. A. Whitelock, S. Casertano, et al. (2018).  A near-infrared period–luminosity relation for miras in ngc 4258, an anchor for a new distance ladder.  The Astrophysical Journal 857(1), 67.
  * Javadi et al. (2015) Javadi, A., M. Saberi, J. T. van Loon, H. Khosroshahi, N. Golabatooni, and M. T. Mirtorabi (2015).  The uk infrared telescope m33 monitoring project–iv. variable red giant stars across the galactic disc.  Monthly Notices of the Royal Astronomical Society 447(4), 3973–3991.
  * Jordan et al. (1999) Jordan, M. I., Z. Ghahramani, T. S. Jaakkola, and L. K. Saul (1999, Nov).  An introduction to variational methods for graphical models.  Machine Learning 37(2), 183–233.
  * Lomb (1976) Lomb, N. R. (1976).  Least-squares frequency analysis of unequally spaced data.  Astrophysics and space science 39(2), 447–462.
  * Long et al. (2016) Long, J. P., E. C. Chi, and R. G. Baraniuk (2016).  Estimating a common period for a set of irregularly sampled functions with applications to periodic variable star data.  The Annals of Applied Statistics 10(1), 165–197.
  * Macri et al. (2001) Macri, L., K. Stanek, D. Sasselov, M. Krockenberger, and J. Kaluzny (2001).  Direct distances to nearby galaxies using detached eclipsing binaries and cepheids. vi. variables in the central part of m33.  The Astronomical Journal 121(2), 870.
  * Macri (2001) Macri, L. M. M. (2001).  The extragalactic distance scale.  Ph. D. thesis, Harvard University.
  * Mowlavi et al. (2018) Mowlavi, N., I. Lecoeur-Taïbi, T. Lebzelter, L. Rimoldini, D. Lorenz, M. Audard, J. De Ridder, L. Eyer, L. Guy, B. Holl, et al. (2018).  Gaia data release 2-the first gaia catalogue of long-period variable candidates.  Astronomy & Astrophysics 618, A58.
  * Pellerin and Macri (2011) Pellerin, A. and L. M. Macri (2011).  The m 33 synoptic stellar survey. i. cepheid variables.  The Astrophysical Journal Supplement Series 193(2), 26\.
  * Pietrzyński et al. (2013) Pietrzyński, G., W. Gieren, D. Graczyk, I. Thompson, B. Pilecki, N. Nardetto, R.-P. Kudritzki, F. Bresolin, G. Bono, P. P. Moroni, P. Konorski, M. Gorski, J. Storm, R. Smolec, and P. Karczmarek (2013, February).  A precise and accurate distance to the large magellanic cloud from late-type eclipsing-binary systems.  In R. de Grijs (Ed.), Advancing the Physics of Cosmic Distances, Volume 289 of IAU Symposium, pp. 169–172.
  * Rasmussen and Williams (2005) Rasmussen, C. E. and C. K. I. Williams (2005).  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning).  The MIT Press.
  * Riess et al. (2019) Riess, A. G., S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic (2019).  Large magellanic cloud cepheid standards provide a 1% foundation for the determination of the hubble constant and stronger evidence for physics beyond $\lambda$cdm.  The Astrophysical Journal 876(1), 85.
  * Riess et al. (2011) Riess, A. G., L. Macri, S. Casertano, H. Lampeitl, H. C. Ferguson, A. V. Filippenko, S. W. Jha, W. Li, and R. Chornock (2011).  A 3% solution: determination of the hubble constant with the hubble space telescope and wide field camera 3.  The Astrophysical Journal 730(2), 119.
  * Riess et al. (2016) Riess, A. G., L. M. Macri, S. L. Hoffmann, D. Scolnic, S. Casertano, A. V. Filippenko, B. E. Tucker, M. J. Reid, D. O. Jones, J. M. Silverman, R. Chornock, P. Challis, W. Yuan, P. J. Brown, and R. J. Foley (2016, July).  A 2.4% Determination of the Local Value of the Hubble Constant.  ApJ 826, 56.
  * Sarajedini et al. (2006) Sarajedini, A., M. Barker, D. Geisler, P. Harding, and R. Schommer (2006).  Rr lyrae variables in m33. i. evidence for a field halo population.  The Astronomical Journal 132(3), 1361.
  * Scargle (1982) Scargle, J. D. (1982, December).  Studies in astronomical time series analysis. II - Statistical aspects of spectral analysis of unevenly spaced data.  The Astrophysical Journal 263, 835–853.
  * Scowcroft et al. (2009) Scowcroft, V., D. Bersier, J. Mould, and P. R. Wood (2009).  The effect of metallicity on cepheid magnitudes and the distance to m33.  Monthly Notices of the Royal Astronomical Society 396(3), 1287–1296.
  * Soszyński et al. (2009) Soszyński, I., A. Udalski, M. K. Szymański, M. Kubiak, G. Pietrzyński, Ł. Wyrzykowski, O. Szewczyk, K. Ulaczyk, and R. Poleski (2009, March).  The Optical Gravitational Lensing Experiment. The OGLE-III Catalog of Variable Stars. III. RR Lyrae Stars in the Large Magellanic Cloud.  ACTA ASTRONOMICA 59, 1–18.
  * Udalski et al. (2008) Udalski, A., M. K. Szymański, I. Soszyński, and R. Poleski (2008).  The optical gravitational lensing experiment. final reductions of the ogle-iii data.  ACTA ASTRONOMICA 58, 69–87.
  * Vanderplas and Ivezic (2015) Vanderplas, J. T. and Ž. Ivezic (2015).  Periodograms for multiband astronomical time series.  The Astrophysical Journal 812(1), 18.
  * Wainwright et al. (2008) Wainwright, M. J., M. I. Jordan, et al. (2008).  Graphical models, exponential families, and variational inference.  Foundations and Trends® in Machine Learning 1(1–2), 1–305.
  * Wang and Blei (2019) Wang, Y. and D. M. Blei (2019).  Frequentist consistency of variational bayes.  Journal of the American Statistical Association 114(527), 1147–1161.
  * Whitelock et al. (2014) Whitelock, P., J. Menzies, M. Feast, F. Nsengiyumva, and N. Matsunaga (2014).  Vizier online data catalog: Jhks photometry of agb stars in ngc 6822 (whitelock+, 2013).  VizieR Online Data Catalog 742.
  * Whitelock and Feast (2014) Whitelock, P. A. and M. W. Feast (2014, July).  Gaia, Variable Stars and the Distance Scale.  In EAS Publications Series, Volume 67 of EAS Publications Series, pp. 263–269.
  * Yuan et al. (2017) Yuan, W., L. M. Macri, S. He, J. Z. Huang, S. M. Kanbur, and C.-C. Ngeow (2017).  Large magellanic cloud near-infrared synoptic survey. v. period–luminosity relations of miras.  The Astronomical Journal 154(4), 149.
  * Yuan et al. (2018) Yuan, W., L. M. Macri, A. Javadi, Z. Lin, and J. Z. Huang (2018).  Near-infrared mira period–luminosity relations in m33.  The Astronomical Journal 156(3), 112.
  * Zechmeister and Kürster (2009) Zechmeister, M. and M. Kürster (2009).  The generalised lomb-scargle periodogram. a new formalism for the floating-mean and keplerian periodograms.  Astronomy and Astrophysics 496(2), 577–584.
  * Zhang et al. (2019) Zhang, C., J. Butepage, H. Kjellstrom, and S. Mandt (2019).  Advances in variational inference.  IEEE Transactions on Pattern Analysis and Machine Intelligence 41(8), 2008–2026.

Supplementary Materials to “Simultaneous inference of periods and period-
luminosity relations for Mira variable stars”

##  S.1 Canonical Exponential Family

This section reviews the notation for the exponential family used in the
paper. The general form for the canonical exponential family is

$f(\mathbf{x}|\eta)=h(\mathbf{x})\exp\\{\eta^{T}t(\mathbf{x})-A(\eta)\\}.$

where $\eta$ is the canonical parameter, $r(\mathbf{x})$ is the sufficient
statistics, and $A(\eta)$ is the log-partition function. It is known that
$\mathbb{E}(t(\mathbf{x}))=\nabla_{\eta}A(\eta)$.

For a $p$ dimensional multivariate Gaussian distribution with mean
$\boldsymbol{\mu}$ and precision matrix $\boldsymbol{\Omega}$, its density is

$\displaystyle\exp\Big{\\{}-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Omega}(\mathbf{x}-\boldsymbol{\mu})+\frac{1}{2}\log\det\boldsymbol{\Omega}+\frac{1}{2}\log
2\pi\Big{\\}}$ $\displaystyle=$
$\displaystyle\exp\Big{\\{}\langle-\frac{1}{2}\boldsymbol{\Omega},\mathbf{x}\mathbf{x}^{T}\rangle+\langle\boldsymbol{\Omega}\boldsymbol{\mu},\mathbf{x}\rangle-\frac{1}{2}\boldsymbol{\mu}^{T}\boldsymbol{\Omega}\boldsymbol{\mu}+\frac{1}{2}\log\det\boldsymbol{\Omega}+\frac{1}{2}\log
2\pi\Big{\\}}$

with $t(\mathbf{x})=(\mathbf{x}\mathbf{x}^{T},\mathbf{x})$,
$\eta=(\eta^{1},\eta^{2})=(-\frac{1}{2}\boldsymbol{\Omega},\boldsymbol{\Omega}\boldsymbol{\mu})$
and
$A(\eta)=\frac{1}{2}\boldsymbol{\mu}^{T}\boldsymbol{\Omega}\boldsymbol{\mu}-\frac{1}{2}\log\det\boldsymbol{\Omega}$.
In addition, it holds that

$\displaystyle\mathbb{E}(\mathbf{x})=-\frac{1}{2}(\eta^{1})^{-1}\eta^{2},\;\mathrm{Var}(\mathbf{x})=-\frac{1}{2}(\eta^{1})^{-1},\;\mathbb{E}(\mathbf{x}\mathbf{x}^{T})=-\frac{1}{2}(\eta^{1})^{-1}+\frac{1}{4}(\eta^{1})^{-1}\eta^{2}(\eta^{2})^{T}(\eta^{1})^{-1}.$

The expression for its entropy is

$-\int\log
f(\mathbf{x}|\eta)f(\mathbf{x}|\eta)\mathrm{d}\mathbf{x}=\frac{1}{2}\log\det\left(2\pi
e\boldsymbol{\Omega}^{-1}\right)=-\frac{1}{2}\log\det\left(-\eta^{1}/(\pi
e)\right)\,.$

For gamma distribution with shape $\alpha$ and rate $\beta$, its density is

$\exp\Big{\\{}-\beta x+(\alpha-1)\ln
x+\alpha\ln\beta-\ln\Gamma(\alpha)\Big{\\}}$

with $t(x)=(x,\ln x)$ and canonical parameter
$\eta=(\eta^{1},\eta^{2})=(-\beta,\alpha-1)$. We also have
$\mathbb{E}(x)=\alpha/\beta=-(\eta^{2}+1)/\eta^{1}$.

For $p$ dimensional Wishart distribution with $n(>p-1)$ degrees of freedom and
scale matrix $\mathbf{V}\succ 0$, its density is

$\exp\Big{\\{}-\frac{1}{2}\mathrm{tr}\left(\mathbf{V}^{-1}\mathbf{X}\right)+\frac{n-p-1}{2}\log\det\mathbf{X}-\frac{np}{2}\log
2-\frac{n}{2}\log\det\mathbf{V}-\log\Gamma_{p}(n/2)\Big{\\}}$

with canonical parameter
$\eta=(\eta^{1},\eta^{2})=(-\frac{1}{2}\mathbf{V}^{-1},(n-p-1)/2)$. The mean
is
$\mathbb{E}\mathbf{X}=n\mathbf{V}=-0.5\cdot(2\eta^{2}+p+1)\cdot(\eta^{1})^{-1}$.

##  S.2 Proof of Theorem 3.1

With fixed
$q(\boldsymbol{\alpha}),q(\boldsymbol{\Omega}),q(\boldsymbol{\gamma})$,
updating
$q(\boldsymbol{\theta}_{i},f_{i})=q(\boldsymbol{\theta}_{i}|f_{i})q(f_{i})$ to
maximize ELBO in (17) can be equivalently expressed by maximizing

$\displaystyle\int\log\frac{p(f_{i},\boldsymbol{\theta}_{i}|\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma},\mathcal{D}_{i})}{q(f_{i},\boldsymbol{\theta}_{i})}\times
q(f_{i},\boldsymbol{\theta}_{i})q(\boldsymbol{\alpha})q(\boldsymbol{\Omega})q(\boldsymbol{\gamma})\
\,\mathrm{d}\boldsymbol{\theta}_{i}\,\mathrm{d}f_{i}\,\mathrm{d}\boldsymbol{\alpha}\,\mathrm{d}\boldsymbol{\Omega}\,\mathrm{d}\boldsymbol{\gamma}$
$\displaystyle=$ $\displaystyle\int h(f_{i})q(f_{i})\,\mathrm{d}f_{i}-\int\log
q(f_{i})q(f_{i})\,\mathrm{d}f_{i}\,,$ (S.1)

where

$\displaystyle h(f_{i})=$
$\displaystyle\int\log\frac{p(f_{i},\boldsymbol{\theta}_{i}|\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma},\mathcal{D}_{i})}{q(\boldsymbol{\theta}_{i}|f_{i})}\times
q(\boldsymbol{\theta}_{i}|f_{i})q(\boldsymbol{\alpha})q(\boldsymbol{\Omega})q(\boldsymbol{\gamma})\
\,\mathrm{d}\boldsymbol{\theta}_{i}\,\mathrm{d}\boldsymbol{\alpha}\,\mathrm{d}\boldsymbol{\Omega}\,\mathrm{d}\boldsymbol{\gamma}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f_{i})}\mathbb{E}_{q(\boldsymbol{\alpha})}\mathbb{E}_{q(\boldsymbol{\Omega})}\mathbb{E}_{q(\boldsymbol{\gamma})}\left[\log\frac{p(\boldsymbol{\theta}_{i}|f_{i},\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma},\mathcal{D}_{i})}{q(\boldsymbol{\theta}_{i}|f_{i})}+\log
p(f_{i}|\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma})\right]\,.$
(S.2)

To find the optimal $q(\boldsymbol{\theta}_{i},f_{i})$, we firstly find the
optimal $q(\boldsymbol{\theta}_{i}|f_{i})$ at each fixed $f_{i}$. Then from
equation (S.1), $q(f_{i})$ is readily available as
$q(f_{i})\propto\exp(h(f_{i}))$.

Notice

$\displaystyle\log
p(f_{i},\boldsymbol{\theta}_{i}|\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma},\mathcal{D}_{i})$
$\displaystyle\propto-\frac{1}{2}\sum_{b=1}^{B}(\mathbf{y}_{ib}-\mathbf{C}_{ib}\boldsymbol{\theta}_{ib})^{T}\boldsymbol{\Sigma}_{ib}^{-1}(\mathbf{y}_{ib}-\mathbf{C}_{ib}\boldsymbol{\theta}_{ib})-\frac{1}{2}\log\det\boldsymbol{\Sigma}_{ib}$
$\displaystyle\qquad\qquad-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0})^{T}\boldsymbol{\Theta}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0}),$

where the irrelevant constants have been dropped. Notice in the above, the
basis matrix $\mathbf{C}_{ib}$, the covariance matrix
$\boldsymbol{\Sigma}_{ib}$ as well as the prior mean vector
$\boldsymbol{\theta}_{0}$ all depend on the frequency $f_{i}$. The frequency
has uniform prior. The full conditional for $\boldsymbol{\theta}_{i}$ is

$\displaystyle
p(\boldsymbol{\theta}_{i}|f_{i},\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma},\mathcal{D}_{i})\propto\exp$
$\displaystyle\Big{\\{}\sum_{b=1}^{B}\left\langle-\frac{1}{2}\mathbf{C}_{ib}^{T}\boldsymbol{\Sigma}_{ib}^{-1}\mathbf{C}_{ib},\boldsymbol{\theta}_{ib}\boldsymbol{\theta}_{ib}^{T}\right\rangle+\sum_{b=1}^{B}\left\langle\mathbf{C}_{ib}^{T}\boldsymbol{\Sigma}_{ib}^{-1}\mathbf{y}_{ib},\boldsymbol{\theta}_{ib}\right\rangle\Big{\\}}$
$\displaystyle\times\exp\Big{\\{}-\frac{1}{2}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0})^{T}\boldsymbol{\Theta}(\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{0})\Big{\\}}$
$\displaystyle\propto\exp$ $\displaystyle\Big{\\{}\left\langle\xi_{i}^{\theta
1},\boldsymbol{\theta}_{i}\boldsymbol{\theta}_{i}^{T}\right\rangle+\left\langle\xi_{i}^{\theta
2},\boldsymbol{\theta}_{i}\right\rangle\Big{\\}}.$

This full conditional is Gaussian with canonical parameters

$\displaystyle\xi^{\theta 1}_{i}$
$\displaystyle=-\frac{1}{2}\times\mathrm{diag}\left(\mathbf{C}_{i1}^{T}\boldsymbol{\Sigma}_{i1}^{-1}\mathbf{C}_{i1},\cdots,\mathbf{C}_{iB}^{T}\boldsymbol{\Sigma}_{iB}^{-1}\mathbf{C}_{iB}\right)-\frac{1}{2}\boldsymbol{\Theta}\,,$
$\displaystyle\xi^{\theta 2}_{i}$
$\displaystyle=\left((\mathbf{C}_{i1}^{T}\boldsymbol{\Sigma}_{i1}^{-1}\mathbf{y}_{i1})^{T},\cdots,(\mathbf{C}_{iB}^{T}\boldsymbol{\Sigma}_{iB}^{-1}\mathbf{y}_{iB})^{T}\right)^{T}+\boldsymbol{\Theta}\boldsymbol{\theta}_{0}\,.$

For fixed $f_{i}$, $q(\boldsymbol{\theta}_{i}|f_{i})$ is imposed as a Gaussian
distribution with canonical parameter $(\eta_{i}^{\theta 1},\eta_{i}^{\theta
2})$. Holding
$q(\boldsymbol{\alpha}),q(\boldsymbol{\gamma}),q(\boldsymbol{\Omega})$ and
$f_{i}$ fixed, the optimal $q(\boldsymbol{\theta}_{i}|f_{i})$ maximizing (S.2)
is found by setting its canonical parameter as $\eta^{\theta
1}_{i}=\mathbb{E}_{q(\boldsymbol{\alpha})}\mathbb{E}_{q(\boldsymbol{\gamma})}\mathbb{E}_{q(\boldsymbol{\Omega})}(\xi^{\theta
1}_{i})$ and $\eta^{\theta
2}_{i}=\mathbb{E}_{q(\boldsymbol{\alpha})}\mathbb{E}_{q(\boldsymbol{\gamma})}\mathbb{E}_{q(\boldsymbol{\Omega})}(\xi^{\theta
2}_{i})$. This is the updating equation (18).

In order to compute $q(f_{i})\propto\exp(h(f_{i}))$, we now need to evaluate
$h(f_{i})$, which is

$h(f_{i})=\
\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f)}\mathbb{E}_{q(\boldsymbol{\alpha})}\mathbb{E}_{q(\boldsymbol{\gamma})}\log
p(f_{i},\boldsymbol{\theta}_{i}|\mathbf{y}_{i},\boldsymbol{\alpha},\boldsymbol{\Omega},\boldsymbol{\gamma})-\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f)}q(\boldsymbol{\theta}_{i}|f)\,.$

Then, with irrelevant constants dropped, we have

$\displaystyle h(f_{i})\propto$ $\displaystyle\ \langle\eta_{i}^{\theta
1},\,\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f_{i})}(\boldsymbol{\theta}_{i}\boldsymbol{\theta}_{i}^{T})\rangle+\langle\eta_{i}^{\theta
2},\,\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f_{i})}(\boldsymbol{\theta}_{i})\rangle-\frac{1}{2}\langle\mathbb{E}_{q(\boldsymbol{\Omega})}\mathbb{E}_{q(\boldsymbol{\gamma})}\boldsymbol{\Theta},\,\mathbb{E}_{q(\boldsymbol{\alpha})}\boldsymbol{\theta}_{0}\boldsymbol{\theta}_{0}^{T}\rangle$
$\displaystyle\qquad\qquad-\frac{1}{2}\sum_{b=1}^{B}\left[\mathbf{y}_{ib}^{T}\boldsymbol{\Sigma}_{ib}^{-1}\mathbf{y}_{ib}+\log\det\boldsymbol{\Sigma}_{ib}\right]-\frac{1}{2}\ln\det(-\eta_{i}^{\theta
1}/(\pi e)).$ (S.3)

The last term is the entropy of the Gaussian distribution
$q(\boldsymbol{\theta}_{i}|f_{i})$. In addition, the first and second moments
of $\boldsymbol{\theta}_{i}$ can be evaluated,

$\displaystyle\langle\eta_{i}^{\theta
1},\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f_{i})}(\boldsymbol{\theta}_{i}\boldsymbol{\theta}_{i}^{T})\rangle+\langle\eta_{i}^{\theta
2},\mathbb{E}_{q(\boldsymbol{\theta}_{i}|f_{i})}(\boldsymbol{\theta}_{i})\rangle$
$\displaystyle=$ $\displaystyle\langle\eta_{i}^{\theta
1},-\frac{1}{2}(\eta_{i}^{\theta 1})^{-1}+\frac{1}{4}(\eta_{i}^{\theta
1})^{-1}\eta_{i}^{\theta 2}(\eta_{i}^{\theta 2})^{T}(\eta_{i}^{\theta
1})^{-1}\rangle+\langle\eta_{i}^{\theta 2},-\frac{1}{2}(\eta_{i}^{\theta
1})^{-1}\eta_{i}^{\theta 2}\rangle$ $\displaystyle=$
$\displaystyle-\frac{3B}{2}-\frac{1}{4}(\eta_{i}^{\theta
2})^{T}(\eta_{i}^{\theta 1})^{-1}\eta_{i}^{\theta 2}.$ (S.4)

Now, (S.3) and (S.4) together imply (19).

##  S.3 Discussion of parameterization of the non-periodic component of a
light curve

In our formulation of the non-periodic component of a light curve (Section
2.4), $h_{ib}(\cdot)$ is a function of phase $u=f_{i}\cdot t$. This differs
from the work of He et al. (2016) and Yuan et al. (2018), where the term
$h_{ib}(\cdot)$ is chosen as a function of time $t$. This section provides
some numerical results supporting the new formulation.

Our choice of the model is guided by the goal of using the same kernel length-
scale parameter $\tau_{b2}$ for all Miras so that the computationally
expensive process of finding the optimal length-scale parameter individually
for each Mira can be avoided. Equation (10) reads

$s_{ib}(t)=m_{ib}+\mathbf{b}(f_{i}\cdot
t)^{T}\boldsymbol{\beta}_{ib}+h_{ib}(f_{i}\cdot t),$ (S.5)

where $\mathbf{b}(u)=(\cos(2\pi u),\sin(2\pi u))^{T}$. Note that the argument
in the periodic component $\mathbf{b}(f_{i}\cdot
t)^{T}\boldsymbol{\beta}_{ib}$ is $f_{i}\cdot t$. Our choice of using
$h_{ib}(f_{i}\cdot t)$ makes the function argument of the non-period component
consistent with the period component. This choice is supported by two
numerical results given below.

First result. We fit the Gaussian process model to the $I$\- band light curves
of OGLE MLC O-rich Miras under the two parameterizations ($h_{ib}(t)$ versus
$h_{ib}(f_{i}t)$) of the non-period component. We used the kernel given in
Section 2.4 of the paper and compared the optimal fitted length-scale
parameters. Because these MLC light curves are densely sampled with high
signal-to-noise ratio, we can fit the Gaussian process model directly to each
Mira, without using the hierarchical model structure developed in our work
(each Mira has its own length-scale parameter).

Figure S.1: The fitted length-scale parameters $\tau_{b2}$ (on the $\log_{10}$
scale) versus the logarithm of the periods for O-rich LMC Mira $I$-band light
curves. The model fitting is carried out independently to 469 OGLE LMC O-rich
Miras. The left penal shows the fitted $\tau_{b2}$ of the unscaled model
$h_{ib}(t)$. The right penal shows the fitted $\tau_{b2}$ of the scaled model
$h_{ib}(f_{i}\cdot t)$.

Figure S.2: The frequency estimation comparison of using unscaled version
$h_{ib}(t)$ (the left penal) versus the scaled version $h_{ib}(f_{i}\cdot t)$
(the right penal) in the Beyesian hierarchical model. In each penal, the
horizontal axis is the estimated frequency, and the vertical axis is the true
frequency.

The left and right panels of Figure S.1 shows respectively the results of
using the unscaled model $h_{ib}(t)$ and the scaled model $h_{ib}(f_{i}\cdot
t)$. The vertical axis (on the $\log_{10}$ scale) is the optimal fitted
$\tau_{b2}$, and the horizontal axis is the logarithm of the period. The left
panel shows that $\tau_{b2}$ varies in a large range, invalidating the idea of
letting all Miras share the same length-scale parameter, while the right panel
shows much smaller scattering of $\tau_{b2}$. The positive correlation between
$\tau_{b2}$ and $\log_{10}(p)$ is evident on the left panel but much less so
on the right panel. This result suggests that the scaled model
$h_{ib}(f_{i}\cdot t)$ is preferred if we let all Miras share the same length-
scale parameter in the Gaussian process model of the non-periodic component.

Second result. Using the simulated data for 5,000 O-rich Miras presented in
Section 4, we tested the unscaled version $h_{ib}(t)$ with the shared kernel
parameter for all Miras, and keeping the rest of procedures the same in our
Bayesian hierarchical model. It turns out the model performance deteriorates
considerably. More precisely, when $h_{ib}(t)$ is used instead of
$h_{ib}(f_{i}\cdot t)$, the recovery rate (RR) as defined in the paper drops
from $98.04\%$ to $28.04\%$, and the absolute deviation error (ADE) increases
from $0.52\times 10^{-4}$ to $8.51\times 10^{-4}$.

Figure S.2 depicts the estimation result. The left penal shows the result of
using $h_{ib}(t)$, while the right penal presents the result of using
$h_{ib}(f_{i}\cdot t)$. In each penal, the horizontal axis is the estimated
frequency, and the vertical axis is the true frequency. It is obvious that the
left panel shows a large scattering, while on the right panel most points
reside on the diagonal. This result strongly supports our use of the scaled
model $h_{ib}(f_{i}\cdot t)$.

##  S.4 Additional Simulation Study

In addition to the simulation in the main text, we have conducted a more
extensive simulation. We generated 90 datasets of O-rich Mira light curves,
containing observations in the $I$ and $K_{s}$ bands. In a given dataset, each
star has exactly $n_{I}$ and $n_{K}$ data points in the $I$ and $K_{s}$ bands,
respectively. The number of data points $(n_{K},n_{I})$ takes one of the
following ten choices: $(5,5)$, $(5,10)$, $(5,20)$, $(5,30)$, $(10,10)$,
$(10,20)$, $(10,30)$, $(20,20)$, $(20,30)$ and $(30,30)$.

Each dataset chooses one of three patterns (cadence) to generate time points.
The observation time points are generated by the cadence of the M33 dataset
(C1), the cadence of the OGLE LMC dataset (C2), or simply by uniform
distribution (C3). For the first cadence C1, the $K_{s}$-band data is obtained
at a later time than the $I$-band data. For the second cadence C2, the two
bands are observed roughly at the same time but have a seasonal pattern, such
that there are no observations for a range of months each year when the Sun is
too close to the galaxy. As for the last cadence, the two bands have similar
uniform coverage over the whole observational time interval.

Given the observation time points, the light curve magnitudes are also
generated from the template of Yuan et al. (2017). In addition, Gaussian noise
is added to the magnitudes. According to real data from the actual survey, the
noise level increases as the object becomes fainter (with larger magnitude
values). A commonly-used function to fit the magnitude-noise relation is
$\sigma(m)=\exp\left(a\cdot m^{c}-b\right)$. We adopted the values
$a=1.82\times 10^{-6}$, $b=5.84$ and $c=5$ that closely match the OGLE survey.
This is our first noise setting (N1). The other two settings increase the
noise level by using $\sigma(m)+0.05$ (N2) and $\sigma(m)+0.1$ (N3),
respectively. Figure S.3 illustrates the three noise levels. In the figure,
the grey points represent the actual noise-mag relation from the OGLE survey.
They are 20,000 randomly-selected observation points from OGLE. For our second
simulation, the signal-to-noise ratio is lower than that of the OGLE survey,
as the curves are slightly higher than the points. Therefore, it is more
difficult to recover the correct periods.

Figure S.3: The magnitude noise relation for Simulation II. The three noise
levels (N1, N2, and N3) are plotted as solid, dashed and dotted curves,
respectively.

In summary, we have 10 choices for the number of data points, 3 choices of
sampling cadence, and 3 choices of noise level. By all possible combinations,
the above procedure produces 90 datasets. Figure S.4 illustrates a few
generated light curves for the choice of sample size $n_{K}=30$ and
$n_{I}=20$. The first, second and third rows correspond to sampling cadences
C1, C2 and C3, respectively, while the first, second, and third columns
correspond to noise levels N1, N2 and N3, respectively. The black points and
blue squares denote $I$\- and $K_{s}$-band light curves, respectively.

Figure S.4: Randomly selected synthetic Mira data for Simulation 2. The first,
second and third rows correspond to sampling cadence C1, C2 and C3,
respectively. The first, second, and third columns correspond to noise level
N1, N2 and N3, respectively. Black points are $I$-band light curves, and blue
squares are $K_{s}$-band light curves.

We apply the methods SP, MSP, GLS, MGLS and SVI to the simulated datasets. In
these datasets, the single-band methods are applied to the $I$-band data. The
same metrics as in the simulation Section 4 of the paper, RR and ADE, are used
for model comparison. Figure S.5 shows the recovery rate for the 90 datasets.
The horizontal axis represents the varying sample sizes $(n_{K},n_{I})$ for
$K_{s}$\- and $I$-band data. Its columns correspond to different noise levels
(N1–3 from left to right), and its rows for various sampling cadences (C1–3
from up to bottom). Similarly, Figure S.6 shows the result for the ADE, with
its vertical axes on the logarithm scale. It is remarkable that the proposed
method SVI has the best performance, and consistently estimates period with
higher accuracy for all of the datasets. For this simulation, the single-band
methods SP and GLS have a similar performance. However, MGLS outperforms MSP
in both metrics. This suggests the MSP might not be robust for a dataset with
slightly different configurations from its model setting.

##  S.5 Visual assessment of light curve fitting for Miras in M33

The quality of light curve fitting using Gaussian process can be visually
assessed by plotting the MAP estimator. Figure S.7 plots the data and fitted
curves for six randomly selected Miras in M33. The light curves are fitted by
our model with the double exponential kernel. The model provides reasonable
fits to the data. In particular, for the sparse M33 light curves, the
$J,H,K_{s}$ bands show smaller variability beyond the sinusoidal periodical
component. The $I$-band light curves have larger variability beyond the
sinusoidal periodical component, and the corresponding kernel parameter
$\tau_{b1}$ is larger than that of the other three bands. The kernel we used
is

$k_{b}(u,u^{\prime})=\tau_{b1}\cdot\exp\big{\\{}-(u-u^{\prime})^{2}/\tau_{b2}\big{\\}}+\tau_{b3}\cdot\mathrm{I}(u=u^{\prime}).$

The table below summaries the kernel parameter fitted to each band of the real
data. The values of $\tau_{b2}$ (or $\tau_{b3}$) among different bands are
close, the minor difference may not be easily explained.

band $b$ | $\tau_{b1}$ | $\tau_{b2}$ | $\tau_{b3}$
---|---|---|---
$I$ | $e^{-3.33}$ | $e^{-4.20}$ | $e^{-4.21}$
$J$ | $e^{-10.48}$ | $e^{-3.00}$ | $e^{-5.28}$
$H$ | $e^{-10.00}$ | $e^{-2.57}$ | $e^{-3.90}$
$K_{s}$ | $e^{-10.59}$ | $e^{-4.39}$ | $e^{-4.18}$
Figure S.5: The recovery rate for Simulation II.  Figure S.6: The absolute
deviation error for Simulation II.

Figure S.7: Light curve fitting for six randomly selected Miras in M33. The
double exponential kernel is used for the Gaussian process.