  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.
