# A New Dataset for Amateur Vocal Percussion Analysis

Alejandro Delgado  Roli Ltd.2, Glebe RoadLondonEngland, UKE8 4BD
<EMAIL_ADDRESS>,  SKoT McDonald  Roli Ltd.2, Glebe RoadLondonEngland, UKE8
4BD<EMAIL_ADDRESS>,  Ning Xu  Roli Ltd.2, Glebe RoadLondonEngland, UKE8 4BD
<EMAIL_ADDRESS>and  Mark Sandler  Queen Mary University of LondonMile End
RoadLondonEngland, UKE1 4NS<EMAIL_ADDRESS>
(2019)

###### Abstract.

The imitation of percussive instruments via the human voice is a natural way
for us to communicate rhythmic ideas and, for this reason, it attracts the
interest of music makers. Specifically, the automatic mapping of these vocal
imitations to their emulated instruments would allow creators to realistically
prototype rhythms in a faster way. The contribution of this study is two-fold.
Firstly, a new Amateur Vocal Percussion (AVP) dataset is introduced to
investigate how people with little or no experience in beatboxing approach the
task of vocal percussion. The end-goal of this analysis is that of helping
mapping algorithms to better generalise between subjects and achieve higher
performances. The dataset comprises a total of 9780 utterances recorded by 28
participants with fully annotated onsets and labels (kick drum, snare drum,
closed hi-hat and opened hi-hat). Lastly, we conducted baseline experiments on
audio onset detection with the recorded dataset, comparing the performance of
four state-of-the-art algorithms in a vocal percussion context.

dataset, beatbox, vocal, imitation, percussion

††copyright: acmcopyright††ccs: Applied computing Sound and music
computing††journalyear: 2019††conference: Audio Mostly; September 18–20, 2019;
Nottingham, United Kingdom††booktitle: Audio Mostly (AM’19), September 18–20,
2019, Nottingham, United Kingdom††price: 15.00††doi:
10.1145/3356590.3356844††isbn: 978-1-4503-7297-8/19/09

##  1\. Introduction

Music Information Retrieval (MIR), which uses Digital Signal Processing (DSP)
and machine learning techniques to analyse music recordings, has had a growing
influence on the music industry in recent decades. Results in tasks like music
genre recognition, chord estimation and source separation (Downie and Hao,
2018) bring optimism to the field in this regard. Some of the main sub-
disciplines in MIR are also receiving a significant amount of attention from
independent musicians. Automatic music transcription, for instance, enables
artists to learn and compose musical pieces more comfortably, while source
separation helps them practice more efficiently by deriving individual
instrumental tracks from the original audio mix.

In this context of innovation, a set of musical sounds worth exploring is that
of vocalised percussion. It includes vocal utterances that are articulated so
to communicate a rhythmic idea, usually by imitating the sound of percussive
instruments like those featured in a drum set. As such, these sounds and their
dynamics could be mapped to real drum samples from a sound library so to
create a realistic drum loop in seconds, making composers save time and effort
prototyping rhythms without actual music knowledge.

However, despite the possibilities that these rhythmic exploration tools
offer, they are seldom used today. This could be due to several factors,
ranging from the limited commercial spread of the already available
applications to the insufficient precision of their algorithms. As a response
to the current situation, we have recorded a dataset of vocal imitations of
percussion instruments, which could help researchers to shed light on the
problem and head towards a reliable and easy way of creating rhythmic
patterns.

A literature review of the work in vocalised percussion is presented in
section 2, alongside a list of the main datasets available in the public
domain. In Section 3, we introduce the Amateur Vocal Percussion (AVP) dataset,
laying out its main features and detailing its production process. In section
4, we perform baseline experiments on onset detection, where we compare the
performance of four well-established algorithms in the context of vocal
percussion. General conclusions are drawn in section 5.

##  2\. Previous Work

Between the recording of a vocal percussion performance and its mapping to a
drum pattern, there are usually two main steps: onset detection and utterance
classification. The former aims at localising the utterances in time, while
the latter deals with the association of these with percussion instruments.

There are two main frameworks in vocal percussion analysis for classification
purposes. Match and adapt tries to find the most similar spectrum to the query
one in a database of rhythmic performances (Nakano et al., 2004), while onset-
wise separation divides the audio file in utterances and analyse them
separately (Hazan, 2005). Both of them use MIR analytic routines, based on
frame-wise extraction and statistical aggregation of DSP descriptors (spectral
centroid, zero-crossing rate, mel-frequency cepstral coefficients…) and
machine learning techniques (support vector machines, decision trees…).
Studies are usually carried out focusing on three classes of imitated sounds:
kick drum, snare drum and hi-hat.

In most cases, a modest vocal percussion dataset is recorded for evaluation
purposes, which could be oriented to general rhythmic vocal percussion (Hazan,
2005) (Sinyor et al., 2005) (Ramires, 2017) (Gillet and Richard, 2005) or just
beatboxing (Kapur et al., 2004) (Stowell and Plumbley, 2010). Some of the
mentioned work target user-specific vocal percussion analysis rather than
general methods that could work for everyone, as vocal imitation styles seem
to change significantly from person to person (Hipke et al., 2014) (Ramires,
2017). Nevertheless, the idea of a universal model for vocal percussion
classification may also be plausible, as some studies have suggested that
humans, when given the task of imitating non-vocal environmental sounds,
naturally give realism to these imitations by focusing on certain
characteristic perceptual features of the original sound queries (Lemaitre et
al., 2016).

There are four widely adopted vocal percussion datasets in the public domain.
The first one is the Beatbox dataset by Stowell et al. (Stowell and Plumbley,
2010). It gathered experienced beatboxers to record 14 audio files (one per
participant) with a mean duration of 47 seconds that resulted in a total of
7460 annotated utterances. Both typical drum sounds and beatbox-specific ones
were labelled, being the audio files recorded in several environmental
conditions (different microphones, equipment, noise levels…). The authors also
discovered that a better classification performance could be achieved by
delaying the start of the first analysis frame around 23 milliseconds from the
actual onset, effectively discarding the transient information of the sounds.
The second one is the Live Vocalised Transcription (LVT) dataset by Ramires et
al. (Ramires, 2017). It was recorded by 20 participants (from experienced to
amateur beatboxers) using three microphones with different noise levels and
sound qualities. Two files per participant are provided, one of them
containing the imitation of a simple drum loop and the other a free rhythmic
improvisation. A software plugin for drum transcription by vocal percussion
was presented afterward (Ramires et al., 2018). The third dataset by Mehrabi
et al. (Mehrabi et al., 2018) is aimed at both transcribing the instruments
being imitated and recognising their specific model (e.g. type of snare). In
it, 14 participants with experience in music production are asked to imitate
30 percussion sounds, including different kick drums, snare drums, hi-hats,
toms and cymbals. Finally, the fourth dataset is the Vocal Imitation Set by
Kim et al. (Kim et al., 2018). This is directed to vocal imitation of sound
events in general, containing a total of 11242 crowd-sourced imitations of 302
different classes. It features imitations of several percussion instruments
from approximately twenty participants, including opened hi-hat, snare drum
and kick drum.

##  3\. THE AVP DATASET

As audio recording quality is becoming better and more accessible with time,
aspiring musicians are also moving away from the recording studios to a new
workspace consisting of just a relatively quiet room and a laptop. Enabling
independent artists to prototype rhythmic ideas without musical knowledge and
in an immediate way would make the creative workflow in this setting more
natural, quick and spontaneous. This is the main reason the Amateur Vocal
Percussion (AVP) dataset was created.

All audio files and annotations can be found in the link
https://doi.org/10.5281/zenodo.3250230

###  3.1. DESCRIPTION

The AVP dataset is characterised by the following attributes:

  * •

28 participants

  * •

A total of 9780 utterances in 280 audio files

  * •

Annotated onsets and labels (kick, snare, closed hi-hat and opened hi-hat)

  * •

Recorded with one microphone (MacBook Pro’s built-in mic)

  * •

2 modalities: personal imitations and fixed imitations

  * •

5 files for each modality each (four files with utterances of the same class
and one with an improvisation)

In contrast with the previously cited datasets, the AVP focuses exclusively on
people with little or no experience in beatboxing, making it suitable for the
goal mentioned before as well as functioning as an extension of past datasets.
It is also directed to query-by-example applications, allowing algorithms to
be trained on several isolated example utterances (around twenty-five per
label) before they are tested on the improvisation files.

The AVP is the largest vocal percussion dataset both in number of participants
and number of utterances. It also incorporates a second subset of fixed drum
imitations, i.e., four specific sounds given to all participants to imitate
each one of the four instruments. This subset was recorded in order to
investigate how practical a fixed technique of performing vocal percussion
could be, as deriving accurate algorithms to classify the personal subset as a
whole is relatively challenging. Phoneme recognition systems, for instance,
are likely to help in this particular endeavour; though, in comparison with
regular speech, different articulations and expressive dynamics are likely to
be encountered in rhythmic vocal percussion.

Table 1. Dataset content summary (utterances) Instrument - ‘label’ | Personal | Fixed | Improvisation
---|---|---|---
Kick Drum - ‘kd’ | 799 | 818 | 1201
Snare Drum - ‘sd’ | 813 | 839 | 811
Closed Hi-Hat - ‘hhc’ | 799 | 833 | 673
Opened Hi-Hat - ‘hho’ | 816 | 830 | 548

The main weakness of the AVP dataset is that, due to the lack of available
resources at the time, it was recorded without using any extra microphones.
However, when it comes to the dataset’s goal, this particularity does not
affect it in a critical way, as it is directed to drum pattern query by vocal
percussion in an indoor setting. Also, the realism of the recording context,
featuring quiet environmental noises as well, could be proven helpful for the
algorithms to perform well in non-controllable external conditions.

The total number of vocal percussion utterances belonging to each of the four
classes are gathered in table 1.

###  3.2. RECORDING METHODOLOGY

The materials used to record the dataset were a MacBook Pro laptop, GarageBand
software and a closed room of approximately 40 m3. The experiment took an
average of 15 minutes per participant.

