with $\|\bm{\alpha}_{ij}\|_{\mathrm{F}}\leq 1$, then for each $\mu_{i}$ it
should satisfy $\mu_{i}\geq\|\bm{\Theta}_{ij}\|_{\mathrm{F}}$ for any $j$ with
$C(j)\neq C(i)$. By further plugging the definition of $\mu_{i}$ we get
$\bar{\mu}\geq\max\left\\{\bar{\mu}_{1},\bar{\mu}_{2}\right\\}$ with
$\bar{\mu}_{1},\bar{\mu}_{2}$ defined in Lemma B.2, and we set
$\bar{\mu}=\max\left\\{\bar{\mu}_{1},\bar{\mu}_{2}\right\\}$. This completes
our guess of dual variables. ∎

###### Proof of Lemma B.3.

We follow the same route as in Lemma 3.4 by first rewriting $\bm{A}$,
$\bm{\Lambda}$ and $\bm{\Theta}$ into four blocks as (63). Also, let
$\widetilde{\bm{\mu}}_{1}=\text{diag}(\\{\mu_{i}\bm{I}_{d}\\}_{i=1}^{m_{1}})\in\mathbb{R}^{m_{1}d\times
m_{1}d}$ and
$\widetilde{\bm{\mu}}_{2}=\text{diag}(\\{\mu_{i}\bm{I}_{d}\\}_{i=m_{1}+1}^{n})\in\mathbb{R}^{m_{2}d\times
m_{2}d}$. Then the four blocks of $\bm{\Theta}$ satisfy

$\displaystyle\widetilde{\bm{\Theta}}_{11}=\frac{\widetilde{\bm{\mu}}_{1}\widetilde{\bm{M}}_{1}^{*}}{\sqrt{d}},\quad\quad\quad\widetilde{\bm{\Theta}}_{22}=\frac{\widetilde{\bm{\mu}}_{2}\widetilde{\bm{M}}_{2}^{*}}{\sqrt{d}},$
$\displaystyle\widetilde{\bm{\Theta}}_{12}+\widetilde{\bm{\Theta}}_{21}^{\top}=\frac{1}{m_{1}}\widetilde{\bm{M}}_{1}^{*}\widetilde{\bm{A}}_{21}+\frac{1}{m_{2}}\widetilde{\bm{A}}_{12}\widetilde{\bm{M}}_{2}^{*}-\frac{1}{m_{1}m_{2}}\widetilde{\bm{M}}_{1}^{*}\widetilde{\bm{A}}_{12}\widetilde{\bm{M}}_{2}^{*}.$

The remaining proof is very similar to the one of Lemma 3.4, and we do not
repeat. ∎

###### Proof.

According to Lemma B.3, $\lambda_{\text{min}}((p-\epsilon)\bm{I}_{nd})$
satisfies

$\displaystyle\lambda_{\text{min}}((p-\epsilon)\bm{I}_{nd})=p-\epsilon=p-\max\\{\bar{\mu}_{1},\;\bar{\mu}_{2}\\}-\sigma^{(1)}-\sigma^{(2)}$
$\displaystyle=p\\!+\\!\min\bigg{\\{}\underbrace{2\min_{i\in
C_{1}}\bigg{\\{}\sigma_{i}\\!-\\!\frac{m_{1}\gamma}{\sqrt{d}}\max_{j\in
C_{2}}\|\widetilde{\bm{\alpha}}_{ij}\|_{\mathrm{F}}\bigg{\\}}\\!-\\!\sigma^{(1)}}_{=:\epsilon_{1}},\;\underbrace{2\min_{i\in
C_{2}}\bigg{\\{}\sigma_{i}\\!-\\!\frac{m_{2}(1\\!-\\!\gamma)}{\sqrt{d}}\max_{j\in
C_{1}}\|\widetilde{\bm{\alpha}}_{ji}\|_{\mathrm{F}}\bigg{\\}}\\!-\\!\sigma^{(2)}}_{=:\epsilon_{2}}\bigg{\\}}$
$\displaystyle=p+\min\\{\epsilon_{1},\epsilon_{2}\\}.$

Then we bound $\epsilon_{1}$ and $\epsilon_{2}$ separately as

$\displaystyle\epsilon_{1}$ $\displaystyle\geq 2\min_{i\in
C_{1}}\sigma_{i}-\frac{2m_{1}\gamma}{\sqrt{d}}\max_{i\in C_{1},j\in
C_{2}}\|\widetilde{\bm{\alpha}}_{ij}\|_{\mathrm{F}}-\sigma^{(1)},$
$\displaystyle\epsilon_{2}$ $\displaystyle\geq 2\min_{i\in
C_{2}}\sigma_{i}-\frac{2m_{2}(1-\gamma)}{\sqrt{d}}\max_{i\in C_{2},j\in
C_{1}}\|\widetilde{\bm{\alpha}}_{ji}\|_{\mathrm{F}}-\sigma^{(2)}.$

For $\min_{i\in C_{1}}\sigma_{i}$ and $\min_{i\in C_{2}}\sigma_{i}$, we can
use the bounds derived in (74) and (75) respectively. For $\max_{i\in
C_{1},j\in C_{2}}\|\widetilde{\bm{\alpha}}_{ij}\|_{\mathrm{F}}$, by using the
result in (71) we obtain

$\max_{i\in C_{1},j\in
C_{2}}\|\widetilde{\bm{\alpha}}_{ij}\|_{\mathrm{F}}\leq\sqrt{2(1+c)\beta}\bigg{(}\frac{1}{\sqrt{\rho}}+\frac{1}{\sqrt{1-\rho}}\bigg{)}\bigg{(}\frac{\log
n}{n}\bigg{)}+o\bigg{(}\frac{\log n}{n}\bigg{)}$

with probability $1-n^{-c}$. For $\sigma^{(1)}$, since $r_{ij}=r_{ji}$ for
$i,j\in C_{1}$, it satisfies

$\sum\limits_{s_{1}\in C_{1}}\sum\limits_{s_{2}\in
C_{1}}r_{s_{1}s_{2}}=2\sum\limits_{s_{1},\;s_{2}\in
C_{1},\;s_{1}<s_{2}}r_{s_{1}s_{2}},$

and $\sum\limits_{s_{1},s_{2}\in C_{1},\;s_{1}<s_{2}}r_{s_{1}s_{2}}$ follows a
binomial distribution $\mathrm{Binom}\left(m_{1}(m_{1}-1)/2,\alpha\log
n/n\right)$. By applying Hoeffding’s inequality [6] we obtain

$\displaystyle\mathbb{P}\Bigg{\\{}\sum\limits_{s_{1},s_{2}\in
C_{1},\;s_{1}<s_{2}}r_{s_{1}s_{2}}-\frac{(m_{1}-1)\alpha\rho\log n}{2}\geq
t\Bigg{\\}}\leq\exp\left(-\frac{4t^{2}}{m_{1}^{2}}\right),$

for any $t>0$. This further indicates $\sigma^{(1)}\leq\alpha\rho\log n+o(\log
n)$ with high probability. Similarly we can obtain
$\sigma^{(2)}\leq\alpha(1-\rho)\log n+o(\log n)$ with high probability. By
putting all the results together and let $c\rightarrow 0$ to be close to zero,
we get the bounds for $\epsilon_{1},\epsilon_{2}$ and further
$\lambda_{\text{min}}((p-\epsilon)\bm{I}_{nd})$. ∎