The first step of the process was a small survey about the role that music
played in the participant’s life. This informal one-minute questionnaire was
oriented towards relaxing the participants and get them ready for the task.
Right after this, a standard loop featuring kick and snare was presented and
the participant was asked to both reproduce it vocally and write down the
onomatopoeias of the sounds he/she used in a notebook page. This performance
was intended as a first contact with the process and it did not get recorded,
while the annotation of the onomatopoeias was done to facilitate the recalling
of imitations and to better stick to them. As it turned out to be hard for
participants to imitate complex beats with hi-hat included in the preliminary
tests, five isolated utterances of both closed hi-hat and open hi-hat sounds
were presented instead. The participants decided their imitations and wrote
down their onomatopoeias in the same way as with the kick drum and snare drum.

Once the participants were familiarised with the task, they were asked to sit
naturally in front of the computer, as they would usually do, and the
recording of the dataset started. There were two recording modalities taking
place: personal and fixed. For the personal modality, participants used their
own vocal imitations to record around twenty-five utterances of kick drum,
snare drum, closed hi-hat and opened hi-hat sounds in four separated audio
files. The utterances followed a simple rhythmic loop (one crotchet and two
quavers) and a 90-BPM metronome track was provided through headphones to the
participants in order to record them. An audio file of improvised vocal
percussion featuring all or most instruments was recorded afterward. For the
fixed modality, the procedure above was repeated once more, but now four
specific sounds were given to the subjects. These fixed imitations were based
on speech syllables so to feel natural for participants to reproduce, and
their timbral characteristics were intended to mimic the imitated percussion
instruments. A “pm” syllable would correspond to the kick drum, “ta” to the
snare drum, “ti” to the closed hi-hat and a “chi” to the opened hi-hat. The
articulation of these utterances, as they were performed in a percussive
setting, generally resulted in brighter transient signals and more inharmonic
steady state signals compared to usual speech utterances.

###  3.3. CLEANING AND ANNOTATION

Once the raw audio files were recorded, two post-processing stages took place
to prepare them for analysis: trim silent regions at the beginning and the end
of each file to reduce their size and remove in-between passages where the
participants made accidental mistakes or notable pauses. The annotation of the
files was manually carried out by the author right after this cleaning
process, using Sonic Visualiser (Cannam et al., 2010) to write down both onset
locations and class labels. As illustrated in table 1, the tag ‘kd’ was used
for kick drum, ‘sd’ for snare drum, ‘hhc’ for closed hi-hat and ‘hho’ for
opened hi-hat.

The choice of the exact starting point of a sound event is generally
considered to be dependent on the task at hand and, thus, a matter of
convenience. For instance, if the goal is preparing an utterance with a long
transient to be used as a sound triggered by a MIDI message, the onset would
sometimes be preferably placed near the point of maximum energy in the signal.
In our case, we decided to place the onsets at the very beginning of the
utterance, where the percussive transient starts to build the sound. This is
due to the fact that our primary goal is to classify the utterance, and its
transient region could be informative as well.

An important note to make regarding the annotation process is that, in a few
occasions, two utterances were very close to each other or even appeared to
merge in one sound. A joint approach of waveform visualisation, spectrogram
visualisation and listening at quarter speed was employed to solve the
ambiguity in these cases. An example of how these specific annotations were
approached is illustrated in figure 1.

Figure 1. Annotation of two onsets, corresponding to two /s/ phonemes. Figure
1a) displays the utterances with their onsets marked in red. The waveform is
plotted in blue and the magnitude spectrogram in green. Figures 1b) and 1c)
are close captions of the first and the second onset respectively.

###  3.4. FURTHER OBSERVATIONS

A list of observations worth commenting on was made while recording and
listening back the audio files in the AVP dataset. Some of them are the
following:

  * •

A small set of recorded utterances exhibit a form of double percussion (such
as “suc” or “brrr”), making them challenging for the onset detector to output
the events as a whole without splitting them.

  * •

Various participants got mistaken in few occasions when recording the
improvisation files, using other speech-like phonemes like “cha” instead of
“ta” when imitating the snare drum or “dm” instead of “pm” when imitating the
kick drum. These passages were omitted in the final version of the dataset.

  * •

Several participants reported that the given fixed sounds, despite their
resemblance with the original drum sounds, felt unnatural for them to perform.

  * •

Participants generally improvised non-complex and predictable loops, which
could make the classification routines benefit from a rhythmic pattern
analyser.

Finally, there have been three participants whose utterances within the
personal dataset were either not consistent with each other, unintelligible or
practically indistinguishable from the rest. The audio files and annotations
pertaining to these cases are stored in the “Discarded” folder. Despite this,
they could still be used for classification purposes as long as the
improvisation files, where the ambiguities occur, are excluded.

##  4\. EXPERIMENT ON ONSET DETECTION

The problem of onset detection is a recurrent one in MIR literature. Most of
its subfields, in one way or another, need to face it at some point when
trying to make music analysis fully automatic. Vocal percussion resembles
speech when it comes to the articulation of utterances, making the job
slightly more challenging than with regular percussion, which exhibits short
and well-defined attacks. All vocal percussion utterances are usually composed
by one fricative phoneme that is sometimes followed by a vowel phoneme. This
is illustrated in figure 2.

Figure 2. Piece of waveform pertaining to a single vocal percussion utterance.
One can appreciate the plosive phoneme /t/ and the vowel phoneme /a/ composing
the “ta” syllable as distinct components in the audio file.

Here we describe a baseline study on onset detection for vocal percussion
using the AVP dataset. It aims at localising the time instances when the
utterances begin up to a tolerance margin of 50 ms (Downie and Hao, 2018).
This section is merely intended to be a brief analysis of the task, provoking
further discussions and laying down foundations for future work.

###  4.1. MODELS

In this section, we present the onset detection methods that will be evaluated
in this study. We will first do a quick overview of the spectral flux
descriptor and its relevance to the task. Then, we will take a look at four
state-of-the-art algorithms for audio onset detection. The first two,
Convolutional Neural Networks (CNN) and Recurrent Neural Networks RNN, are
based on deep machine learning methods while the other two, High-Frequency
Content (HFC) and Complex, are based on traditional signal processing
routines.

####  4.1.1. Spectral Flux (for onset refinement)

Despite this descriptor will not be used for onset detection but rather to
refine the locations of the detected onsets, all four algorithms investigated
here are either based on or inspired by it. For this reason, we will describe
it before we introduce the rest of the methods, considering it a theoretical
reference point.

The spectral flux is a frame-wise measure of how quickly the magnitude in each
frequency bin of the spectrum changes over time (Dixon, 2006). The sharper the
change in a particular region, the more probable a percussive onset is
occurring there. If $x$ is the audio waveform and $X$ is its Fourier
transform, the spectral flux is expressed in the following way:

(1) $SF(n)=\sum_{k=-\frac{N}{2}}^{\frac{N}{2}-1}H(|X(n,k)|-|X(n-1,k)|)$

where $k$ is a frequency bin of the $n$th frame and $H(x)=\frac{x+|x|}{2}$ is
called the half-wave rectifier function, which restricts the spectral flux to
only output positive changes in time, i.e., when the flux increases. The
higher a peak in the spectral flux function, the more likely an onset is
taking place in that location.

####  4.1.2. Convolutional Neural Network (CNN)

CNNs (Krizhevsky et al., 2012) are feed-forward artificial neural networks
composed of convolutional layers. The neurons in these layers compose a set of
small local filter kernels to analyse the input, creating several feature maps
as a result. A convolutional layer could be followed by a pooling layer, which
subsamples these resulting feature maps following a certain set of rules, or
by a fully connected layer to introduce non-linearities in the process. These
last layers are specially useful for classification purposes.

CNNs are especially popular in computer vision, and they achieve high
accuracies when detecting edges in images. The main idea behind their use for
audio onset detection is that they can effectively detect edges in the
magnitude spectrograms as well, which usually correspond to audio onsets. In
this way, the usual input for the CNN is composed by three 80-band mel
spectrograms with constant hop size but different window sizes so to have both
high frequency and time resolution. The details of the original network’s
architecture can be found in (Schluter and Böck, 2014). It was discovered in
the same study that the CNN, like spectral flux based techniques, computes
spectral differences over time and uses them to estimate the likelihood of an
onset happening in a certain time region. The implementation of the CNN-based
onset detection method, which this study follows, is contained in Madmom’s
CNNOnsetProcessor (Böck et al., 2016).

####  4.1.3. Recurrent Neural Network (RNN)

The functioning of an RNN, in contrast with CNNs, is dependent on the
information acquired from past inputs, allowing these networks to model
complex data sequences in a non-linear way. Long Short-Term Memory (LSTM)
units in RNNs (Hochreiter and Schmidhuber, 1997) allow the networks to better
handle dependencies in time and avoid the exploding/vanishing gradient
problem.

When spectrograms are fed to RNNs, these networks can detect and model changes
in time by contrasting observations from current frames with information from
previous ones. Their routines, then, would be analogous to those of the
spectral flux onset detector, with RNNs taking information from more past
frames into account. The original RNN model for onset detection can be found
in (Böck et al., 2012). This time, the authors use three Bark spectrograms
with constant hop size and different window sizes as inputs, and the resulting
model is also optimised to handle less data than other onset detection
algorithms. This makes it suitable for real-time onset detection, using
unidirectional (causal) RNNs without LSTM cells. The implementation of the
RNN-based onset detection method, which this study follows, is contained in
Madmom’s RNNOnsetProcessor (Böck et al., 2016). We use the non-real-time
version, with bidirectional RNNs.

####  4.1.4. High Frequency Content Method (HFC)

In general, most of the energy of a sound’s steady state part is located in
the lower frequencies, being the higher frequencies more represented in its
transient part. When a percussive onset occurs, the sudden increase in energy
tends to be particularly more pronounced at higher frequencies, indicating the
beginning of its transient.

The HFC method for onset detection (Masri, 1996) exploits this insight making
use of the high-frequency content descriptor, i.e., the summation of the bin
magnitudes of the spectrum multiplied by their own bin indices. This can be
formalised in the following way:

(2) $\mathrm{HFC}=\sum_{i=0}^{N-1}i|X(i)|$