## References

  * [1] Emmanuel Abbe.  Community detection and stochastic block models: recent developments.  The Journal of Machine Learning Research, 18(1):6446–6531, 2017\.
  * [2] Emmanuel Abbe, Afonso S Bandeira, and Georgina Hall.  Exact recovery in the stochastic block model.  IEEE Transactions on Information Theory, 62(1):471–487, 2015.
  * [3] Greg W Anderson, Alice Guionnet, and Ofer Zeitouni.  An introduction to random matrices.  Cambridge university press, 2010.
  * [4] Chandrajit Bajaj, Tingran Gao, Zihang He, Qixing Huang, and Zhenxiao Liang.  SMAC: simultaneous mapping and clustering using spectral decompositions.  In International Conference on Machine Learning, pages 324–333, 2018.
  * [5] Afonso S Bandeira and Ramon Van Handel.  Sharp nonasymptotic bounds on the norm of random matrices with independent entries.  Annals of Probability, 44(4):2479–2506, 2016.
  * [6] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration inequalities: A nonasymptotic theory of independence.  Oxford university press, 2013.
  * [7] Nicolas Boumal.  Nonconvex phase synchronization.  SIAM Journal on Optimization, 26(4):2355–2377, 2016.
  * [8] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [9] Aurelien Decelle, Florent Krzakala, Cristopher Moore, and Lenka Zdeborová.  Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications.  Physical Review E, 84(6):066106, 2011.
  * [10] Patrick Doreian, Vladimir Batagelj, and Anuska Ferligoj.  Generalized blockmodeling, volume 25.  Cambridge university press, 2005.
  * [11] Martin E. Dyer and Alan M. Frieze.  The solution of some random NP-hard problems in polynomial expected time.  Journal of Algorithms, 10(4):451–489, 1989.
  * [12] Yifeng Fan, Tingran Gao, and Zhizhen Zhao.  Unsupervised co-learning on $G$-manifolds across irreducible representations.  In Advances in Neural Information Processing Systems, pages 9041–9053, 2019.
  * [13] Yifeng Fan, Tingran Gao, and Zhizhen Zhao.  Representation theoretic patterns in multi-frequency class averaging for three-dimensional cryo-electron microscopy.  Information and Inference: A journal of IMA, page accepted, 2021\.
  * [14] Yifeng Fan and Zhizhen Zhao.  Multi-frequency vector diffusion maps.  In International Conference on Machine Learning, pages 1843–1852. PMLR, 2019.
  * [15] Stephen E Fienberg, Michael M Meyer, and Stanley S Wasserman.  Statistical analysis of multiple sociometric relations.  Journal of the american Statistical association, 80(389):51–67, 1985.
  * [16] Santo Fortunato.  Community detection in graphs.  Physics reports, 486(3-5):75–174, 2010.
  * [17] Joachim Frank.  Three-Dimensional Electron Microscopy of Macromolecular Assemblies: Visualization of Biological Molecules in Their Native State.  Oxford University Press, New York, 2nd edition, 2006.
  * [18] Tingran Gao and Zhizhen Zhao.  Multi-frequency phase synchronization.  In International Conference on Machine Learning, pages 2132–2141. PMLR, 2019.
  * [19] Bruce Hajek, Yihong Wu, and Jiaming Xu.  Achieving exact cluster recovery threshold via semidefinite programming.  IEEE Transactions on Information Theory, 62(5):2788–2797, 2016\.
  * [20] Bruce Hajek, Yihong Wu, and Jiaming Xu.  Achieving exact cluster recovery threshold via semidefinite programming: Extensions.  IEEE Transactions on Information Theory, 62(10):5918–5937, 2016\.
  * [21] Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  Social networks, 5(2):109–137, 1983.
  * [22] Brian Karrer and Mark EJ Newman.  Stochastic blockmodels and community structure in networks.  Physical review E, 83(1):016107, 2011.
  * [23] Rafał Latała, Ramon van Handel, and Pierre Youssef.  The dimension-free structure of nonhomogeneous random matrices.  Inventiones mathematicae, 214(3):1031–1080, 2018.
  * [24] Roy R Lederman and Amit Singer.  A representation theory perspective on simultaneous alignment and classification.  Applied and Computational Harmonic Analysis, 2019.
  * [25] Jing Lei and Alessandro Rinaldo.  Consistency of spectral clustering in stochastic block models.  The Annals of Statistics, 43(1):215–237, 2015.
  * [26] Xiaodong Li, Yudong Chen, Jiaming Xu, et al.  Convex relaxation methods for community detection.  Statistical Science, 36(1):2–15, 2021.
  * [27] Shuyang Ling.  Near-optimal performance bounds for orthogonal and permutation group synchronization via spectral methods.  arXiv preprint arXiv:2008.05341, 2020.
  * [28] Laurent Massoulié.  Community detection thresholds and the weak ramanujan property.  In Proceedings of the forty-sixth annual ACM symposium on Theory of computing, pages 694–703, 2014.
  * [29] Frank McSherry.  Spectral partitioning of random graphs.  In Proceedings 42nd IEEE Symposium on Foundations of Computer Science, pages 529–537. IEEE, 2001.
  * [30] Elchanan Mossel, Joe Neeman, and Allan Sly.  Reconstruction and estimation in the planted partition model.  Probability Theory and Related Fields, 162(3):431–461, 2015.
  * [31] Elchanan Mossel, Joe Neeman, and Allan Sly.  A proof of the block model threshold conjecture.  Combinatorica, 38(3):665–708, 2018.
  * [32] Mark EJ Newman.  The structure and function of complex networks.  SIAM review, 45(2):167–256, 2003.
  * [33] Andy Nguyen, Mirela Ben-Chen, Katarzyna Welnicka, Yinyu Ye, and Leonidas Guibas.  An optimization approach to improving collections of shape maps.  In Computer Graphics Forum, volume 30, pages 1481–1491. Wiley Online Library, 2011.
  * [34] Jiming Peng and Yu Wei.  Approximating $k$-means-type clustering via semidefinite programming.  SIAM journal on optimization, 18(1):186–205, 2007.
  * [35] Amelia Perry and Alexander S Wein.  A semidefinite program for unbalanced multisection in the stochastic block model.  In 2017 International Conference on Sampling Theory and Applications (SampTA), pages 64–67. IEEE, 2017.
  * [36] John William Strutt Baron Rayleigh.  The theory of sound, volume 2.  Macmillan, 1896.
  * [37] Lord Rayleigh.  XII. on the resultant of a large number of vibrations of the same pitch and of arbitrary phase.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 10(60):73–78, 1880.
  * [38] Yunpeng Shi and Gilad Lerman.  Message passing least squares framework and its application to rotation synchronization.  arXiv preprint arXiv:2007.13638, 2020.
  * [39] Amit Singer.  Angular synchronization by eigenvectors and semidefinite programming.  Applied and computational harmonic analysis, 30(1):20–36, 2011\.
  * [40] Amit Singer and H-T Wu.  Vector diffusion maps and the connection laplacian.  Communications on pure and applied mathematics, 65(8):1067–1144, 2012.
  * [41] Amit Singer, Zhizhen Zhao, Yoel Shkolnisky, and Ronny Hadani.  Viewing angle classification of cryo-electron microscopy images using eigenvectors.  SIAM Journal on Imaging Sciences, 4(2):723–759, 2011.
  * [42] Gilbert W Stewart.  Perturbation theory for the singular value decomposition.  Technical report, University of Maryland, College Park, 1998.
  * [43] Terence Tao.  Topics in random matrix theory, volume 132.  American Mathematical Soc., 2012.
  * [44] Joel A… Tropp.  User-friendly tail bounds for sums of random matrices.  Foundations of computational mathematics, 12(4):389–434, 2012.
  * [45] Joel A… Tropp.  An introduction to matrix concentration inequalities.  Foundations and Trends® in Machine Learning, 8(1-2):1–230, 2015.
  * [46] Lieven Vandenberghe and Stephen Boyd.  Semidefinite programming.  SIAM review, 38(1):49–95, 1996.
  * [47] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [48] Bowei Yan, Purnamrita Sarkar, and Xiuyuan Cheng.  Provable estimation of the number of blocks in block models.  In International Conference on Artificial Intelligence and Statistics, pages 1185–1194. PMLR, 2018.
  * [49] Zhizhen Zhao and Amit Singer.  Rotationally invariant image representation for viewing direction classification in cryo-EM.  Journal of structural biology, 186(1):153–166, 2014.


11in

Vol.0 (20xx) No.0, 000–000

11institutetext: Department of Physics, M. K. Bhavnagar University, Bhavnagar
- 364002, Gujarat, India<EMAIL_ADDRESS>

Received 2019 March 24; accepted 2019 April 09

# Time-resolved CCD photometry and time-series analysis of RR Lyrae type RR
Gem star

V. A. Gohil and S.P. Bhatnagar

###### Abstract

We present the results of a time-resolved photometric and time-series analysis
of a RR Lyrae type star RR Gem. The main results are as follows: We found RR
Gem’s pulsation period, 0.39689 d, and its V and I mean magnitudes, 11.277 (V)
and 11.063 (I) mag respectively. We confirm its variability type as RRab/BL
because of its showing the Bla$\breve{z}$ko effect and it also shows
asymmetric light curves (steep ascending branches), periods from 0.3 to 1.0
days, and amplitudes from 0.3 to 2 mag. in V. They are fundamental mode
pulsators.

###### keywords:

Stars: individual: RR Gem — stars: variables: RR Lyr — stars: oscillations —
stars: variables: RRAB/BL — stars: variables: Blazhko effect — techniques:
photometric.

##  1 Introduction