for a STFT spectrum $X$ with bin indices $[0...n]$. The result, which is
linearly biased towards high frequencies, is then fed to a detection function
that returns the predicted onset location. This detection function is based on
both the high frequency energy flux between the two adjacent frames and the
normalised high frequency content of the current frame. The implementation
that this study follows is included in Aubio’s AubioOnset (Böck et al., 2016).
It has also been reported that the HFC method works especially well for vocal
percussion (Ramires, 2017).

####  4.1.5. Complex Method (Complex)

The phases of the partials in the steady state of a sound can be accurately
predicted from past frames, as their frequencies and amplitudes remain
constant in this region. In the case of its transient part, however, these
phases tend to evolve in a non-linear, unpredictable way. We can use this
particularity of the transient’s STFT phase spectrum to locate it within the
waveform.

The complex method for onset detection effectively combines both energy-based
and phase-based approaches simultaneously by detecting onsets in the complex
domain (Duxbury et al., 2003). The complex detection function is then pre-
processed with a weighted moving average and fed to a peak-picking algorithm
that outputs the predicted onset locations. The implementation that this study
follows is included in Aubio’s AubioOnset (Böck et al., 2016).

###  4.2. EVALUATION AND RESULTS

We carry out three different studies to retrieve the parameters that produce
the best performances for each algorithm in a vocal percussion context. Grid
search is performed for each parameter using ten linearly spaced values. We
test the algorithms on all the utterances in the AVP dataset and the accuracy
of each performance is measured using the F1-Score.

The first analysis aims at retrieving the best values for the method-specific
parameters. These parameters are the peak picking threshold of the onset
probability functions for all methods and the frame size for the HFC and
Complex ones. It has been found that the best results for each algorithm are
reached by defining a threshold of 0.55 for CNN, 0.30 for RNN, 0.8 for HFC and
0.7 for complex. A frame size of 11 ms gave the best results for both the HFC
and the Complex methods.

The second study explores how a minimum separation between onsets could
further benefit the algorithms’ performance. This means that in case two of
the predicted onsets happen to be separated in time by less than the value of
a certain threshold, the second onset will be automatically discarded. This
appears to be relevant to our case for two reasons. The first one is that two
vocal percussion utterances occur seldom or never simultaneously, i.e., we are
dealing with a monophonic transcription problem. Therefore, there will always
exist a minimum separation between onsets. The second reason is that these
algorithms can easily confuse the transient of a vowel phoneme with the onset
of a new utterance, lowering their precision score significantly. This
inconvenience could be avoided with the inclusion of a fixed minimum
separation between onsets, as the onset from the vowel phoneme could be
effectively discarded.

Results for this second study are displayed in figure 3. The main observation
is that both the CNN and RNN methods benefit from the inclusion of a minimum
onset separation, reaching the peak of their performance at the 90 ms mark,
whilst the HFC and the Complex methods show a detriment in performance when
this separation is applied. A reason for this phenomenon could be that both
the HFC and the Complex algorithms are based on distinctive features of
percussive transients, being them the high-frequency content and the
predictability of the phase spectrum respectively. For these methods, the
peaks in the onset probability function corresponding to plosive phonemes (see
figure 2) would be appreciably higher than the peaks relative to the vowel
phonemes. Hence, setting an appropriate value for the peak-picking threshold
would be enough to separate these types of phonemes, and a minimum onset
separation would only make the algorithm discard correct onsets in situations
when two real onsets are significantly close to each other. Deep learning
based algorithms, on the other hand, are optimised so to detect musical onsets
in general and do not explicitly rely on these characteristics of percussive
transients, which makes them more likely to detect the onsets of the vocal
phonemes as well. Thus, setting a minimum onset separation is a simple yet
effective way for them to ignore onsets coming from vowel phonemes, at the
expense of discarding correct onsets that are close in time.

Figure 3. Effect of different minimum separation values on the performance of
the four onset detection algorithms.

Finally, we attempt to refine these onsets in time using the spectral flux
descriptor. The goal here is to relocate each predicted onset in a
neighbouring point which is closer to the real onset. More technically, it
consists in minimising the mean absolute deviation of the predicted onsets
with respect to the real onsets without significantly affecting the
performance (F1-Score). We approach this task by centring an analysis window
in the predicted onset and then finding the maximum of the spectral flux
function in that region, where the real onset is more likely to be. We
optimised the length of this analysis window selecting the value that achieved
the best balance in the performance vs. precision trade-off. Specifically, we
set a maximum allowed drop in performance of 1% of its original F1-Score,
considering only the parameter values that did not make the algorithms exceed
that limit.

The final results are gathered in table 2. The HFC and the Complex methods
perform the best when it comes to accuracy and the onsets derived from deep
learning methods are the closest in time to the real onsets. Results also
present CNNs and RNNs as slow compared to HFC and Complex, which are 38 and 27
times faster respectively.

Method | F1-Score | Deviation | Duration
---|---|---|---
CNN | 0.89 |  9 ms | 267 s
RNN | 0.89 | 10 ms | 266 s
HFC | 0.94 | 17 ms |  7 s
Complex | 0.95 | 21 ms | 10 s
CNN-SF | 0.89 |  6 ms | -
RNN-SF | 0.88 | 9 ms | -
HFC-SF | 0.93 | 11 ms | -
Complex-SF | 0.94 | 11 ms | -
Table 2. Final results of the onset detection study for vocal percussion. The
’-SF’ suffix indicates the inclusion of spectral flux refinement in the
process. From left to right, it is displayed the F1-Score, the mean absolute
deviation of the predicted onsets with respect to the real ones, and the
duration of the analysis for all onsets using the same laptop.

##  5\. CONCLUSION

In this piece of work, we have presented the AVP dataset, comprising 9780
utterances of vocal percussion recorded by 28 participants. This aims at
improving algorithms for drum pattern query by vocal imitation so that
independent musicians can accelerate their creative routines by sketching
realistic percussion rhythms on the fly. Baseline experiments on utterance
onset detection were carried out in the final section. It was shown that
methods relying on specific DSP features outperformed deep learning based
techniques in the context of vocal percussion, laying the groundwork for novel
approaches to improve upon.

###### Acknowledgements.

$\begin{array}[]{l}\includegraphics[scale={0.1}]{EU_Low.jpg}\end{array}$ This
project has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Sk$\l{}$odowska-Curie grant agreement
No. 765068 A warm and special thank you to the people that participated in the
creation of this dataset.

## References

  * (1)
  * Böck et al. (2012) Sebastian Böck, Andreas Arzt, Florian Krebs, and Markus Schedl. 2012.  Online real-time onset detection with recurrent neural networks. In _Proceedings of the 15th International Conference on Digital Audio Effects DAFx-12, York, UK_.
  * Böck et al. (2016) Sebastian Böck, Filip Korzeniowski, Jan Schlüter, Florian Krebs, and Gerhard Widmer. 2016\.  madmom: a new Python Audio and Music Signal Processing Library.  _arXiv:1605.07008_ (2016).  https://doi.org/10.1145/2964284.2973795
  * Cannam et al. (2010) Chris Cannam, Christian Landone, and Mark Sandler. 2010\.  Sonic visualiser: An open source application for viewing, analysing, and annotating music audio files. In _Proceedings of the 18th ACM international conference on Multimedia_. ACM, 1467–1468.  https://doi.org/10.1145/1873951.1874248
  * Dixon (2006) Simon Dixon. 2006\.  Onset Detection Revisited. In _Proceedings of the 9th International Conference on Digital Audio Effects DAFx-06, Montreal, Canada_.
  * Downie and Hao (2018) J Stephen Downie and Yun Hao. 2018.  MIREX 2018 Evaluation Results.  (2018).
  * Duxbury et al. (2003) Chris Duxbury, Juan Pablo Bello, Mike Davies, and Mark Sandler. 2003\.  Complex domain onset detection for musical signals. In _Proceedings of the 6th International Conference on Digital Audio Effects DAFx-03, London, UK_.
  * Gillet and Richard (2005) Olivier Gillet and Gaël Richard. 2005.  Drum Loops Retrieval from Spoken Queries.  _Journal of Intelligent Information Systems_ 24, 2-3 (March 2005).  https://doi.org/10.1007/s10844-005-0321-9
  * Hazan (2005) Amaury Hazan. 2005\.  Towards automatic transcription of expressive oral percussive performances. In _Proceedings of the 10th international conference on Intelligent user interfaces - IUI ’05_. ACM Press, San Diego, California, USA.  https://doi.org/10.1145/1040830.1040904
  * Hipke et al. (2014) Kyle Hipke, Michael Toomim, Rebecca Fiebrink, and James Fogarty. 2014.  BeatBox: end-user interactive definition and training of recognizers for percussive vocalizations. In _Proceedings of the 2014 International Working Conference on Advanced Visual Interfaces - AVI ’14_. ACM Press, Como, Italy.  https://doi.org/10.1145/2598153.2598189
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  _Neural computation_ 9, 8 (1997), 1735–1780.  https://doi.org/10.1162/neco.1997.9.8.1735
  * Kapur et al. (2004) Ajay Kapur, Manj Benning, and George Tzanetakis. 2004.  Query-By-Beat-Boxing: Music Retrieval for the DJ. In _5th International Conference on Music Information Retrieval (ISMIR), Barcelona, Spain_.
  * Kim et al. (2018) Bongjun Kim, Madhav Ghei, Bryan Pardo, and Zhiyao Duan. 2018\.  Vocal Imitation Set: a dataset of vocally imitated sound events using the AudioSet ontology.  _Proceedings of the Detection and Classification of Acoustic Scenes and Events 2018 Workshop (DCASE2018)_ (Nov. 2018).  https://doi.org/10.5281/zenodo.1340763
  * Krizhevsky et al. (2012) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton. 2012\.  Imagenet classification with deep convolutional neural networks. In _Advances in neural information processing systems_. 1097–1105.  https://doi.org/10.1145/3065386
  * Lemaitre et al. (2016) Guillaume Lemaitre, Olivier Houix, Frédéric Voisin, Nicolas Misdariis, and Patrick Susini. 2016\.  Vocal Imitations of Non-Vocal Sounds.  _PLOS ONE_ 11, 12 (Dec. 2016).  https://doi.org/10.1371/journal.pone.0168167
  * Masri (1996) Paul Masri. 1996\.  Computer Modelling of Sound for Transformation and Synthesis of Musical Signals.  _University of Bristol_ (1996).
  * Mehrabi et al. (2018) Adib Mehrabi, Keunwoo Choi, Simon Dixon, and Mark Sandler. 2018\.  Similarity measures for vocal-based drum sample retrieval using deep convolutional auto-encoders.  _arXiv:1802.05178_ (Feb. 2018).  https://doi.org/10.1109/icassp.2018.8461566
  * Nakano et al. (2004) Tomoyasu Nakano, Jun Ogata, Masataka Goto, and Yuzuru Hiraga. 2004.  A Drum Pattern Retrieval Method by Voice Percussion. In _5th International Conference on Music Information Retrieval (ISMIR), Barcelona, Spain_.
  * Ramires et al. (2018) António Ramires, Rui Penha, and Matthew E. P. Davies. 2018\.  User Specific Adaptation in Automatic Transcription of Vocalised Percussion.  _arXiv:1811.02406_ (2018).
  * Ramires (2017) António Filipe Santana Ramires. 2017\.  Automatic Transcription of Drums and Vocalised percussion.  _Universidade do Porto_ (2017).
  * Schluter and Böck (2014) Jan Schluter and Sebastian Böck. 2014.  Improved Musical Onset Detection with Convolutional Neural Networks.  _Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ (2014).  https://doi.org/10.1109/ICASSP.2014.6854953
  * Sinyor et al. (2005) Elliot Sinyor, Cory Mckay Rebecca, Daniel Mcennis, and Ichiro Fujinaga. 2005.  Beatbox classification using ACE. In _Proceedings of the International Conference on Music Information Retrieval_.
  * Stowell and Plumbley (2010) Dan Stowell and Mark D. Plumbley. 2010.  Delayed decision-making in real-time beatbox percussion classification.  _Journal of New Music Research_ (2010).  https://doi.org/10.1080/09298215.2010.512979