RR Lyrae variables are older pulsating white giants with low metallicity. They
are common in globular clusters-dense groups of old stars in the halos of
galaxies. Like Cepheids, their pulsations are periodic. RR Lyraes have 0.5
solar mass and have a short pulsation period of 0.2 to 0.85 days and amplitude
variations of 0.3 to 2 magnitudes. RR Lyrae stars are usually spectral class
A. The main sequence and the HB(Horizontal Branch) do not intersect. RR Lyrae
stars are on the HB, indeed where the instability strip intersects the HB. The
HB stars have left the red giant branch and are characterized by helium fusion
in their cores surrounded by a shell of hydrogen fusion. All of the RR Lyrae
stars in a cluster have the same average apparent magnitude. In different
clusters, the average apparent magnitude is different. This is because all RR
Lyraes have about the same average absolute magnitude of +0.75.

RR111https://www.aavso.org/vsx/index.php?view=help.vartype&nolayout=1&abbrev=RR
means that this is the RR variable in the Lyrae constellation, which are
radially-pulsating giant A-F stars having amplitudes from 0.2 to 2 mag. in V
band. Cases of variable light-curve shapes as well as variable periods are
known. If their amplitudes are modified and also their periods, then they are
called the ”Bla$\breve{z}$ko effect”(Bla$\breve{z}$ko 1907)(denoted by the sub
type BL). The majority of these stars belong to the spherical component of the
Galaxy; they are present, sometimes in large numbers, in some globular
clusters, where they are known as pulsating horizontal-branch stars. Like
Cepheids, maximum expansion velocities of surface layers for these stars
practically coincide with maximum light.

RRAB/BL222https://www.aavso.org/vsx/index.php?view=help.vartype&nolayout=1&abbrev=RRAB%2FBL
type variable stars are RR Lyrae variables with asymmetric light curves (steep
ascending branches), periods from 0.3 to 1.0 days, and amplitudes from 0.5 to
2 mag. in V band. They are fundamental mode pulsators. BL indicates RR Lyrae
stars showing the Bla$\breve{z}$ko effect.

J. Jurcsik and her team (Jurcsik 1990, 1998; Jurcsik et al. 2005, 2006a,
2006b, 2006c, 2008a, 2008b; Jurcsik 2009a; Jurcsik et al. 2009b, c) has done a
remarkable work on the Blazhko behavior of RR Gem Star and in particular RR
Lyr type stars showing Bla$\breve{z}$ko effects. Most RR Lyrae stars repeat
their light curves with remarkable regularity. About 30 % of the known
galactic RR Lyrae stars, however, display cyclic modulation in shape and
amplitude of the light curve over tens of pulsation cycles.

##  2 Methodology

###  2.1 Instrumentation

M. K. Bhavnagar University’s Kumari Aanya Binoy Gaardi
Observatory333https://www.google.com/maps/place/Bhavnagar+Observatory/@21.7537555,72.1302003,
18z/data=!4m5!3m4!1s0x395f509a719b3a2b:0x30504b9860ad6057!8m2!3d21.7542085!
4d72.1303511?hl=en was used for this study. The Fully automated Observatory
contains a 14 inch Celestron Schmidt-Cassegrain reflector type Telescope
mounted on a German Equatorial mount. Telescope has an SBIG STF-7 CCD Camera
along with SBIG CFW-8 filter wheel with Johnson-Cousins _UBVRI_
444http://www.aip.de/en/research/facilities/stella/instruments/data/johnson
-ubvri-filter-curves photometric filters inbuilt into the CCD are attached on the telescope. Telescope is manual as well as fully computer controlled and telescope has two stepper motors for RA and DEC axis which is controlled by stepper motor controllers indigenously developed in collaboration with IUCAA, Pune(Inter University Center for Astronomy and Astrophysics). Open source software SCOPE is used to control the stepper motor controller and Cartes du Ciel(Sky Chart) was used for the CCD control as well as to give positions commands to SCOPE software which send command to stepper motor controller. Stepper motor controller performs the PWM and send electrical signals to RA and DEC stepper motors. There is an USB connection between CCD and PC for data as well as commands. Astronomy data is locally stored in Windows XP based PC. A GPS receiver is also connected with Windows PC which periodically updates time of the PC which in turn updates the DOC based PC’s clock and thus the whole system’s time gets updated periodically. Manual sliding rooftop is connected with a motor through pulley for opening and closing down the Observatory’s roof.

###  2.2 Observations and data reduction

####  2.2.1 Differential Photometry

MaxIm DL555http://diffractionlimited.com/product/maxim-dl/ photometric image
analysis software was used for our data reduction work. Each and every time-
series astronomical image was first calibrated for basic noise removal and
then every continuous sequence was aligned and then photometry was performed
on them.

For doing differential photometry(Rodrigez 2005) we needed our object star(RR
Gem)666https://www.aavso.org/vsx/index.php?view=detail.top&oid=14304, one or
more check stars and one reference star(standard star). Standard observatory
procedures were used and used with minimum air-mass considerations.

For data reduction(calibration), flats field frames, dark frames and bias
frames as well as light frames(See Figure 1.) were required so we programmed
it into our MaxIm DL software and then it’s done automatically according to
our exposure time and other Johnson-Cousins ubvri filter value settings. For
flat fielding we used daily twilight images and after all these we made master
files of each of these and used them in our final calibration part.

In MaxIm DL we combined all the flats, darks and bias files and combined them
into single standard master files according to different photometry filers.
For calibration we used these master files. We first opened all the object
star’s light images and loaded the master calibration files and performed
calibrate using calibrate all menu. By this we get all the files calibrated,
then we again opened the calibrated light images and align them using the
align tool.

Now for performing Photometry we again opened all the light images and using
the photometry menu’s differential photometry tool we performed the photometry
on the light images by selecting each individual object star, reference star
and check stars. Finally it gave us results in a Magnitude Vs JD(Jullian Day)
graph in CSV or AAVSO Text format file.

####  2.2.2 Time-series analysis

Vstar(AAVSO)777https://www.aavso.org/vstar program was used for performing
time-series analysis (Rodrigez 2005). We used DC DFT with Period Range tool
and fed the proper values in it like low period as 0.1 d, high period as 0.5 d
and resolution as 0.0001 and we got the periodogram with period 0.39689 d.

We downloaded the standard published data for our star from AAVSO data
tool888https://www.aavso.org/data-download for the period of 10 months to get
the bigger view of our star and performed the same above mentioned process and
got the results(See Figures 2 and 3) then we also downloaded about 1 months
data near to our observing period and analyzed the same and got this
results(See Figures 4 and 5). By comparing this one month’s data with our
data(See Figures 6, 7 for (V) and 8, 9 for (I) respectively) we found it to be
nearly perfectly matching with the magnitude as well as with the period. This
ascertains that our observations and data reduction process is standardized.

Figure 1: Clear Filter Object, Comparison and Check Star Field-CLEAR Filter:
16-04-2014 Figure 2: RR GEM 10 Months Data 2014 Johnson V Periodogram Figure
3: RR GEM 10 Months Phase Data 2014 Johnson V Phase Graph Figure 4: RR GEM 1
Months Data 2014 Johnson V Periodogram Figure 5: RR GEM 1 Months Phase Data
2014 Johnson V Phase Graph Figure 6: RR GEM 15 Days Our Data 2014 Johnson V
Periodogram Figure 7: RR GEM 15 Days Phase Our Data 2014 Johnson V Phase Graph
Figure 8: RR GEM 15 Days Our Data 2014 Johnson I Periodogram Figure 9: RR GEM
15 Days Phase Our Data 2014 Johnson I Phase Graph

According to AAVSO(VSX)999https://www.aavso.org/vsx/index.php?view=search.top
its spectral type is A9-F6. It was discovered by L. Ceraski (Ceraski, L. 1903)
and all these information is based on 09 Feb 1972 (HJD 2441357.205). Its AAVSO
UID : 000-BBM-940. Its other names are AAVSO 0715+31, AN 13.1903, HIP 35667.
Our data is from JD: 2456764.16483 (16 APR 2014) to JD: 2456801.11748 (23 MAY
2014) totaling 15 nights. During this time the mean magnitudes in visual
filter and infrared filters were 11.277 (V) and 11.063 (I) respectively.

We also compared our results with other astronomical databases. The details
are as shown in Table 1.