| $SB1+p(m_{1}|\theta)$ | $57.142_{50.329}^{61.911}$ | $1990.277_{1989.567}^{1990.943}$ | $0.395_{0.369}^{0.442}$ | $0.315_{0.285}^{0.336}$ | $59.717_{42.404}^{70.114}$ | $261.708_{257.161}^{267.265}$ | $30.198_{16.674}^{36.155}$ | $-22.869_{-23.137}^{-22.710}$ | $18.463_{16.757}^{19.099}$ | $20.127_{17.627}^{29.597}$ | $0.955_{0.926}^{0.985}$ | $0.591_{0.506}^{0.992}$
| $SB1+p(\varpi)+p(m_{1}|\theta)$ | $53.035_{50.393}^{54.951}$ | $1989.717_{1989.475}^{1990.167}$ | $0.421_{0.405}^{0.448}$ | $0.293_{0.284}^{0.301}$ | $49.694_{41.413}^{56.121}$ | $261.708_{256.664}^{267.486}$ | $19.877_{15.208}^{23.874}$ | $-22.957_{-23.153}^{-22.814}$ | $16.843_{16.488}^{17.234}$ | $29.109_{26.896}^{30.299}$ | $0.958_{0.934}^{0.990}$ | $0.962_{0.859}^{1.000}$
115126 | Orb6 | 6.321 | 2012.301 | 0.173 | 0.189 | 28.3 | 314.9 | 44.5 | — | — | — | — | —
| SB9 | 6.292 | 2006.032 | 0.1620 | — | 212.13 | — | — | 9.9 | — | — | — | —
| $SB1$ | $6.325_{6.321}^{6.328}$ | $1980.738_{1980.724}^{1980.755}$ | $0.160_{0.158}^{0.163}$ | $0.191_{0.189}^{0.193}$ | $212.995_{212.552}^{213.564}$ | $341.293_{340.775}^{341.807}$ | $49.712_{48.803}^{50.741}$ | $-1.613_{-1.622}^{-1.594}$ | — | $8.695_{8.503}^{8.856}$ | — | —
| $SB1+p(\varpi)$ | $6.325_{6.322}^{6.328}$ | $1980.737_{1980.724}^{1980.755}$ | $0.161_{0.158}^{0.163}$ | $0.191_{0.189}^{0.193}$ | $213.001_{212.543}^{213.586}$ | $341.154_{340.807}^{341.829}$ | $49.900_{48.771}^{50.711}$ | $-1.606_{-1.622}^{-1.595}$ | $44.749_{43.846}^{46.018}$ | $8.658_{8.512}^{8.859}$ | $1.188_{1.055}^{1.282}$ | $0.633_{0.609}^{0.676}$
| $SB1+p(m_{1}|\theta)$ | $6.325_{6.322}^{6.328}$ | $1980.738_{1980.724}^{1980.754}$ | $0.160_{0.158}^{0.163}$ | $0.190_{0.189}^{0.192}$ | $213.114_{212.569}^{213.577}$ | $341.212_{340.773}^{341.790}$ | $49.573_{48.802}^{50.741}$ | $-1.610_{-1.623}^{-1.595}$ | $44.299_{43.910}^{45.266}$ | $8.726_{8.503}^{8.853}$ | $1.217_{1.154}^{1.253}$ | $0.630_{0.615}^{0.649}$
| $SB1+p(\varpi)+p(m_{1}|\theta)$ | $6.326_{6.322}^{6.328}$ | $1980.739_{1980.724}^{1980.755}$ | $0.160_{0.158}^{0.163}$ | $0.191_{0.189}^{0.193}$ | $213.227_{212.546}^{213.570}$ | $341.296_{340.764}^{341.813}$ | $49.887_{48.835}^{50.741}$ | $-1.611_{-1.623}^{-1.596}$ | $44.758_{44.095}^{45.270}$ | $8.647_{8.507}^{8.837}$ | $1.195_{1.151}^{1.241}$ | $0.631_{0.615}^{0.648}$

###  4.4 Concluding remarks

The experimental results presented in this section show the adequacy and
usefulness of the proposed Bayesian inference methodology to characterize and
visualize the posterior uncertainty of $SB1$ visual-spectroscopic binary
systems. Interesting results are obtained in some of the binary systems
evaluated. For example, HIP 3504 shows an extremely high uncertainty on its
orbit and RV curve, due to the few available observations. However, the
incorporation the priors allows us to radically reduce the (posterior)
uncertainty of the estimates, reaching completely different solutions to the
ones obtained without the prior information. This is reflected in a tighter
orbit and in a RV curve with a lower period in Fig.14. The drastic change in
period, reflects also in a large difference in the predicted systemic velocity
$V_{0}$ (see Table 6).

The results obtained for the system HIP 99675 show that the effect of the
prior incorporation is mostly concentrated in the trio of soft-identifiable
parameters $\varpi$, $m_{1}$ and $q$, affecting only slightly the other
orbital parameters. This behavior express the robustness of the estimation
when abundant and precise observations are available. On the soft-identifiable
parameters, we observe that the posterior distributions of the
$SB1+p(m_{1}|\theta)$ and the $SB1+p(\varpi)+p(m_{1}|\theta)$ cases are
similar, but usually quite different to the case $SB1+p(\varpi)$. This
behavior can be interpreted as the prior information of the primary object
mass $p(m_{1}|\theta)$ being more reliable than the prior on the parallax
$p(\varpi)$ for this system. Finally, the results obtained for the system HIP
109951 show the usefulness of providing a good characterization of the
uncertainty. In this case, we observe that the uncertainty in the half portion
of the orbit with no observations is considerably higher than the other half,
but most importantly, the uncertainty of the RV curve visibly increases with
time. This behavior is attributed to the large uncertainty of the period of
the system, making the trajectories of the RV get out of phase with time.

On the prior information incorporated, we observe that the posterior
distribution of the cases $SB1+p(\varpi)$ and $p(m_{1}|\theta)$ are very
similar to the one obtained in the $SB1$ case. However the mixed case
$SB1+p(\varpi)+p(m_{1}|\theta)$ shows a visible uncertainty reduction in all
the orbital parameters with a very similar MAP value. Regarding the posterior
distribution of the mass ratio $q$ for HIP 109951, we observe that the cases
$SB1+p(\varpi)$ and $p(m_{1}|\theta)$ provide a poor estimation, reflected in
a very high uncertainty. However, the mixed case
$SB1+p(\varpi)+p(m_{1}|\theta)$ exhibits a much more constrained posterior
distribution than the other cases, considerably differing in their MAP
estimates. It is important to note that, as the posterior distribution on the
mass ratio $q$ of the $SB1+p(\varpi)$ and $p(m_{1}|\theta)$ cases present a
very high uncertainty (with the shape of an almost uniform distribution), its
corresponding MAP estimates are not very reliable, but the mixed case
$SB1+p(\varpi)+p(m_{1}|\theta)$ does. This last result shows the relevance of
incorporating both sources of prior information to obtain more robust
estimates, and it explains the large variation in the predicted $q$ in the
three cases with priors. Indeed, a similar behavior was already seen in the
case of the benchmark HIP 89000 (see the bottom right panel on Figure 7).

By looking at Table 6 we can see that, besides HIP 109951 discussed above, two
other cases have large variation in the predicted $q$ depending on the prior
used, namely HIP 65982 and 79101, and we discuss them in turn. For HIP 65982
the mixed case induces a large bias in some orbital parameters and a
considerable increase in the uncertainty for $q$ (see the corresponding PDFs
in our web site). For these reasons, in this case, one should probably favor
the $SB1+p(\varpi)$ solution (which is similar to the $SB1++p(m_{1}|\theta)$
case), and a rather small value of $q\sim 0.3$