Database | ID / Name | J2000.0 | Type | V mag | V mag | Period | Epoch | Spectral Class101010https://www.handprint.com/ASTRO/specclass.html
---|---|---|---|---|---|---|---|---
|  | ($RA$) | ($DEC$) |  | ($Max$) | ($Min$) | ($d$) |  |
VSX111111Variable Star Index(AAVSO). | 000-BBM-940 | 07:21:33.53 | 30:52:59.50 | RRab/BL | 10.62 | 11.99 | 0.3973106 | 2441357.205(HJD) | A9-F6
ASAS-SN121212The All Sky Automated Survey for SuperNovae. | ASASSN-VJ072133.66+305259.7 | 07:21:33.66 | 30:52:59.70 | RRab | - | 11.34 | 0.3992869 | 2457690.06(HJD) | -
CDS131313Strasbourg astronomical Data Center(CDS-Strasbourg). | V*̂ RR Gem | 07:21:33.53 | 30:52:59.46 | Variable Star of RR Lyr type | - | 11.92 | - | - | A8
GCVS141414General Catalog of Variable Stars. | RR Gem | 07:21:33.53 | 30:52:59.50 | RRab | 10.62 | 11.99 | 0.3973106 | 41357.21 | A9-F6
GVSSG151515General Variable Star Search Gateway. | RR Gem / 380010 | 07:21:33.50 | 30:52:59.50 | RRab | 10.62 | 11.99 | 0.3973106 | 41357.21(JD) | A9-F6

Table 1: RR Gem in different Astronomical databases.

As shown in Table 1, we compared our results with many astronomical databases
like VSX, ASAS-SN161616https://asas-
sn.osu.edu/database/light_curves/195828,CDS171717http://cdsportal.u-strasbg.fr/?target=rr%20gem,GCVS181818http://www.sai.msu.su/gcvs/cgi-
bin/search.cgi?search=RR+Gem and General Variable Star Search
Gateway191919http://var2.astro.cz/gsg/index.php?star=rr+gem&all=yes&alldata=yes&oejv=yes&gcvs
=yes&nsv=yes&brka=yes&meka=yes&czev=yes&bcvs=yes&dssplate=yes&usecoords=GCVS&rezim=
search_now. We found our mean magnitude in V band: 11.277 mag and in I band:
11.063 mag, which is quite close to the published databases values like 11.185
mag in V band with the difference between our results and VSX(AAVSO) with the
minor difference of 0.092, which is quite close.

We also compared our period with the astronomical databases and found that our
period 0.39689 d is very close to the published value 0.3973106 d VSX(AAVSO)
with the minor difference of 0.0004206 d, which is also quite close.

####  2.2.3 Comparision and Check Stars

We used Check and Comparison
Stars202020https://www.aavso.org/apps/vsp/photometry/?fov=18.5&scale=F&star=rr+gem
&orientation=ccd&maglimit=18.5&resolution=150&north=up&east=left&type=chart
from Field
photometry212121https://www.aavso.org/apps/vsp/photometry/?fov=18.5&scale=F&star=rr+gem
&orientation=ccd&maglimit=18.5&resolution=150&north=up&east=left&type=photometry
for RR Gem from AAVSO Variable Star Database:

Star type | AUID | J2000.0 | Magnitudes
---|---|---|---
|  | ($RA$) | ($DEC$) | ($B$) | ($V$) | ($B-V$) | ($Rc$) | ($Ic$)
Standard / Comparison Star | 000-BJR-395 | 07:21:40.41 | 30:52:23.80 | 9.560(0.020)*̂ | 9.412(0.016)*̂ | 0.148(0.026) | 9.322(0.020)*̂ | 9.219(0.021)*̂
Check Star | 000-BBM-938 | 07:21:33.34 | 30:54:43 | 10.652(0.016)*̂ | 10.271(0.014)*̂ | 0.381(0.021) | 10.043(0.017)*̂ | 9.826(0.019)*̂

Table 2: Comparison and Check Star Table for RR Gem.

The comparison / standard star and check stars were selected from the AAVSO’s
field photometry section. They were selected because they were very near, in
the same CCD field. They were chosen because they were about the same
magnitude as well as were within the CCD field. All the details are given in
the table 2. Practically they should be near to your object star and also
their magnitude as well as spectral class should also match. Comparison star
should be a standard non-variable star and check star could be non-variable or
variable for the aperture(differential) photometry. So, they were chosen
accordingly and we found them to be of satisfactory nature as per our
requirements.

##  3 Results

RR Gem’s mean magnitude for V was 11.277 (V) mag and for I was 11.063 (I) mag.
We had also compared star’s data from AAVSO Variable star database for near
about the same time period and found its visible filter’s mean magnitude value
to be 11.185 V, which gives the difference of 0.092 V which is also quite
negligible.

After performing time-series analysis it was found that the period of RR Gem
star published in the VSX(AAVSO) database was 0.3973106 d and after our
analysis we found it to be 0.39689 d, the difference between the two is
0.0004206 d which is quite low. Thus it establishes our analysis matches with
the published literature values.

##  4 Conclusion

We have established that our observations reported in this paper are in very
good agreement with the published data of the star RR Gem and that it is an RR
Lyrae type star and its variability type is RRab/BL and is clearly showing the
Blazhko effect. We also conclude that our observatory can be now used for
future studies of variable stars and in particular will be useful for
performing optical light curve data and related analysis for the upcoming
transient alerts like GW and Gamma-ray bursts which require optical followups
from ground based telescopes during the bright phases (at various geographic
longitudes).

###### Acknowledgements.

We would like to acknowledge the M. K. Bhavnagar University for providing
Astronomical research facility at Kumari Aanya Binoy Gardi Observatory and We
would also like to thank, IUCAA(Inter University Center For Astronomy And
Astrophysics), AAVSO(American Association for Variable Star Observers),
VSX(Variable Star Index)(AAVSO).

## References

  * Bla$\breve{z}$ko (1907) Bla$\breve{z}$ko, 1907, AN, 175, 325
  * Jurcsik (1990) Jurcsik J., Kovacs G., 1990, A & A, 312, 111
  * Jurcsik (1998) Jurcsik J., Johanna, 1998, A & A, 333, 571
  * Jurcsik et al. (2005) Jurcsik J., Sódor, $\acute{A}$., V$\acute{a}$radi, M., Szeidl, B., Washuettl, A., Weber, M., D$\acute{e}$k$\acute{a}$ny, I., Hurta,
Zs., Lakatos, B., Posztob$\acute{a}$nyi, K., 2005, A &A, 430, 1049

  * Jurcsik et al. (2006a) Jurcsik J., Szeidl, B., Sódor, $\acute{A}$., D$\acute{e}$k$\acute{a}$ny, I., Hurta, Zs., Posztob$\acute{a}$nyi, K., Vida, K., V$\acute{a}$radi, M.,
Szing, A., 2006, AJ, 132, 61

  * Jurcsik et al. (2006b) Jurcsik J., Sódor, $\acute{A}$., V$\acute{a}$radi, M., Vida, K., Posztobanyi, K., Szing, A., Hurta, Zs., Dekany, I.,
Washuettl, A., Vityi, N., 2006, , 5709, 1

  * Jurcsik et al. (2006c) Jurcsik J., Szeidl, B., V$\acute{a}$radi, M., Henden, A., Hurta, Zs., Lakatos, B., Posztob$\acute{a}$nyi, K.,
Klagyivik, P., Sódor, $\acute{A}$., 2006, A & A, 445, 617

  * Jurcsik et al. (2008a) Jurcsik J., Sódor, $\acute{A}$., Hurta, Zs., Kovari, Zs., Vida, K., Hajdu, G., Nagy, I., Dekany, I.,
Posztob$\acute{a}$nyi, K., Koponyas, B., 2008, , 5844, 1

  * Jurcsik et al. (2008b) Jurcsik J., Sódor, $\acute{A}$., Hurta, Zs., V$\acute{a}$radi, M., Szeidl, B., Smith, H. A., Henden, A., Dékány, I.,
Nagy, I., Posztob$\acute{a}$nyi, K., 2008, MNRAS, 391, 164

  * Jurcsik (2009a) Jurcsik J., 2009, CoAst, 159, 53
  * Jurcsik et al. (2009b) Jurcsik J., Sódor, $\acute{A}$., Szeidl, B., Koll$\acute{a}$th, Z., Smith, H. A., Hurta, Zs., V$\acute{a}$radi, M., Henden, A.,
D$\acute{e}$kány, I., Nagy, I., 2009, MNRAS, 393, 1553

  * Jurcsik et al. (2009c) Jurcsik J., Sódor, $\acute{A}$., Szeidl, B., Hurta, Zs., V$\acute{a}$radi, M., Posztob$\acute{a}$nyi, K., Vida, K., Hajdu, G.,