For HIP 79101, the most notable ”feature” is a bi-modal distribution of the
PDF (see our web site) for the auxiliary parameter $f/\pi$ in the $SB1+p(\pi)$
and $SB1+p(m1|\theta)$ scenarios (the MAP value is the left peak) which is
”resolved” by the mixed priors case (uni-modal), in which the MAP value
coincides now with the right peak (which shows also that the value is not
biased), all of which is also reflected in the distribution of $q$ values.
Therefore, in this case it would seem that the mixed case is preferred since
it has less uncertainty in $q$ (for the other scenarios the distribution of
$q$ is very broad), and it is not biased. A final note regarding this system:
Skiff´s catalogue of spectral types (Skiff, 2014) gives a broad range of
possibilities between B8V to A0V, B9IV, B9III and even A0II, from fifteen
different sources, our adopted value being a B9V. Considering its reported
$V=4.27$ on SIMBAD, and its adopted trigonometric parallax in Table 5, the
predicted $M_{V}=+0.14$, this corresponds indeed to a B8V-B9V from Abushattal
et al. (2020). On the other hand, a B9IV should have $M_{V}=-0.05$, a B9II
should have $M_{V}=-0.50$, and an A0II should have $M_{V}=-3.4$, all of them
far from our observed value. We thus conclude that our adopted spectral type
is indeed reasonable.

##  5 Conclusions and final comments

The Bayesian methodology for the inference of orbital parameters in $SB1$
binary systems with a visual orbit proposed in this paper allows us to provide
a computationally-efficient, robust, and precise estimation of the
corresponding joint posterior distributions of these parameters. This
inference is implemented through the No-U-Turn sampler MCMC algorithm, which
allows the incorporation of prior information of the stellar system to
constrain the inference in scenarios with imprecise or missing data. The
flexibility of this sampling scheme to add prior information is very useful
for an estimation of individual component masses in $SB1$ visual-spectroscopic
binaries.

An exhaustive experimental analysis has been carried out for the validation of
the proposed methodology. We study the quality of the inference by comparing
the estimated posterior distribution of well-studied $SB2$ visual-
spectroscopic binaries with their $SB1$ visual spectroscopic counterparts by
omitting the RV observations of the companion object. Our results show a
negligible difference between the estimated posterior distributions of the
orbital parameters (and their uncertainties) when compared to the benchmark
(full observation) case, in which the RV of both components are considered.
This is a very promising result, showing that partial observations ($SB1$
case) offers a good estimation performance.

Our empirical results indicate that the incorporation of prior information of
the system (through the trigonometric parallax and the mass of the primary
object) allows an estimation of the mass ratio of the system (and hence the
individual component masses) with good precision. The incorporation of the
prior distributions makes those parameters identifiable, where the derived
estimations — position, dispersion, and shape of the marginal posterior
distribution — strongly depend on the prior chosen. This prior knowledge has
also an influence on the estimated posterior distribution of the other orbital
parameters. The impact and relevance of incorporating priors on the inference
of previous-identifiable orbital parameters (that are already identifiable
without the incorporation of the priors) depends on the abundance, precision,
and orbital coverage of the available observations. In particular, it is
observed that if the system is precisely determined, the impact of the prior
on the estimation of those set of parameters is negligible, as expected.

Our numerical results show that the lowest estimation error (from the optimal
MAP estimator of all the systems analyzed) on the system’s mass ratio, with
respect to the full-information scenario (with both RV observations), was
achieved by the mixed case that incorporates prior information of
trigonometric parallax and a mass for the primary object simultaneously
($4.92\%$), while the highest error was obtained by incorporating a prior only
on the mass of the primary alone ($7.44\%$), achieving a percentage error
lower than $8\%$. It is shown that the closest marginal posterior distribution
to the full-information scenario —in the KLD sense— was achieved by the
incorporation of a prior on the system parallax alone. The lowest similitude
was obtained by the incorporation of a prior on the primary object mass alone,
which is attributed to the fact that the parallax prior information is
probably more reliable and precise than the prior on the mass of the primary
star, at least for our benchmark systems. Overall, the incorporation of both
priors was the most beneficial to the accuracy of the MAP estimates, where,
when more information is provided, better estimation can be obtained. Large
differences in the posterior distribution of parameters depending on the prior
imposed could signal that one of them is biased, and should thus be taken with
caution.

Taking advantage of the flexibility and richness of our sample-based
methodology, the differences between the estimated posterior distributions in
all the studied cases was also analyzed in the corresponding observations
space. This novel analysis provides a better understanding of the effect of
the different sources of information on the shape and uncertainty of the orbit
and RV curves of the stellar systems. Finally, we applied the proposed
Bayesian framework to twelve previously unstudied $SB1$s with astrometric
data, providing a complete analysis of the obtained results.

The present work addresses the classical orbital parameters estimation of
binary stellar systems through a Bayesian perspective, emphasizing the
importance of providing not only an estimation, but also a complete
characterization of the posterior distribution of the orbital parameters. This
approach allows us to provide an uncertainty quantification of the inference
process (as many classical optimization-fitting methods roughly provide), but
it also allows us to visualize the uncertainty of the orbit itself in the
observation space. This last dimension in our analysis is fundamental to
decide what (and when) future measurements are required to improve the
precision of the estimation. In this way we show how an in-depth statistical
analysis can provide important insights from an observational planing
perspective.

We are very grateful to the referee for his/her careful reading of our
manuscript and, in particular, for pointing out the relevance of verifying the
reliability of the adopted spectral types for the priors. MV and RAM
acknowledge support from FONDECYT/ANID grant No. 1190038. MV and JS
acknowledge support from FONDECYT/ANID grant No. 1210315. MO & JS acknowledge
support from the Advanced Center for Electrical and Electronic Engineering,
AC3E, Basal Project FB0008, ANID. We are very grateful for the continuous
support of the Chilean National Time Allocation Committee under programs
CN2018A-1, CN2019A-2, CN2019B-13, CN2020A-19, CN2020B-10 and CN2021B-17. This
research has made use of the Washington Double Star Catalog maintained at the
U.S. Naval Observatory and of the SIMBAD database, operated at CDS,
Strasbourg, France. This work has made use of data from the European Space
Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the
Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

##  Appendix A Orbital model

###  A.1 Visual binary systems

A visual binary (or visual double system) corresponds to a gravitationally
bound binary system in which the relative position of both components are
observable. The positional or astrometric observations in binary systems
measures the relative position of the fainter or companion object with respect
to the brighter or primary object. Depending on the technique used to obtain
the astrometric observations, these can be classified roughly into
micrometric, photographic, or interferometric positional measurements.

The solution of the differential equation of the motion is described by
Kepler´s laws (Equation (A2)), assuming both objects behave as point masses
(detached binaries), where in the case of binary stars, it corresponds to an
elliptical orbit where the primary star is in the focus and the area swept by
the radius vector is constant per unit time. This elliptical orbit, denoted as
real (relative) orbit, is characterized by four orbital parameters:

  * •

Period ($P$): The revolution period in years.

  * •

Time of periastron passage ($T$): One epoch of passage through periastron
(minimum true distance between the components) in years and fraction of a
year.

  * •

Semi-major axis ($a$): The major semiaxis of the elliptical true orbit in
seconds of arc.

  * •

Eccentricity ($e$): The numerical eccentricity.

The astrometric observations are position measurements of the projection of
the real orbit in the plane of the sky relative to the observer (plane of
reference), denoted as apparent orbit. Three additional parameters are
necessary to project the real orbit into the apparent orbit:

  * •

Longitude of the ascending node ($\Omega$): The position angle from a
reference direction to the ascending node777Point where the real orbit of the
object passes through the plane of reference. in the plane of reference
(ranging from 0° to 360°).

  * •

Argument of periapsis ($\omega$): The angle from the node to the periastron in
the real orbit, following the direction of motion (ranging from 0° to 360°).

  * •

Inclination ($i$): The angle between the plane of projection and that of the
true orbit (ranging from 0° to 180°).

It is worth pointing out that two values for the longitude of the ascending
node ($\Omega$ and $\Omega+\varpi$) results in identical apparent orbits.
Therefore, the ascending node cannot be identified by positional observations.
By convention in astronomy, if the ascending node is undetermined, the value
of $\Omega$ is placed in the first two quadrants, i.e., from 0o to 180o.

On the specifics, the position on the apparent orbit $(\rho,\theta)$ at a
certain time $t$ (the ephemerides formulae) involves the determination of the
position in the real orbit and its projection to the apparent orbit. The
position on the real orbit involves the determination of the three orbital
anomalies: the true anomaly $\nu(t)$, the eccentric anomaly $E(t)$ and the
mean anomaly $M(t)$, in terms of the orbital parameters $P$, $T$, $a$ and $e$.

The true anomaly $\nu(t)$ is defined as the angle between the periapsis and
the current position of the companion object in the orbit, as seen from its
main focus (the position of the primary object). This angular parameter can be
determined by the following geometrical identity:

$\tan{\frac{\nu(t)}{2}}=\sqrt{\frac{1+e}{1-e}}\tan{\frac{E(t)}{2}}.$ (A1)

where $E(t)$ is the eccentric anomaly. The eccentric anomaly is defined as the
angle between the periapsis and the intersection of a perpendicular line to
the semi-major axis of the orbit and the position of the companion object in
the orbit, as seen from its central point. This angular parameter can be
determined by the numerical resolution of Kepler´s equation888The Kepler
equation is commonly resolved using the Newton-Raphson method Ypma (1995).:

$M(t)=\frac{2\pi(t-T)}{P}=E(t)-e\sin E(t),$ (A2)

where $M(t)$ is the mean anomaly of the orbit. The mean anomaly represents the
angular movement of the companion object in the orbit (similar to the true
anomaly) but at a uniform rate. The uniform rate of movement is represent by a
circle circumscribed to the orbit. Therefore, the mean anomaly is defined as
the angle between the periapsis and the point in the circle circumscribed to
the orbit, as seen from its central point, that covers the same area per unit
time as the true anomaly. This angular parameter corresponds to the revolution
period of the companion object in the orbit, i.e., $M(t)=2\pi(t-T)/P$.

Finally, the position on the real orbit is projected to the apparent orbit
through the angular parameters $\omega$, $\Omega$ and $i$:

$\begin{split}\tan(\theta(t)-\Omega)&=\tan(\nu(t)+\omega)\cos(i),\\\
\rho(t)&=r(t)\cos(\nu(t)+\omega)\sec(\theta(t)+\Omega),\end{split}$ (A3)

with $r(t)=a(1-e^{2})/(1+e\cos(\nu(t)))$ is the radius vector.

The procedure to compute the position on the apparent orbit in rectangular
coordinates ($X,Y$) at a certain time $t$ involves determining the normalized
rectangular coordinates in the true orbit $x,y$:

$\begin{split}x(t)&=\cos{E(t)}-e,\\\
y(t)&=\sqrt{1-e^{2}}\sin{E(t)},\end{split}$ (A4)

with $E(t)$ the eccentric anomaly determined in (A2). Therefore, the position
on the true orbit is computed by a ponderation of the normalized coordinates:

$\begin{split}X(t)&=Ax(t)+Fy(t),\\\ Y(t)&=Bx(t)+Gy(t),\end{split}$ (A5)

with $A,B,F$ and $G$ are the so-called the Thiele-Innes elements, defined as:

$\begin{split}A&=a(\cos{\omega}\cos{\Omega}-\sin{\omega}\sin{\Omega}\cos{i}),\\\
B&=a(\cos{\omega}\sin{\Omega}+\sin{\omega}\cos{\Omega}\cos{i}),\\\
F&=a(-\sin{\omega}\cos{\Omega}-\cos{\omega}\sin{\Omega}\cos{i}),\\\
G&=a(-\sin{\omega}\sin{\Omega}+\cos{\omega}\cos{\Omega}\cos{i}).\end{split}$
(A6)

The terms $(A/a,B/a)$ and $(F/a,G/a)$ are interpreted as the direction cosines
of the major and minor axis, respectively, of the orbit in the rectangular
coordinate system formed by the tangential plane and the North direction (more
specifically, $X$ (or $x$) and $Y$ (or $y$) point in the North and East
directions respectively). The Thiele-Innes elements form a one-to-one
correspondence with the elements $a,\Omega,\omega,i$.

###  A.2 Spectroscopic binary systems

A spectroscopic binary system corresponds to a binary system in which the
spectral lines of the light emitted by its components are observable. The
movement of the stars in the orbit produces variations on the spectral lines
observed as a consequence of the Doppler’s effect: Blue or red-shifted lines
are measured when the stars moves towards or away from the observer. The
Doppler shift of the components’ spectral lines measured through a
spectrometer results in RV observations of the objects.

The RV $V$ at a certain time $t$ can be calculated as the sum of the RV of the
system´s center of mass $V_{CoM}$ (a constant since the system is assumed to
be free from external forces), and the radial part of the orbital velocity of
the observed component relative to the center of mass of the system,
$\dot{z}=dz/dt$:

$V(t)=V_{0}+\dot{z}(t).$ (A7)

Following the resolution of the two-body problem, the radial component of the
system can be expressed as:

$z(t)=r\sin(\nu(t)+\omega)\sin(i),$ (A8)

and therefore, by taking the first temporal derivative, we have that:

$\dot{z}(t)=\frac{2\pi
a\sin(i)}{P\sqrt{1-e^{2}}}[e\cos\omega+\cos(\nu(t)+\omega)].$ (A9)

Denoting $K=2\pi a\sin i/(P\sqrt{1-e^{2}})$ as the semi-amplitude of the RV
curve, the above expression becomes:

$V(t)=V_{0}+K[e\cos\omega+\cos(\nu(t)+\omega)].$ (A10)

Therefore, the RV $V(t)$ is characterized by six orbital parameters:

  * •

Period ($P$): The revolution period in days.

  * •

Time of periastron passage ($T$): One epoch of passage through periastron,
typically expressed in Julian Date (J.D.).

  * •

Eccentricity ($e$): The numerical eccentricity.

  * •

Argument of periapsis ($\omega$): The periastron longitude, counted from the
maximum of the RV curve.

  * •

Semi-amplitude ($K$): The semi-amplitude of the RV curve in km/s.

  * •

Velocity of the center of mass, ($V_{0}$): The RV of the center of mas of the
system in km/s (sometimes referred to as the systemic velocity as well).

Although the parameters $P$ and $T$ do not appear directly in the expression
of the RV Equation (A10), they are implicit in the determination of the true
anomaly $\nu$ (A1). According to the convention for the units of the orbital
parameters involved, the semi-amplitude $K$ is measured in km/s, whereby,
$a\sin(i)$ must be measured in km and $P$ must be converted to seconds through
$P[s]=86400\cdot P[days]$.

Expression (A10) is valid for the relative orbit
$\vec{r}=\vec{r}_{2}-\vec{r}_{1}$, however, the RV observations are relative
to the center of mass of the system. To correct this discrepancy, the semi-
major axis $a$ of the relative orbit must be replaced by their counterparts
$a_{1}$ and $a_{2}$ of the components relative to the center of mass of the
system.

By the definition of the center of mass of a system composed by two particles
of mass $m_{1}$ and $m_{2}$, the following relation is directly obtained:

$\begin{split}m_{1}a_{1}=m_{2}a_{2},\end{split}$ (A11)

and noticing that $a=a_{1}+a_{2}$, the following expressions for $a_{1}$ and
$a_{2}$ are obtained:

$\begin{split}a_{1}&=a\cdot\frac{m_{2}}{m_{1}+m_{2}}=a\cdot\frac{q}{1+q},\\\
a_{2}&=a\cdot\frac{m_{1}}{m_{1}+m_{2}}=a\cdot\frac{1}{1+q},\end{split}$ (A12)

where $q=m_{2}/m_{1}$ is defined as the mass ratio. Consequently, the semi-
amplitude of each component becomes:

$\begin{split}K_{1}&=\frac{1}{86400}\frac{2\pi
a_{1}\sin(i)}{P\sqrt{(1-e^{2})}},\\\ K_{2}&=\frac{1}{86400}\frac{2\pi
a_{2}\sin(i)}{P\sqrt{(1-e^{2})}}.\end{split}$ (A13)

Finally, noting that the argument of periapsis of the primary $\omega_{1}$ and
companion $\omega_{2}$ components of the system differs by 180°, we have that
$\omega=\omega_{1}=\omega_{2}+\pi$, the expression for the RV of each
component of the binary system becomes:

$\begin{split}V_{1}(t)&=V_{0}+K_{1}[e\cos(\omega)+\cos(\nu(t)+\omega)],\\\
V_{2}(t)&=V_{0}-K_{2}[e\cos(\omega)+\cos(\nu(t)+\omega)].\end{split}$ (A14)

It is important to note that the terms $a_{1}\sin(i)$ and $a_{2}\sin(i)$ in
the definition of the semi-amplitudes (A13) can not be separated through RV
observations.

When the spectra of both components are distinguishable, i.e., the RV of the
primary and the companion object are observable, the system is denoted as
$SB2$ and is characterized by the set of orbital parameters
$\vartheta_{SB2}=\\{P,T,e,K_{1},K_{2},V_{0}\\}$. However, this is an
infrequent case since the most of the spectroscopic binary systems
observations are only from the primary (brighter) object ($\sim 80\%$). When
the primary object spectra is the only observable, the system is denoted as
$SB1$ and is characterized by the set of orbital parameters
$\vartheta_{SB2}=\\{P,T,e,K_{1},V_{0}\\}$, since the parameter $K_{2}$ is
undetermined.

###  A.3 Visual-spectroscopic binary systems

Visual-spectroscopic binary systems corresponds to binary systems in which the
relative position and the RV of its components are observable. Since four
orbital parameters $P,T,e,\omega$ are common in the visual and spectroscopic
binary systems, the dynamical equations are coupled and a joint modeling that
describes both sources of information (positional and RV observations) allows
to lift the ambiguities and indeterminacy of each individual set of equations.
Furthermore, the joint modeling allows to determine the individual masses of
the system, and its parallax, called orbital parallax.

RV observations allows us to solve the indeterminacy of the longitude of the
ascending node $\Omega$ in the visual binary case, since the maximum/minimum
of the RV curve of each component is reached in the ascending/descending node.
Conversely, the positional observations allow us to decouple the term
$a\sin(i)$ in the spectroscopic binary case, since we can determine the
inclination $i$ of the orbit.

The individual masses of a binary system can be computed through the
determination of the total mass of the system $m_{1}+m_{2}$ and the mass ratio
$q=m_{2}/m_{1}$. According to the Third Law of Kepler, the total mass of the
system is obtained through:

$m_{1}+m_{2}=\frac{a[AU]^{3}}{P^{2}},$ (A15)

where $a[\text{AU}]$ is the relative semi-major axis of the system (in
astronomical units), $P$ the period of the system (in seconds) and
$m_{1},m_{2}$ the mass of the primary and the companion object (in solar
masses), respectively. Since positional observations in the plane of the sky
only allow us to determine the semi-major axis $a$ in angular units (seconds
of arc), the conversion to linear distance units (AU) is determined by the
following expression:

$a[AU]=\frac{a[^{\prime\prime}]}{\varpi},$ (A16)

where $\varpi$ is the system parallax (in seconds of arc), which becomes then
an additional orbital parameter required to determine the individual masses.

For the computation of the mass ratio $q$, the combination of positional and
RV observations are required. Considering that the RV observations can
determine the terms $a_{1}\sin(i)$ and $a_{2}\sin(i)$ and the positional
observations can determine the inclination $i$, both sources of observations
allow to determine the individual semi-major axis of each component ($a_{1}$
and $a_{2}$) and the mass ratio $q=a_{1}/a_{2}$ (A11).

Considering Equation (A16), the RV expression in (A14) becomes999The units
conversion between the common orbital parameters of Section A.1 and Section
A.2 is omitted for simplicity.:

$\begin{split}V_{1}(t)&=V_{0}+\dfrac{2\pi
a_{1}\sin{i}}{P\sqrt{1-e^{2}}}[\cos(\omega+\nu(t))+e\cos{\omega}],\\\
V_{2}(t)&=V_{0}-\dfrac{2\pi
a_{2}\sin{i}}{P\sqrt{1-e^{2}}}[\cos(\omega+\nu(t))+e\cos{\omega}],\end{split}$
(A17)