Kõv$\acute{a}$ri, Zs., Nagy, I., 2009, MNRAS, 400, 1006

  * Rodrigez (2005) P. Rodrigez-Gil, 2005, A & A, 431, 289


# Neural classifiers based Monte Carlo simulation

Elouan Argouarc’h, François Desbouvries,
SAMOVAR
Télécom SudParis, Institut Polytechnique de Paris
91120 Palaiseau France
{elouan.argouarch<EMAIL_ADDRESS>

###### Abstract

Acceptance-rejection (AR), Independent Metropolis Hastings (IMH) or importance
sampling (IS) Monte Carlo (MC) simulation algorithms all involve computing
ratios of probability density functions (pdfs). On the other hand, classifiers
discriminate labellized samples produced by a mixture density model, i.e., a
convex linear combination of two pdfs, and can thus be used for approximating
the ratio of these two densities. This bridge between simulation and
classification techniques enables us to propose (approximate) pdf-ratios-based
simulation algorithms which are built only from a labellized training data
set.

##  1 Introduction

If $a$ and $b$ are two positive numbers,

$\displaystyle
r=\frac{a}{a+b}\in(0,1)\Leftrightarrow\frac{r}{1-r}=\frac{a}{b}>0.$ (1)

This humble algebraic identity has interesting consequences in Bayesian
classification, machine learning and stochastic simulation.

Indeed, if $a$ and $b$ are probabilities of a sample in a binary mixture
context, then ratio $\frac{a}{a+b}$ is the classical formula for the posterior
probability. This posterior provides with the class probabilities given a
sample, and can be approximated by a parametric classifier $r^{\star}\approx
r$ trained to distinguish between the two probability distributions. On the
other hand, positive ratios $\frac{a}{b}$ play a key role in AR, IMH or
sampling importance resampling (SIR) techniques. The r.h.s. of (1) relates $r$
to such positive ratios, and tells us that ratio $\frac{a}{b}$ can be computed
exactly from $r$, or, in practice, approximately from $r^{\star}$, without
necessarily knowing $a$ nor $b$. This simple observation will enable us to
propose versions of these algorithms which rely on relaxed weaker hypotheses.

Let $\lambda\in(0,1)$ and $(1-\lambda)$ be the prior probabilities of two
categories $k=1,0$, distributed resp. $\sim$ $p_{1}$ and $p_{0}$. Binary
classification tries to distinguish samples from mixture $\lambda
P_{1}+(1-\lambda)P_{0}$ by identifying the pdf which generated them. The
appropriate way to classify relies on the posterior probability: $x$ is a
sample $\sim$ $P_{1}$ rather than $\sim$ $P_{0}$ with probability

$\mathbbm{p}(k=1|x,\lambda,P_{0},P_{1})=\frac{\lambda p_{1}(x)}{\lambda
p_{1}(x)+(1-\lambda)p_{0}(x)}.$ (2)

As is well known (see e.g. [1, Chap. 11]), assigning a sample to the label
with highest posterior probability is the optimal decision rule in the sense
that it minimizes the probability of misclassification.

To compute this posterior probability, one needs to evaluate the pdfs
$p_{1},p_{0}$ associated with $P_{1}$ and $P_{0}$ and know the prior
probability $\lambda$. Unfortunately, $\lambda$ is often unknown so (2) is
untractable. If however we dispose of a set
$\mathcal{D}=\\{(x_{i},k_{i})\\}_{i=1}^{N_{0}+N_{1}}$ of labelled
observations, $\lambda$ can be estimated by $\frac{N_{1}}{N_{1}+N_{0}}$ where
$N_{1}$ and $N_{0}$ are respectively the number of samples from $P_{1}$ and
$P_{0}$. This leads to the (approximate) probability:

$\mathbbm{p}(k=1|x,\mathcal{D},P_{0},P_{1})=\frac{N_{1}p_{1}(x)}{N_{1}p_{1}(x)+N_{0}p_{0}(x)}.$
(3)

However in most cases, $p_{0}$ and $p_{1}$ are unknown too, so (3) cannot be
computed either. When we only dispose of $\mathcal{D}$, we can make use of a
parametric classifier (in this paper we call classifier any parameterized
function $r_{\phi}(x)$ which mimics the unknown posterior pdf). So let us
assume that we have at our disposal a function $r_{\phi}$ such that

$r_{\phi}(x)\approx\frac{N_{1}p_{1}(x)}{N_{1}p_{1}(x)+N_{0}p_{0}(x)}.$ (4)

Our paper is based on the observation that (4) is equivalent to

$\frac{N_{0}}{N_{1}}\frac{r_{\phi}(x)}{1-r_{\phi}(x)}\approx\frac{p_{1}(x)}{p_{0}(x)},$
(5)

which implies that (typically neural network based) classifiers can also be
used for approximating density ratios.

Equation (5) has already been observed, and exploited in contexts where
estimating a ratio of pdfs is relevant. First, classifiers are at the core of
adversarial training techniques in which divergence measures involving a ratio
are replaced by an approximation based on a classifier [2]. This enables to
learn implicit generative models (i.e., with untractable pdfs) [3] [4] [5].
Moreover, classifier based density ratio approximation have been applied to
estimation of such metrics as Mutual Information [6]. Finally, classifiers
based ratios have been applied successfully in statistical hypothesis testing
procedures [7], which heavily rely on likelihood-ratio tests.

If now $p_{0}$ is an instrumental distribution, then (5) can be turned into an
(approximately normalized) approximation of target pdf $p_{1}$. So classifiers
can be used for density estimation, conditional density estimation, or
likelihood-to-evidence ratio estimation, making them especially relevant in a
likelihood-free inference setting [8][9][10].

However, the question of sampling the corresponding model remains open, and
this is precisely the point we discuss in this paper. We realize that pdf
ratios also play a key role in such simulation techniques as the AR or Markov
Chain Monte Carlo (MCMC) methods, in which an acceptance probability needs to
be computed, or in importance sampling, which corrects samples from a
surrogate distribution by a weight function which is also a ratio of two
densities. This establishes a connection between classification and MC
sampling, and will enable us to relax the assumptions underlying those
sampling algorithms, at the price of approximate sampling.

The rest of this paper is organized as follows. In §2 we recall classical
ratio-based stochastic simulation algorithms, i.e. the AR, IMH and IS
techniques. In §3 we show that classifiers computed via the Binary Cross
Entropy (BCE) criterion indeed provide with a suitable approximation of the
posterior (3). Finally in §4 we propose classification based sampling methods,
and illustrate our method via simulations. We end the paper with a conclusion.

##  2 Classical ratio-based sampling algorithms

Stochastic simulation includes a variety of techniques, see e.g. [11]-[16]. In
this section we focus on AR, IMH and IS which share in common that they all
compute a ratio of pdfs.

###  2.1 The AR algorithm

####  2.1.1 A brief reminder of AR Sampling

AR Sampling [13, chap. 2] is a simulation algorithm that yields samples
distributed according to a target distribution $P$ via samples from a proposal
distribution $Q$, which are accepted or rejected as valid samples from $P$ via
some acceptance probability. More precisely, let the support of $P$ be inside
that of $Q$. This means that there exists a constant $C\geq 1$ s.t. for all
$y\in\mathbbm{R}^{d},p(y)\leq Cq(y)$. Let $Y\sim Q$, and let $k$ a Bernoulli
random variable (r.v.) with parameter $\alpha(y)=\frac{p(y)}{Cq(y)}$. AR
sampling is based on the fact that $Y|k=1$ is distributed according to $P$.

Note that $\mathbbm{p}(k=1)=\frac{1}{C}$, so the lower the value of $C$, the
higher the acceptance rate.

In order to use the algorithm in practice, we thus need to know pdf $p$, and
build $Q$ s.t. one can sample easily from $Q$, pdf $q$ is available, there
exists $C$ such that $p(y)\leq Cq(y)$ for all $y$, we can compute one such
value of $C$, and $C$ is as small as possible. Note finally that the algorithm
can easily be relaxed to the cases where $p$ and/or $q$ are known up to a (non
necessarily common) constant, see e.g. [15, Th. 4.5].

####  2.1.2 Revisiting AR sampling as optimal binary classification

As we shall now see, AR sampling is indeed nothing but a binary classification
procedure (see also [17, §6] for an application of this principle in the
context of a mixture model for blind source separation). Starting from the
target pdf $p(x)$, we find an easy-to-sample distribution $Q$ and constant
$C>1$ s.t. $Cq(x)$ envelopes $p(x)$. Since $Cq(x)-p(x)$ is non negative, we
write $Cq(x)$ as $p(x)$ plus a positive reminder which, up to a constant, is
also a pdf; so enveloping $p(x)$ with $Cq(x)$ (or indeed $\frac{1}{C}p(x)$
with $q(x)$) is nothing but building the implicit binary mixture distribution
(see also figure 1 below)

$\underbrace{q(x)}_{\rm proposal}=\frac{1}{C}\;\underbrace{p(x)}_{\rm
target}\;+\;(1-\frac{1}{C})\;\underbrace{\frac{q(x)-\frac{1}{C}p(x)}{1-\frac{1}{C}}}_{\rm
reminder}\;$ (6)

with a priori classes probabilities $\frac{1}{C}$ and $1-\frac{1}{C}$. The
first component of the mixture (that of interest) is the target density $p$,
and the second one is the law of the rejected samples.

Figure 1: Envelopping target pdf builds an implicit mixture.

The magic of the AR algorithm consists in drawing samples from mixture $Q$
without needing to sample from its two components (see the r.h.s. of (6)).
Accepting (or rejecting) a sample depending on the ratio probability

$\alpha(x)=\frac{p(x)}{Cq(x)}=\frac{\frac{1}{C}p(x)}{\frac{1}{C}p(x)+(1-\frac{1}{C})\frac{q(x)-\frac{1}{C}p(x)}{1-\frac{1}{C}}}$
(7)

then amounts to classifying the samples with the (optimal [1, Chap. 11])
posterior pdf (compare (7) to (2)).

###  2.2 IMH

MCMC algorithms build a Markov chain whose invariant distribution is the
target distribution $P$; so simulating the chain yields samples asymptotically
distributed $\sim$ $P$. The Metropolis-Hastings (MH) algorithm [13] [14] is a
particular MCMC method which constructs the Markov Chain $x_{t}$ as a two-step
procedure: given a current state $x_{t}$, the algorithm draws a candidate
$x^{*}$ from a proposal distribution $q(.|x_{t})$, and then calculates the
acceptance probability
$\alpha(x^{*},x_{t})=\min(1,\frac{p(x^{*})q(x_{t}|x^{*})}{p(x_{t})q(x^{*}|x_{t})})$.
$x^{*}$ is accepted as the new state $x_{t+1}$ with probability
$\alpha(x^{*},x_{t})$; if $x^{*}$ is rejected then the chain remains in the
current state $x_{t}$. In practice, $q(.|x_{t})$ plays a crucial role in the
performance of the MH algorithm: if not well-tuned, the acceptance rate may be
too low, leading to slow mixing of the chain, or too high, leading to poor
exploration of the target distribution.

The IMH algorithm is a simplified version of MH which considers an independent
transition. The new point $x^{*}$ is hence proposed independently of the
current state $x_{t}$, according to an independent proposal $q(.)$. In this
case, the acceptance probability simplifies to
$\alpha(x^{*},x_{t})=\min(1,\frac{p(x^{*})q(x_{t})}{p(x_{t})q(x^{*})})$.

###  2.3 IS

In many signal processing problems we want to compute the expectation of some
function $f$ with respect to density $p$: $\mu=\int f(x)p(x){\rm
d}x=\mathbbm{E}_{P}\mathopen{}\left[{f(x)}\,\negthickspace\right]\mathclose{}.$
In practice $\mu$ can be very difficult to compute, so one often needs to
resort to approximations. IS is a variance reduction technique for integral MC
estimates which can be traced back to the 1950’s [18] [19] [16, §5.4].

The crude MC estimate of $\mu$ reads $\hat{\mu}^{\rm
MC}=\frac{1}{N}\sum_{i=1}^{N}f(x_{i})$ with $x_{i}\stackrel{{\scriptstyle\rm
iid}}{{\sim}}P$. However it is generally difficult to sample directly from
$P$, moreover $\hat{\mu}^{\rm MC}$ can be a poor estimate, particularly when
the regions where $p$ is large do not coincide with those where $f$ is large.
Rewriting $\mu$ as $\mu=\int f(x)\frac{p(x)}{q(x)}q(x){\rm d}x,$ where $Q$ is
some importance distribution, leads to the IS estimator $\hat{\mu}^{\rm
IS}(q)=\frac{1}{N}\sum_{i=1}^{N}\frac{p(x_{i})}{q(x_{i})}f(x_{i}),\,\,\,x_{i}\stackrel{{\scriptstyle\rm
iid.}}{{\sim}}Q.$ As far as variance reduction is concerned, one can easily
show that the importance pdf which minimizes
$\mathbbm{V}\textrm{ar}(\hat{\mu}^{\rm IS}(q))$ is $q_{\rm opt}^{\rm
IS}(x)\propto|f(x)|p(x)$. Even if in practice $\hat{\mu}^{\rm IS}(q_{\rm
opt}^{\rm IS})$ cannot be computed, this tells us that the regions where it is
important to sample from (whence the term "importance distribution") are not
those where $p$ is large, but rather those where $|f|p$ is large. Note that
$\hat{\mu}^{\rm IS}$ can be computed only if $p$ and $q$ are known exactly, or
known up to a common constant; if this is not the case one can resort to self-
normalized IS [20].

Besides being a variance reduction technique, IS can also be seen as a two
step sampling procedure for producing samples (approximatively) drawn from
$p$, out of samples originally drawn from $q$. The technique is known as
Rubin’s SIR mechanism [21], [22], [23], [24, §9.2]: Let
$\\{x_{i}\\}_{i=1}^{N}$ be $N$ iid. samples from $q(x)$, and given
$\\{x_{i}\\}_{i=1}^{N}$, let $\\{\tilde{x}^{i}\\}_{i=1}^{M}$ be $M$ iid.
samples from
$\sum_{i=1}^{N}\frac{p(x_{i})/q(x_{i})}{\sum_{i=1}^{N}p(x_{i})/q(x_{i})}\delta_{x_{i}}({\rm
d}x)$ (in other words, we draw samples from $q$, weight each of them with
weight $w_{i}$ proportional to $p(x_{i})/q(x_{i})$, and resample $M$ iid.
points from this random discrete mass probability function). Then
$\\{\tilde{x}^{i}\\}_{i=1}^{M}$ are dependent and are not $p$-distributed, but
become iid. samples from $p({x})$ if $N\rightarrow\infty$.

##  3 Parametric classifier by minimizing the BCE

From now on we consider the setting where $\lambda$, $p_{1}$ and $p_{0}$ are
unknown, and we only have a set ${\cal D}$ of samples from $P_{0}$ and
$P_{1}$, see §1. In this context, we should build a parametric function
$r_{\phi}(x)$ that approximates the posterior pdf from the recorded samples.
The aim of this section is to show that minimizing a BCE criterion indeed
yields such a suitable approximation, since the BCE, up to constants, is
nothing but an MC approximation of a Kullback-Leibler Divergence
($D_{\mathrm{KL}}$) between the classifier and the unavailable posterior pdf.

To see this, let us first recall the BCE criterion:

$\mathcal{L}_{\mathrm{BCE}}(\phi)=-\sum_{i=1}^{N_{1}}\log(r_{\phi}(x^{(1)}_{i}))-\sum_{i=1}^{N_{0}}\log(1-r_{\phi}(x^{(0)}_{i})),$
(8)

where $r_{\phi}(x)=\mathbbm{p}_{\phi}(k=1|x)$ is the probability (under model
$\phi$) that the label associated to an observation $x$ is $1$. Let $h(x,k)$
be the joint distribution over observations and labels:

$\displaystyle
h(x,k)\\!=\\!\underbrace{\frac{N_{1}^{k}N_{0}^{1-k}}{N_{1}+N_{0}}}_{h(k)}\underbrace{p_{1}(x)^{k}p_{0}(x)^{1-k}}_{h(x|k)},x\in\mathbbm{R}^{d},k=0,1.$
(9)

On the other hand, using
$r_{\phi}(x)\in\mathopen{}\left[{0,1}\,\negthickspace\right]\mathclose{}$, we
construct another joint probability distribution
$h_{\phi}(x,k)=h(x)r_{\phi}(x)^{k}(1-r_{\phi}(x))^{1-k}$, where $h(x)$ is the
$x$-marginal in (9). As the name Cross-Entropy suggests, the BCE loss is, up
to additive and multiplicative constants, nothing but an MC approximation of
$D_{\mathrm{KL}}\mathopen{}\left({h(x,k)||h_{\phi}(x,k)}\,\negthickspace\right)\mathclose{}=\mathbbm{E}_{h(x)}(D_{\mathrm{KL}}(h(k|x)||r_{\phi}(x)^{k}(1-r_{\phi}(x))^{1-k})$,
see appendix.

The interest of this interpretation is that, as is well known, a
$D_{\mathrm{KL}}$ reaches zero if and only if the two joint distributions are
equal almost surely. So minimizing
$D_{\mathrm{KL}}\mathopen{}\left({h(x,k)||h_{\phi}(x,k)}\,\negthickspace\right)\mathclose{}$
(over an arbitrary function $r_{\phi}$) would ensure that
$\mathbbm{p}_{\phi}(k|x)=h(k|x)$ for all $x\in\mathbbm{R}^{d}$ and for
$k=1,0$, i.e. that the classifier reaches the target posterior pdf.

Of course, in practice, minimizing the BCE does not ensure that this
$D_{\mathrm{KL}}$ decreases to zero. First, since we only dispose of a finite
number of observations and associated labels, minimizing an MC approximation
of the $D_{\mathrm{KL}}$ does not minimize the $D_{\mathrm{KL}}$ itself. Next,
the parametric family of classifiers
$\mathopen{}\left\\{{r_{\phi}|\phi\in\Phi}\,\negthickspace\right\\}\mathclose{}$
may not contain the optimal $r^{*}$, and therefore we may only reach a
positive minimum of the $D_{\mathrm{KL}}$. Lastly, even if the family of
classifiers contained $r^{*}$, standard optimization techniques would only
guarantee convergence to a local minimum of the $D_{\mathrm{KL}}$. Finally in
practice, minimizing the BCE loss only provides with $\phi$ such that
$r_{\phi}$ satisfies (4), and thus approximates the unknown $r^{*}$.

##  4 Using a binary classifier for (approximate) Sampling

We now come to the heart of this paper. If $P_{1}$ is a distribution of
interest in an MC sampling setting, and $P_{0}$ a suitable easy-to sample
instrumental distribution \- be it the proposal distribution in AR, the
independent Markov transition Kernel in IMH, or the importance distribution in
IS, then the three sampling algorithms involve the density ratio
$\frac{p_{1}(x)}{p_{0}(x)}$. As explained in section 3, a parametric binary
classifier trained from a set of recorded labeled observations computes an
approximation of the unknown posterior distribution. However, remember that
(4) is equivalent to (5); we thus see that classifiers can also be used for
approximating density ratios of interest, which enables us to propose
approximate versions of the sampling algorithms based on this classifier-ratio
approximation, and thus to relax the tractable densities requirement, but at
the cost of approximate sampling. Of course, the closer $P_{0}$ is to $P_{1}$,
the more efficient the algorithms will be. However, here $P_{0}$ is supposed
to be given and hence, our problem is not (as usual) to adjust $P_{0}$ from a
given $P_{1}$, but to make the most of dataset $\mathcal{D}$ for fixed $P_{0}$
and $P_{1}$.

Assumptions.

$P_{1}$ is the distribution of interest and $P_{0}$ a fixed instrumental
distribution. We neither know $p_{1}(x)$ nor $p_{0}(x)$, even up to a
constant, but we dispose of a dataset $\mathcal{D}$ made of samples
$x^{(0)}_{1},...,x^{(0)}_{N_{0}}\sim P_{0}$ and
$x^{(1)}_{1},...,x^{(1)}_{N_{1}}\sim P_{1}$, and assume that we can train a
binary classification model $r_{\phi}$ which minimizes (8).

###  4.1 A classifier-based AR algorithm

Remember from §2.1.1 that the AR algorithm is all the more efficient that
$\frac{1}{C}P_{1}$ is close to $P_{0}$, and constant $C$ (which is $\geq 1$)
is as small as possible.

Indeed the minimum value of admissible $C$ is $C_{\rm
min}=\sup_{y\in\mathbbm{R}^{d}}\frac{p_{1}(y)}{p_{0}(y)}$. Now, a key
ingredient in the algorithm is the acceptance probability
$\alpha(x)=\frac{1}{C}\times\frac{p_{1}(x)}{p_{0}(x)},$ but in our setting
$p_{1}(x)$ and $p_{0}(x)$ are unavailable, so we cannot compute neither ratio
$\frac{p_{1}(x)}{p_{0}(x)}$, nor its supremum $C_{\rm min}$.

Following the idea expressed in (5), we can however make use of a classifier
for approximating ratio $\frac{p_{1}(x)}{p_{0}(x)}$, and approximate its
supremum over an uncountable set by using the data at hand:

$\alpha(x)\leftarrow\alpha_{\phi}(x)$ (10)
$\alpha_{\phi}(x)=\frac{1}{\widetilde{C}}\frac{r_{\phi}(x)}{1-r_{\phi}(x)}\text{
where }\widetilde{C}=\max_{y\in\mathcal{D}}\frac{r_{\phi}(y)}{1-r_{\phi}(y)}.$
(11)

The supremum is only defined on the finite set of samples in $\mathcal{D}$ and
is therefore easy to compute. Note that the approximation $\widetilde{C}$ can
also be updated effortlessly each time we propose a new sample from $P_{0}$.

Finally, if moreover proposal $P_{0}$ has a tractable pdf $p_{0}(x)$, then the
algorithm provides as a byproduct an estimate $p_{\phi}$ (up to a normalizing
constant) of the pdf of the accepted samples:

$p_{\phi}(x)=\frac{p_{0}(x)\alpha_{\phi}(x)}{\int_{\mathbbm{R}^{d}}p_{0}(z)\alpha_{\phi}(z)\mathrm{d}z}=\frac{p_{0}(x)\frac{r_{\phi}(x)}{1-r_{\phi}(x)}}{\int_{\mathbbm{R}^{d}}p_{0}(z)\frac{r_{\phi}(z)}{1-r_{\phi}(z)}\mathrm{d}z}.$
(12)

If required by the application, the normalizing constant can be estimated via
MC, i.e. by replacing the denominator with
$\frac{1}{M}\sum_{i=1}^{M}\frac{r_{\phi}(z_{i})}{1-r_{\phi}(z_{i})},z_{i}\sim
P_{0}$.

###  4.2 A classifier based IMH algorithm

Remember from section 2.2 that the IMH algorithm requires computing the
acceptance probability

$\alpha(x^{*},x_{t})=\min\mathopen{}\left({1,\frac{p_{1}(x^{*})/p_{0}(x^{*})}{p_{1}(x_{t})/p_{0}(x_{t})}}\,\negthickspace\right)\mathclose{},$
(13)

but in our setting we cannot compute $\alpha(x^{*},x_{t})$ since $p_{1}$ and
$p_{0}$ are unavailable. However, using (5) we can make use of the classifier
in order to compute an approximated acceptance probability:

$\displaystyle\alpha(x^{*},x_{t})$
$\displaystyle\leftarrow\alpha_{\phi}(x^{*},x_{t}),$ (14)
$\displaystyle\alpha_{\phi}(x^{*},x_{t})$
$\displaystyle=\min\mathopen{}\left({1,\frac{r_{\phi}(x^{*})}{(1-r_{\phi}(x^{*}))}\frac{(1-r_{\phi}(x_{t}))}{r_{\phi}(x_{t})}}\,\negthickspace\right)\mathclose{}.$
(15)

###  4.3 A classifier-based IS algorithm

Remember from section 2.3 that SIR relies on the normalized importance weights
$w(x^{(0)}_{i})=w^{\rm u}(x^{(0)}_{i})/(\sum_{j=1}^{N}w^{\rm
u}(x^{(0)}_{j})),$ where $\\{x^{(0)}_{j}\\}_{j=1}^{N}$ are samples from the
importance pdf $P_{0}$, and in which the unnormalized weights read
$p_{1}(x^{(0)}_{i})/p_{0}(x^{(0)}_{i})$. Here we cannot compute $w^{\rm
u}(x^{(0)}_{i})$ since $p_{1}$ and $p_{0}$ are unavailable. However, using (5)
we can make use of the classifier in order to compute an approximated
unormalized weight:

$\displaystyle w^{\rm u}(x^{(0)}_{i})\leftarrow w^{\rm
u}_{\phi}(x^{(0)}_{i})\;\text{where}\;w^{\rm
u}_{\phi}(x^{(0)}_{i})=\frac{r_{\phi}(x^{(0)}_{i})}{1-r_{\phi}(x^{(0)}_{i})}.$
(16)

###  4.4 Illustrating examples

Figure 2: Density ratio (middle-right) via classification of samples from $P$
(left) and $Q$ (middle-left) - approximate samples from $P$ (right) obtained
via a ratio based algorithm.

We illustrate our approach (see fig. 2 below) on reference 2D examples in
order to illustrate the mechanism of (i) obtaining an approximate of the
density ratio from samples, and (ii) sampling from the target distribution via
that density ratio. The instrumental $Q$ was set to be Gaussian with mean and
covariance estimated from the samples from $P$ (even though it can be
computed, pdf $q$ was not used during the procedure). Notice that the
classifier based sampler is indeed compatible with multimodal data, so long as
$Q$ covers the different regions of mass of $P$.

##  5 Conclusion

In this paper we showed that making use of (typically neural network based)
classifiers enables to relax the assumptions underlying some classical, pdfs-
ratio-based MC samplers. The main interest of the versions of the AR, IMH and
SIR samplers which we propose is that they only require a labellized training
data set; in particular, the knowledge of the two pdfs (possibly up to a
constant) involved in these samplers is no longer necessary. On the other
hand, our algorithms produce samples which are only approximately drawn from
the target density. Our methodology is confirmed by simulations.

$\displaystyle D_{\mathrm{KL}}(h(x,k)||h_{\phi}(x,k))$
$\displaystyle\\!\\!\\!=\\!\\!\\!$
$\displaystyle-\mathbbm{H}[(x,k)]-\mathbbm{E}_{h(\\!k,x)}[\log\mathopen{}\left({h(x)}\,\negthickspace\right)\mathclose{}]$
(17)
$\displaystyle-\mathbbm{E}_{h(k,x)}[\log\mathopen{}\left({\mathbbm{p}_{\phi}(k|x)}\,\negthickspace\right)\mathclose{}],$

in which $\mathbbm{H}[(x,k)]=\mathbbm{E}_{h(x,k)}[-\log(h(x,k))]$ is the
entropy of the joint r.v. $(x,k)$. In (17), only the last term depends on
$\phi$. We get the BCE loss by computing an MC approximation of that last term
(or, equivalently, replacing the expectation with one computed on the
empirical joint distribution):

$\displaystyle\mathbbm{E}_{h(k,x)}[\log(\mathbbm{p}_{\phi}(k|x))]\\!=\\!\sum_{k=0}^{1}h(k)\\!\\!\\!\int_{\mathbbm{R}^{d}}\\!\\!\\!\\!\\!\\!\log(\mathbbm{p}_{\phi}(k|x))h(x|k)\mathrm{d}x$
$\displaystyle\stackrel{{\scriptstyle\eqref{joint}}}{{=}}$
$\displaystyle\sum_{k=0}^{1}\frac{N_{1}^{k}N_{0}^{1-k}}{N_{1}+N_{0}}\int_{\mathbbm{R}^{d}}\log(\mathbbm{p}_{\phi}(k|x))h(x|k)\mathrm{d}x$
$\displaystyle\approx$
$\displaystyle\sum_{k=0}^{1}\frac{N_{1}^{k}N_{0}^{1-k}}{N_{1}+N_{0}}[\frac{1}{N_{k}}\sum_{i=1}^{N_{k}}\log(\mathbbm{p}_{\phi}(k|x_{i}^{(k)}))]$
$\displaystyle=$
$\displaystyle\frac{1}{N_{1}+N_{0}}(\sum_{i=1}^{N_{1}}\log(r_{\phi}(x^{(1)}_{i}))+\sum_{i=1}^{N_{0}}\log(1-r_{\phi}(x^{(0)}_{i})))$
$\displaystyle=$
$\displaystyle-\frac{1}{N_{1}+N_{0}}\mathcal{L}_{\mathrm{BCE}}(\phi).$

So
$D_{\mathrm{KL}}\mathopen{}\left({h(x,k)||h_{\phi}(x,k)}\,\negthickspace\right)\mathclose{}\approx
A+\frac{1}{N_{1}+N_{0}}\mathcal{L}_{\mathrm{BCE}}(\phi),$ in which $A$ does
not depend on $\phi$, and consequently
$\arg\min_{\phi}D_{\mathrm{KL}}\mathopen{}\left({h(x,k)||h_{\phi}(x,k)}\,\negthickspace\right)\mathclose{}$
$\approx$ $\arg\min_{\phi}\mathcal{L}_{\mathrm{BCE}}(\phi).$

## References

  * [1] C. Giraud, Introduction to High-Dimensional Statistics, vol. 139 of Monographs on Statistics and Applied Probability.  Chapman and Hall/CRC, Boca Raton, FL, 2014.
  * [2] X. Nguyen, M. J. Wainwright, and M. I. Jordan, “Estimating divergence functionals and the likelihood ratio by convex risk minimization,” IEEE Transactions on Information Theory, vol. 56, no. 11, pp. 5847–5861, 2010.
  * [3] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial networks,” Communications of the ACM, vol. 63, no. 11, pp. 139–144, 2020.
  * [4] A. Creswell, T. White, V. Dumoulin, K. Arulkumaran, B. Sengupta, and A. A. Bharath, “Generative adversarial networks: An overview,” IEEE Signal Processing Magazine, vol. 35, no. 1, pp. 53–65, 2018.
  * [5] A. Creswell and A. A. Bharath, “Denoising adversarial autoencoders,” IEEE transactions on neural networks and learning systems, vol. 30, no. 4, pp. 968–984, 2018.
  * [6] M. I. Belghazi, A. Baratin, S. Rajeshwar, S. Ozair, Y. Bengio, A. Courville, and D. Hjelm, “Mutual information neural estimation,” in International conference on machine learning, pp. 531–540, PMLR, 2018.
  * [7] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola, “A kernel two-sample test,” The Journal of Machine Learning Research, vol. 13, no. 1, pp. 723–773, 2012.
  * [8] C. Durkan, I. Murray, and G. Papamakarios, “On contrastive learning for likelihood-free inference,” in International conference on machine learning, pp. 2771–2781, PMLR, 2020.
  * [9] J. Hermans, V. Begy, and G. Louppe, “Likelihood-free MCMC with amortized approximate ratio estimators,” in International conference on machine learning, pp. 4239–4248, PMLR, 2020.
  * [10] O. Thomas, R. Dutta, J. Corander, S. Kaski, and M. U. Gutmann, “Likelihood-free inference by ratio estimation,” Bayesian Analysis, vol. 17, no. 1, pp. 1–31, 2022.
  * [11] N. Bartoli and P. del Moral, Simulation et algorithmes stochastiques.  Cépaduès éditions, 2001.
  * [12] J. Gentle, Random Number Generation and Monte Carlo Methods.  Springer, 2004.
  * [13] C. Robert and G. Casella, Monte Carlo Statistical Methods (Springer Texts in Statistics).  Secaucus, NJ, USA: Springer-Verlag New York, Inc., 2005.
  * [14] S. Chib and E. Greenberg, “Understanding the Metropolis-Hastings algorithm,” The American Statistician, vol. 49, no. 4, pp. 327–335, 1995\.
  * [15] A. B. Owen, “Monte carlo theory, methods and examples,” 2013.
  * [16] J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods.  London: Methuen and Co., 1964.
  * [17] M. Castella, S. Rafi, P. Comon, and W. Pieczynski, “Separation of instantaneous mixtures of dependent sources using classical ICA methods,” EURASIP Journal on Advances in Signal Processing, p. n/c, Mar. 2013.
  * [18] H. Kahn and A. W. Marshall, “Methods of reducing sample size in Monte Carlo computations,” Journal of the Operations Research Society of America, vol. 1, no. 5, pp. 263–278, 1953.
  * [19] A. W. Marshall, “The use of multi-stage sampling schemes in Monte Carlo computations,” in Symposium on Monte Carlo Methods (M. Meyer, ed.), (New York), pp. 123–140, 1956.
  * [20] J. Geweke, “Bayesian inference in econometric models using Monte Carlo integration,” Econometrica, vol. 57, no. 6, pp. 1317–39, 1989.
  * [21] D. B. Rubin, “Using the SIR algorithm to simulate posterior distributions,” in Bayesian Statistics III (M. H. Bernardo, K. M. Degroot, D. V. Lindley, and A. F. M. Smith, eds.), Oxford: Oxford University Press, 1988.
  * [22] A. E. Gelfand and A. F. M. Smith, “Sampling based approaches to calculating marginal densities,” Journal of the American Statistical Association, vol. 85, no. 410, pp. 398–409, 1990.
  * [23] A. F. M. Smith and A. E. Gelfand, “Bayesian statistics without tears : a sampling-resampling perspective,” The American Statistician, vol. 46, no. 2, pp. 84–87, 1992.
  * [24] O. Cappé, É. Moulines, and T. Rydén, Inference in Hidden Markov Models.  Springer-Verlag, 2005.