with $a_{1}=a^{\prime\prime}/\varpi\cdot q/(1+q)$,
$a_{2}=a^{\prime\prime}/\varpi\cdot 1/(1+q)$ and $a^{\prime\prime}$ the semi-
major axis in arcseconds.

If the RV observations of each component ($V_{1}(t)$ and $V_{2}(t)$) are
available ($SB2$ case), the combined model that describes the positional and
RV observations is characterized by the set of orbital parameters
$\vartheta_{VB+SB2}=\\{P,T,e,a,\omega,\Omega,i,V_{0},\varpi,q\\}$. However, if
the RV observations of only one component are available ($SB1$ case), the
parameters $q$ and $\varpi$ cannot be simultaneously determined.

##  Appendix B The Hamiltonian Monte Carlo and the No-U-Turn sampler

###  B.1 Hamiltonian Monte Carlo

The Hamiltonian Monte Carlo (Neal et al., 2011; Betancourt, 2017), also known
as Hybrid Monte Carlo, is an instance of the Metropolis-Hastings algorithm
that makes use of the geometry of the target probability distribution to guide
the transitions of the Markov chain. This allows to perform the sampling very
efficiently, avoiding a random walk behavior of the solution, and overly
sensitivity to correlated parameters. These features facilitate convergence on
high-dimensional target distributions much more quickly than with other
simpler methods, such as the random walk Metropolis-Hastings (Metropolis et
al., 1953; Hastings, 1970) or the Gibbs sampler (Geman & Geman, 1984).

The core of the Hamiltonian Monte Carlo method is to sample from the zone of
the parameter space that contributes highly to the computation of the
expectation of a target distribution $P(x)$ given a parametrization $f(x)$ of
the parameter space $\mathcal{X}$:

$\mathbb{E}_{P}(f)=\int_{\mathcal{X}}f(x)P(x)dx,$ (B1)

or in other words, to sample from the area of the parameter space of highest
mass $P(x)dx$ given a parametrization $f(x)$. This zone is denoted as the
typical set.

For that purpose, the transitions of the Markov chain must be guided by a
vector field in the direction of the typical set by exploiting the
differential structure of the target distribution. Hence, the vector field is
generated by using the gradient of the target distribution together with
auxiliary parameters of momentum that compensates the attractive force of the
gradient to the target distribution mode, preserving a dynamical equilibrium
that allows to align the generated vector field with the typical set.

A conservative dynamic in physical systems requires that any compression or
expansion in the position space must be compensated with a respective
expansion or compression in the momentum space, preserving the volume in the
joint space of position and momentum. To ensure this conservative dynamic
behavior, the transition probabilities of the chain follows the Hamiltonian
dynamics. The Hamiltonian dynamical system is described by a function over the
position $x$ and momentum $p$ variables, known as the Hamiltonian function
$H(x,p)$.

Let $x_{n}\in\mathbb{R}^{d}$ be the vector of parameters of the space state
and $P(x)$ the target distribution, each dimension of the space state is
complemented by a fictitious momentum variable $x$:

$x_{n}\rightarrow(x_{n},p_{n}),$ (B2)

where $p_{n}\in\mathbb{R}^{d}$. The combined space of the parameters
$(x_{n},p_{n})\in\mathbb{R}^{2d}$ is denoted as the phase space and the
respective induced distribution $P(x,p)$ is denoted as the canonical
distribution.

To mimic the conservative dynamic behaviour of the space variables and the
momentum variables, the canonical distribution is written in terms of the
Hamiltonian function101010It follows the Boltzmann canonical distribution
$P(x)=z^{-1}e^{-E(x)/t}$, with $z$ a normalization constant and $t$ the
temperature variable fixed to one.:

$P(x,p)=e^{-H(x,p)},$ (B3)

which implies that:

$H(x,p)=-\log{P(x,p)},$ (B4)

Hence, the Hamiltonian function captures the probabilistic structure of the
phase space, and consequently, the geometry of its typical set.

A marginalization of the canonical distribution $P(x,p)$ in terms of the state
variable $p$ induces the following decomposition of the Hamiltonian function
$H(x,p)$:

$\begin{split}H(x,p)&=-\log P(p|x)-\log P(x)\\\ &\equiv
K(p,x)+V(x).\end{split}$ (B5)

The decomposition can be interpreted in terms of a kinetic energy $K(p,x)$
function, dependant upon both the spatial and momentum variables, and a
potential energy function $V(x)$, dependant upon the momentum variables only.
The potential function is simply the negative logarithm of the target
distribution, while the kinetic energy is usually expressed as a quadratic
term on $p$:

$K(p)=\frac{1}{2}p^{T}\cdot M^{-1}\cdot p,$ (B6)

where $M$ is a symmetric, positive-definite matrix denoted as mass-matrix. The
mass matrix is typically a scalar multiple of the identity matrix, but can
explicitly depends on $x$ as in (B5).

With these elements, the vector field oriented in the direction of the typical
set can be defined through the Hamiltonian equations:

$\begin{split}\frac{dx}{dt}&=+\frac{\partial H}{\partial p}=\frac{\partial
K}{\partial p}\\\ \frac{dp}{dt}&=-\frac{\partial H}{\partial
x}=-\frac{\partial K}{\partial x}-\frac{\partial V}{\partial x}.\end{split}$
(B7)

Following the vector field (determined by the Hamiltonian equations for a time
$t$) we can generate trajectories $\phi_{t}(x,p)$ that move along the typical
set. To compute these trajectories, the solution of (B7) is obtained by
numerical methods. In particular, the trajectory $\phi_{T}(x,p)$ of a variable
$i$ can be approximated by the leap-frog integration method iterating the
following expressions:

$\begin{split}p_{i}(t+\epsilon/2)&=p_{i}(t)-\frac{\epsilon}{2}\frac{\partial
V}{\partial x_{i}}(x_{i}(t))\\\
x_{i}(t+\epsilon)&=x_{i}(t)+\frac{\epsilon}{m_{i}}p_{i}(t+\epsilon/2)\\\
p_{i}(t+\epsilon)&=p_{i}(t+\epsilon/2)-\frac{\epsilon}{2}\frac{\partial
V}{\partial x_{i}}(x_{i}(t+\epsilon)),\,\,\,\,\forall
i\in{1,...,L}.\end{split}$ (B8)

$L$ is the number of steps, $\epsilon\in\mathbb{R}$ is the step size and
$T=\lfloor L/\epsilon\rfloor$ the integration time. The adequate selection of
the algorithm hyper-parameter $\epsilon$ and $L$ is crucial for a good
sampling performance.

In summary, the Markov chain that samples from the target distribution $P(x)$
will follow the Metropolis-Hastings algorithm defined on the phase space
$(x,p)$ with transition probabilities $T_{(x_{0},p_{0}),(x_{L},p_{L})}$
determined by the solution of the Hamiltonian equations following the leap-
frog integration method for a fixed number of steps $L$ and step size
$\epsilon$. The momentum variables are sampled from a proposal marginal
distribution $P(p|x)$ and the final samples of $P(x)$ are obtained by
projecting the samples of the phase space on the state space $(x,p)\rightarrow
x$, i.e., ignoring the momentum variables.

The transition probabilities of the current Markov Chain must be modified
since the transition ratio in the Metropolis-Hastings acceptance probability
$T_{(x_{L},p_{L}),(x_{0},p_{0})}/T_{(x_{0},p_{0}),(x_{L},p_{L})}=0/1=0$,
because the leap-frog integration does not allow reverse trajectories. Thus,
the transition probabilities are modified to be reversible by augmenting the
numerical integration with a negation step that flips the sign of momentum
$(x,p)\rightarrow(x,-p)$. Thereby, the Metropolis-Hastings acceptance rate
becomes:

$\begin{split}\mathcal{A}((x_{L},-p_{L}),(x_{0},p_{p}))&=\min\left\\{1,\frac{T_{(x_{L},-p_{L}),(x_{0},p_{0})}P(x_{L},-p_{L})}{T_{(x_{0},p_{0}),(x_{L},-p_{L})}P(x_{0},p_{0})}\right\\}\\\
&=\min\left\\{1,\frac{\delta(x_{L}-x_{L})\delta(-p_{l}+p_{L})P(x_{L},-p_{L})}{\delta(x_{0}-x_{0})\delta(p_{0}-p_{0})P(x_{0},p_{0})}\right\\}\\\
&=\min\left\\{1,\frac{P(x_{L},-p_{L})}{P(x_{0},p_{0})}\right\\}\\\
&=\min\left\\{1,e^{-H(x_{L},-p_{L})+H(x_{0},p_{0})}\right\\}.\end{split}$ (B9)

The complete Hamiltonian Monte Carlo sampling procedure is described in
Algorithm LABEL:alg:hmc.

###  B.2 No-U-Turn sampler

The No-U-Turn Sampler algorithm (Hoffman et al., 2014) is an extension of the
Hamiltonian Monte Carlo algorithm that adaptively sets the number of steps $L$
of the trajectories, facilitating the use of the sampling tool by avoiding a
low performance selection of the Hamiltonian Monte Carlo user-defined hyper-
parameters. The selection of the number of steps $L$ follows the computation
of forward and backward exploration trajectories of the Hamiltonian Monte
Carlo algorithm until an end condition is met, where the new sample is
obtained by a random selection of the generated trajectories.

To generate the exploration trajectories, a binary tree is constructed
iteratively. Let $(x_{n}(0),p_{n}(0))$ be an initial particle composed by a
position and momentum of the $n$-th iteration of the Markov chain,
$(x_{n}^{+},p_{n}^{+})$ be a forward in time particle and
$(x_{n}^{-},p_{n}^{-})$ be a backward in time particle. In each iteration $j$,
the binary tree selects at random uniformly to move the $(j-1)$-particle
forwards or the backward particle backwards in time, with $2^{j}$ number of
leap-frog integration steps.

The iterative procedure continues until the following condition (namely the
U-Turn condition) is met:

$(x_{n}^{+}-x_{n}^{-})\cdot p_{n}^{-}<0\,\,\vee\,\,(x_{n}^{+}-x_{n}^{-})\cdot
p_{n}^{+}<0,$ (B10)

or when the Hamiltonian trajectory generated by the leap-frog integration
becomes imprecise in the sense that:

$\begin{gathered}e^{-H(x_{n}^{+},p_{n}^{+})+\Delta_{max}}<U_{n}\,\,\vee\,\,e^{-H(x_{n}^{-},p_{n}^{-})+\Delta_{max}}<U_{n},\end{gathered}$
(B11)

where $U_{n}\sim\mathcal{U}\left(0,e^{-H(q_{n}(0),p_{n}(0))}\right)$ is a
slice random variable sample and $\Delta_{max}$ is a maximum energy hyper-
parameter. The idea of the No-U-Turn condition is avoiding the generation of
redundant trajectories by stopping the exploration when the trajectory begins
to turn back to previous explored zones.

Finally, the new sample $(x_{n+1},p_{n+1})$ is selected by a uniform sampling
of the generated trajectory that satisfies the precision condition
$U_{n}<e^{-H(x_{n+1},p_{n+1})}$. The complete No-U-Turn sampler sampling
procedure is described in Algorithm LABEL:alg:nuts.

###  B.3 Implementation considerations in the binary stars context

The implementation of the proposed Bayesian inference methodology based on the
No-U-Turn sampler algorithm requires to compute the gradient of the posterior
function, which can be analytically derived from the Keplerian model formulae
presented in Section A by taking the partial derivatives with respect to each
orbital parameter that characterizes the binary stellar system. However, a
special consideration must be taken with the partial derivatives of the
eccentric anomaly $E$, since it is not analytically calculated, but rather
numerically approximated. The non-zero partial derivatives of the eccentric
anomaly can be expressed as a function of the variable itself as follows:

$\begin{split}\frac{\partial E}{\partial e}&=\frac{\sin{E}}{1-e\cos{E}},\\\
\frac{\partial E}{\partial T}&=-\frac{2\pi}{P}\cdot\frac{1}{1-e\cos{E}},\\\
\frac{\partial E}{\partial
P}&=-\frac{2\pi(t-T)}{P^{2}}\cdot\frac{1}{1-e\cos{E}},\end{split}$ (B12)

with the value of $E$ previously approximated by any numerical method (e.g.,
the Newton-Raphson method (Ypma, 1995)).

## References

  * Abushattal et al. (2020) Abushattal, A. A., Docobo, J. A., & Campo, P. P. 2020, AJ, 159, 28, doi: 10.3847/1538-3881/ab580a
  * Agati et al. (2015) Agati, J. L., Bonneau, D., Jorissen, A., et al. 2015, A&A, 574, A6, doi: 10.1051/0004-6361/201323056
  * Aller et al. (1996) Aller, L., Appenzeller, I., Baschek, B., et al. 1996, lbor, 3
  * Betancourt (2017) Betancourt, M. 2017, arXiv preprint arXiv:1701.02434
  * Bouffanais & Porter (2019) Bouffanais, Y., & Porter, E. K. 2019, Physical Review D, 100, 104023
  * Brown et al. (2018) Brown, A., Vallenari, A., Prusti, T., et al. 2018, Astronomy & astrophysics, 616, A1
  * Burgasser et al. (2012) Burgasser, A. J., Luk, C., Dhital, S., et al. 2012, The Astrophysical Journal, 757, 110
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., et al. 2017, Journal of statistical software, 76, 1
  * Carroll & Ostlie (2006) Carroll, B. W., & Ostlie, D. A. 2006, An introduction to modern astrophysics and cosmology
  * Cid Palacios (1958) Cid Palacios, R. 1958, AJ, 63, 395
  * Claveria et al. (2019) Claveria, R. M., Mendez, R. A., Silva, J. F., & Orchard, M. E. 2019, PASP, 131, 084502, doi: 10.1088/1538-3873/ab22e2
  * Docobo (1985) Docobo, J. 1985, Celestial Mechanics, 36, 143
  * Docobo et al. (1992) Docobo, J., Ling, J., & Prieto, C. 1992, in International Astronomical Union Colloquium, Vol. 135, Cambridge University Press, 220–222
  * Docobo et al. (2018) Docobo, J. A., Tamazian, V. S., Campo, P. P., & Piccotti, L. 2018, The Astronomical Journal, 156, 85
  * Ford (2005) Ford, E. B. 2005, The Astronomical Journal, 129, 1706
  * Geman & Geman (1984) Geman, S., & Geman, D. 1984, IEEE Transactions on pattern analysis and machine intelligence, 721
  * Gray (2008) Gray, D. F. 2008, The Observation and Analysis of Stellar Photospheres
  * Gregory (2005) Gregory, P. 2005, The Astrophysical Journal, 631, 1198
  * Gregory (2011) Gregory, P. C. 2011, Monthly Notices of the Royal Astronomical Society, 410, 94
  * Habets & Heintze (1981) Habets, G., & Heintze, J. 1981, Astronomy and Astrophysics Supplement Series, 46, 193
  * Hajian (2007) Hajian, A. 2007, Physical Review D, 75, 083525
  * Halbwachs et al. (2016) Halbwachs, J.-L., Boffin, H., Le Bouquin, J.-B., et al. 2016, Monthly Notices of the Royal Astronomical Society, 455, 3303
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Hastings (1970) Hastings, W. K. 1970
  * Hoffman et al. (2014) Hoffman, M. D., Gelman, A., et al. 2014, J. Mach. Learn. Res., 15, 1593
  * Hou et al. (2012) Hou, F., Goodman, J., Hogg, D. W., Weare, J., & Schwab, C. 2012, The Astrophysical Journal, 745, 198
  * Hummel et al. (1994) Hummel, C., Armstrong, J., Quirrenbach, A., et al. 1994, The Astronomical Journal, 107, 1859
  * Hunter (2007) Hunter, J. D. 2007, Computing in science & engineering, 9, 90
  * Ji et al. (2017) Ji, Y., Banks, T., Budding, E., & Rhodes, M. 2017, Astrophysics and Space Science, 362, 1
  * Kumar et al. (2019) Kumar, R., Carroll, C., Hartikainen, A., & Martin, O. 2019, Journal of Open Source Software, 4, 1143, doi: 10.21105/joss.01143
  * Liu & Nocedal (1989) Liu, D. C., & Nocedal, J. 1989, Mathematical programming, 45, 503
  * Lucy (2014) Lucy, L. 2014, Astronomy & Astrophysics, 563, A126
  * Lucy (2018) —. 2018, Astronomy & Astrophysics, 618, A100
  * MacKnight & Horch (2004) MacKnight, M., & Horch, E. 2004, AAS, 204, 07
  * Mason (2015) Mason, B. D. 2015, in IAU General Assembly, Vol. 29, 2300709
  * Mason et al. (2001) Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466, doi: 10.1086/323920
  * Mendez et al. (2021) Mendez, R. A., Clavería, R. M., & Costa, E. 2021, AJ, 161, 155, doi: 10.3847/1538-3881/abdb28
  * Mendez et al. (2017) Mendez, R. A., Claveria, R. M., Orchard, M. E., & Silva, J. F. 2017, The Astronomical Journal, 154, 187
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, The journal of chemical physics, 21, 1087
  * Morbey (1975) Morbey, C. 1975, Publications of the Astronomical Society of the Pacific, 87, 689
  * Moré (1978) Moré, J. J. 1978, in Numerical analysis (Springer), 105–116
  * Muterspaugh et al. (2010) Muterspaugh, M. W., Hartkopf, W. I., Lane, B. F., et al. 2010, The Astronomical Journal, 140, 1623
  * Neal et al. (2011) Neal, R. M., et al. 2011, Handbook of markov chain monte carlo, 2, 2
  * Nelson et al. (2013) Nelson, B., Ford, E. B., & Payne, M. J. 2013, The Astrophysical Journal Supplement Series, 210, 11
  * Pourbaix (1994) Pourbaix, D. 1994, Astronomy and Astrophysics, 290, 682
  * Pourbaix (1998) —. 1998, Astronomy and Astrophysics Supplement Series, 131, 377
  * Pourbaix (2000) Pourbaix, D. 2000, A&AS, 145, 215, doi: 10.1051/aas:2000237
  * Pourbaix et al. (2004) Pourbaix, D., Tokovinin, A. A., Batten, A. H., et al. 2004, A&A, 424, 727, doi: 10.1051/0004-6361:20041213
  * Prusti et al. (2016) Prusti, T., De Bruijne, J., Brown, A. G., et al. 2016, Astronomy & Astrophysics, 595, A1
  * Sahlmann et al. (2013) Sahlmann, J., Lazorenko, P., Ségransan, D., et al. 2013, Astronomy & Astrophysics, 556, A133
  * Schmidt-Kaler et al. (1982) Schmidt-Kaler, T., Schaifers, K., & Voigt, H. 1982, Stars and Star Clusters, vol
  * Shabram et al. (2020) Shabram, M. I., Batalha, N., Thompson, S. E., et al. 2020, The Astronomical Journal, 160, 16
  * Skiff (2014) Skiff, B. A. 2014, VizieR Online Data Catalog, B/mk
  * Söderhjelm (1999) Söderhjelm, S. 1999, A&A, 341, 121
  * Straizys & Kuriliene (1981) Straizys, V., & Kuriliene, G. 1981, Ap&SS, 80, 353, doi: 10.1007/BF00652936
  * Thiele (1883) Thiele, T. N. 1883, Astronomische Nachrichten, 104, 245
  * Tokovinin (1992) Tokovinin, A. 1992, Astronomy and Astrophysics, 256, 121
  * Wang et al. (2009) Wang, Q., Kulkarni, S. R., & Verdú, S. 2009, IEEE Transactions on Information Theory, 55, 2392
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Wenger et al. (2000a) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000a, A&AS, 143, 9, doi: 10.1051/aas:2000332
  * Wenger et al. (2000b) —. 2000b, A&AS, 143, 9, doi: 10.1051/aas:2000332
  * Ypma (1995) Ypma, T. J. 1995, SIAM review, 37, 531
  * Ziegler et al. (2018) Ziegler, C., Law, N. M., Baranec, C., et al. 2018, AJ, 156, 259, doi: 10.3847/1538-3881/aad80a