Estimating Microlensing Parameters from Observables and Stellar Isochrones with pyLIMASS

We present pyLIMASS, a novel algorithm for estimating the physical properties of the lensing system in microlensing events. The main idea of pyLIMASS is to combine all available information regarding the microlensing event, defined as observables, and to estimate the parameter distributions of the system, such as the lens mass and distance. The algorithm is based on isochrones for the stars model and combines the observables using a Gaussian mixture approach. After describing the mathematical formalism and its implementation, we discuss the algorithm’s performance on simulated and published events. Generally, the pyLIMASS estimations are in good agreement (i.e., within 1σ) with the results of the selected published events, making it an effective tool to estimate the lens properties and their distribution. The applicability of the method was tested by using a catalog of realistically simulated events that could be observed by the future Galactic Bulge Time Domain Survey of the Nancy Grace Roman Space Telescope. By solely using constraints from the Roman lightcurves and images, pyLIMASS estimates the masses of the lens of the Roman catalog with a median precision of 20% with almost no bias.


INTRODUCTION
Measuring the mass of astrophysical objects is crucial to understanding their distribution and abundance in the Milky Way.However, this is often challenging due to their large distances and the technical limitations.Gravitational microlensing (Einstein 1936;Paczynski 1986) is ideally suited to probe the mass function of objects in the entire Milky Way, because it observes the change in brightness of a background star, and the relatively faint lens object can be characterized solely by studying the light curve of the background star.Originally, the technique has been used to study lens populations of the the Milky Way halo by surveying the source located in the Magellanic Clouds (Paczynski 1986;Alcock et al. 2000;Tisserand et al. 2007;Wyrzykowski et al. 2011).More recently, microlensing surveys have been focused towards the Galactic Bulge (Bond et al. 2001;Udalski 2003;Kim et al. 2016), where the event rate is highest (Sumi et al. 2013;Mróz et al. 2019), but microlensing events have also been detected in the entire Galactic Disk (Fukui et al. 2019;Wyrzykowski et al. 2023).Based on these observations, more than two hundred planets have been detected 1 , as well as brown dwarfs (Zhu et al. 2016;Chung et al. 2017;Shvartzvald et al. 2019;Bachelet et al. 2019), free-floating planets (Mróz et al. 2017(Mróz et al. , 2018;;Sumi et al. 2023), stellar binaries (see for example Street et al. (2019); Tsapras et al. (2019)) and stellar remnants (Blackman et al. 2021;Sahu et al. 2022).Moreover, the galactic structure can be inferred from the distribution of microlensing event parameters (Calchi Novati et al. 2015;Awiphan et al. 2016;Mróz et al. 2017Mróz et al. , 2019)).Given its potential, a microlensing survey is one of the core missions of the upcoming Nancy Grace Roman Space Telescope (Spergel et al. 2015;Penny et al. 2019;Johnson et al. 2020).
One of the difficulties of the method is to reconstruct the physical properties of the microlens, especially its mass and distance.Generally, the lightcurve of the microlensing event is the main phenomenon observed.From the modeling of the lightcurve, the Einstein ring crossing time t E is usually measured with great precision (Gould 2000): where µ rel = µ L − µ S is the relative proper motion vector between the lens and the source.It is measured in the observer frame (geocentric frame for example), and can be transformed to the heliocentric frame by (Skowron et al. 2011): Here V ⊕,⊥ is the projected observer velocity in the plane of sky at the time of the event maximum and π rel = π L − π S is the relative parallax (Gould 2000).θ E is the angular Einstein ring radius that is a fundamental characteristic of the lensing phenomenon and defined as (Gould 2000): with κ a constant and M L the lens mass in solar units.θ E can be measured if finite-source effects are seen in the lightcurve, commonly parameterized with ρ * = θ S θ E , where θ S is the angular source radius (Witt & Mao 1994;Yoo et al. 2004).We note that this parameter is also sometimes parameterized by t S = ρ * t E , which is a measure of the duration of finite-source effects.The observations of the microlensing astrometric shifts (Dominik & Sahu 2000;Sahu et al. 2022) or interferometric measurements during the peak of the event (Dong et al. 2019;Cassan et al. 2022) also lead to a direct measurement of θ E , but these observations are generally challenging and expensive in telescope time.The measurement of θ E can also be made by observing the lens and source several years after the event peak with high-resolution imagers (Batista et al. 2015;Bennett et al. 2015;Bhattacharya et al. 2018;Vandorou et al. 2020).Indeed, the observed lens and source angular separation δ at a time offset ∆t is directly related to θ E via: where π E is the norm of the microlensing parallax vector (Gould 2004).The microlensing parallax vector π E is crucial to measure because it contains information about the lens mass and distance as well as the direction of the relative proper motion, for example in the equatorial system (Gould 2004): (µ rel,geo,N , µ rel,geo,E ) The microlensing parallax can be measured from a single location if the event is long enough, e.g. the annual parallax effect (Gould 2004), or if the event is observed by several observatories with significant separations, e.g. the space based parallax (Refsdal 1966).This has be done many times, especially with the Spitzer telescope, see for example Calchi Novati et al. (2015) and Zang et al. (2020), and Gaia (Wyrzykowski et al. 2020).
The measurements of the lens brightness m L,λ , preferably in several bands or with spectroscopy, also provide strong constraints on the lens system (Beaulieu 2018): where A L,λ is the absorption towards the lens and M 0 is the absolute magnitude of the lens in a given band λ.The latter can be estimated via spectral templates (Kurucz 1993;Allard & Hauschildt 1995), stellar isochrones libraries (Bressan et al. 2012;Choi et al. 2016) or empirical relation (Pecaut & Mamajek 2013;Benedict et al. 2016;Mann et al. 2019;Rabus et al. 2019).
Another source of constraints can be obtained by using models of the Milky Way.This has been done on many occasions and for different models (Han & Gould 1995, 2003;Dominik 2006;Bennett et al. 2014;Mróz et al. 2020;Specht et al. 2020;Koshimoto et al. 2021).In particular, by comparing the timescale, optical depth and event rate distributions of microlensing events with Galactic Models predictions, constraints can be derived on the various population of objects as well as the kinematics and architecture of the Galaxy (Mróz et al. 2020;Specht et al. 2020;Koshimoto et al. 2021).These models can also be used to study individual event, by generating priors for the modeling of the lightcurves or to estimate the event physical parameters after the modeling.The downside is that the results are inevitably biased because of the different models assumptions.For example, Bachelet et al. (2022a) found a good agreement with the prediction from the Besançon (Robin et al. 2003), the GalMod (Pasetto et al. 2018), Dominik (2006) and Koshimoto et al. (2021) models, but small discrepancies exist.
Another difficulty is to combine the different measurements to reconstruct the properties of the lensing system.A natural way to do is to include any available exterior information with the form of prior during the modeling process.This is routinely done during the reanalysis of event with high-resolution imaging constraints (Bennett et al. 2015;Bhattacharya et al. 2018;Vandorou et al. 2020).However, the modeling of lightcurve can be non-trivial and time consuming, while the inclusion of all available information could also be difficult to add into modeling codes in practice.
In recent years, the microlensing community has made great efforts to publish public codes for the modeling of microlensing events.The first package to be released was pyLIMA (Bachelet et al. 2017), a python-based software that includes multiple tools for analyzing microlensing events and is now widely used by the community.MulensModel by Poleski & Yee (2019) is another public package for modeling microlensing events.Both packages make use of the VBBinaryLensing C++ library (Bozza 2010;Bozza et al. 2018), a contour-integration code for the estimation of the microlenisng magnifcation.This library is also in use in the RTModel infrastructure, which provides real-time modeling for microlensing events since 20132 and has been recently released to the public3 .Similarly, Cheongho Han provides access to the community to real-time models made by his team4 .More recently, the software eesunhong have been made public5 .This code has been used for many years to analyze plenty of microlensing events (Bennett & Rhie 1996;Bennett 2010;Bennett et al. 2023).
As far as we know, none of these software provides tools to estimate the physical properties of the lensing system.While this is sometimes (relatively) straightforward to do so, when the microlensing parallax and finite-source effects are measured in the lightcurve for example, the combination of all the available information to accurately reconstruct the physical properties of the lenses is non trivial.Therefore, we propose in this work a new tool to achieve this goal.The algorithm is solely based on the information in observables, as defined in the next section and provided by the user with, by default, the additional prior of the stellar isochrones.The purpose of this tool is to combine the provided measurements (and stellar isochrones) to estimate the lensing system physical parameters in a timely manner.We note that this estimate is done by default without using any prior from Galactic Models.The algorithm is described in Section 2 as well as its python implementation in pyLIMASS.Section 3 and Section 4 present the algorithm performance on simulated data and published analysis and we present our conclusions and potential upgrades in Section 5.

Definition of the observables
In the following, we define an observable x as a measurable physical quantity, and its associated uncertainty, related to a microlensing event.The list of default observables defined in pyLIMASS can be found in the Appendix, but as an example the source magnitudes in different bands or the relative proper motion vector can be considered as observables.In general, an observable can be any kind of measurable quantity related to the microlensing event and does not have to be directly related to the modeling of the microlensing event data (photometric lightcurve, astrometric time series and/or high resolution imaging).For example, if the parallax of the lens is known directly from the Gaia satellite measurements (Gaia Collaboration et al. 2016, 2023), then this is valuable information that is relevant to include as an observable.

Inverse problem
As discussed in the introduction, the parameters of most interest P , such as the lens mass and distance, are generally not measured directly and need to be reconstructed from a set of observables X = (x 1 , x 2 , ....x i ).This is a classic inverse problem and commonly expressed through Bayes' theorem in terms of probabilities: p(X|P ) is the likelihood, p(P ) is the prior probability and p(X) is a constant often referred as model evidence.The goal is to sample the posterior probabilities p(P |X) ∝ p(X|P )p(P ), i.e. estimate the more plausible parameters P that replicates the set of observables X.The choice of priors is also important to sample the posterior probabilities.This is further discussed in the Appendix A.

Gaussian Mixture
Equation 1 shows that measuring only t E in a microlensing source event, yields weak, but not null, constraints on the lens mass.In fact, studying the parameters distributions of multiple events provides valuable information about the architecture of the Milky Way (Sumi et al. 2011;Calchi Novati et al. 2015;Mróz et al. 2019;Sumi et al. 2023).But the primary goal of our approach is to gather as much information as possible, expressed in terms of the observables as defined previously, about one specific microlensing event to reconstruct the physical parameters of the source and lens system.Ideally, it is desirable that the algorithm could handle constraints from multiple sources and of various forms.For example, t E is generally well measured from the lightcurve and leads to a constraint of the form N (t E , σ t E ).Similarly, the blend flux is also generally well constrained from the lightcurve only.However, the hypothesis that the additional flux comes solely from the lens often turns to be wrong, especially towards the Galactic Bulge where the crowding is extreme.In fact, many high-resolution follow-up studies revealed unrelated companions at ≤ 1" from the target, see for example Blackman et al. (2021); Bhattacharya et al. (2021); Terry et al. (2021).In that case, a safer assumption is to assume a constraint (for example in the V band) on the form V L ≥ V blend .
Therefore, we seek for a mechanism that supports such various distribution of constraints and also provide the likelihood P (X|P ) (or at least an approximation) of a set of parameters P given the set of observables X.To combine these observables, we model them using a multidimensional Gaussian Mixture (GM).GMs are frequently used for this purpose (Kelly 2007;Hao et al. 2010;Alsing et al. 2018;Frühwirth-Schnatter et al. 2019;Ansari et al. 2021), and the empirical distribution of the observables X is then represented by G Gaussian components.The joint probability distribution of a single observation y is then given by (Frühwirth-Schnatter et al. 2019): where θ = (ϕ k , µ k , Σ k ) is a set of G hyperparameters.µ k and Σ k are the mean vector and covariance matrix associated with the multidimensional Gaussian component k.The parameter ϕ k is the weight of the Gaussian components and by definition The goal of the GM model is to approximate the underlying probability distribution of a sample of N observables X = (x 1 , x 2 , ...., x N ) or, in other words, find the parameters θ f that maximize the likelihood: This is a non-trivial problem that is an area of active research (Jin et al. 2016;Lucic et al. 2018), but it can be solved by using the Expectation-Maximization (EM) algorithm (Dempster et al. 1977) for example.Once the best hyperparameters θ f are found, the likelihood of the observables serie X for the parameters P can be approximated via the GM model: It is worth noting that assuming the simplest case with no correlation between the parameters, Equation 10 simplifies to: In this equation µ j,k and σ j,k represent the mean and the standard deviation of the Gaussian component k and the observable j, while d represents the total number of observables (i.e. the number of dimension).If the generated observable is close to the most probable region of a specific Gaussian k (i.e.y ∼ µ k ), the log-likelihood becomes: Here, it is assumed for simplicity that all weights ϕ are zero except for the value of the k component.Although this is an oversimplified model, it effectively demonstrates the general behaviour of the log-likelihood implemented in pyLIMASS.
The great advantage of GM is its simplicity as well as its potential to approximate multi-modal distributions with complex shapes (Frühwirth-Schnatter et al. 2019).Such distributions are quite common in microlensing, where parameters degeneracies or competitive models can be found.To illustrate this, we realized an ad hoc simulation, where two random variables U = 0.1±0.05 and V = 300U 2 are somehow measured.The resulting distribution is visible in Figure 1.We then modeled different GM models using G components.It is clear that the original distribution is better replicate as the number of Gaussian is increased.However, there is a risk of overfitting by increasing arbitrarily the number of components.While the optimal number of components is problem dependent, statistical tools like the Bayesian information criterion (BIC, Schwarz (1978)) or the Akaike Information Criterion (AIC, Akaike (1974)) can be used to assess the significance of additional components.For the example presented, ∼ 20 components seems sufficient according to both metrics.

Implementation in pyLIMASS
We have implemented this approach in pyLIMASS, a new module of the pyLIMA package (Bachelet et al. 2017).First, it is worth noting that most of the relations presented earlier contain products, ratios, and power-laws that are much better represented in logarithmic space.Moreover some of the observables, such as ρ * or t E , are strictly positive.The GM algorithm also performs better in logarithmic space because it represents the actual distribution more accurately and is not affected by floor and ceiling effects.Furthermore, the previously described distributions of parameters, such as ρ * , are transformed into a sum of distributions, and consequently behave much better in the logarithmic space.As a result, the algorithm utilizes logarithmic quantities whenever possible.In pyLIMASS, GM models are generated with the scikit-learn implementation (Pedregosa et al. 2011), and the number of Gaussian functions used can be adjusted by the user, although it is set by default to the number of observable distributions ingested.Eight parameters P are defined for the source and lens, namely the mass M in solar mass units, the distance D in kpc, the effective temperature T ef f in Kelvin, the metallicity F e in dex, the logarithmic surface gravity logg, the geocentric North and East proper motions µ N and µ E in mas/yr and the V-band absorption A V .For the lens body, it is more efficient to sample ϵ D = D L /D S ≤ 1 and ϵ A V = A V,L /A V,S ≤ 1, these quantities are therefore used in pyLIMASS instead of D L and A V,L .Most of the observables can be reconstructed from these parameters and compared with the GM distribution model.
The reconstruction of the observed magnitudes of the source and, possibly, the lens is the main notable exception.To estimate these observables, pyLIMASS uses an interpolation of stellar isochrones.In particular, the models described in Bressan et al. (2012), available online6 , that includes a broad range of age and metallicity values as well as a large catalog of passbands, have been implemented.By default (and for the rest of this study), pyLIMASS considers only stars with an age ≥ 1 Gyr, which can be adjusted to meet the user's needs.We would like to point out that pyLIMASS also supports the option to compute the magnitudes directly from spectral templates using the Spyctres library7 .In short, Spyctres relies on stsynphot and synphot to simulate spectra of stars with the Kurucz (Kurucz 1993) or Phoenix (Allard & Hauschildt 1995) libraries (depending on the user's choice).However, this option is significantly slower and is not used in the following.The blackbody approximation option is also available, as it is much faster and still provides coherent results, except for the cooler stars.To estimate the absorption in the needed bands, pyLIMASS uses the absorption law of Wang & Chen (2019).
The observables x are ingested in pyLIMASS as samples.If one has only knowledge of the mean µ and standard deviation σ of given observables, then a sample of N (µ, σ) is generated independently for each observable.However, if one has some prior knowledge about the correlation between observables, joint distributions can be used instead.For example, if the correlation between the parameters π EN and π EE is known from MCMC modeling of a specific event, then the corresponding samples should be included.
To explore the parameter space, two algorithms are available by default in pyLIMASS.The scipy implementation (Virtanen et al. 2020) of the differential evolution algorithm (Storn & Price 1997) is implemented to locate the region of maximum likelihood and the emcee package (Foreman-Mackey et al. 2013) is available to run Monte-Carlo Marko Chain explorations.The description of the default parameters priors can be found in the Appendix.It is however straightforward for users to use their preferred sampler.
It is useful to know if a given model is plausible.One can compare the model likelihood L to the GM model probability when y = µ k .Indeed, models centered close to the GM means will have a high likelihood (but potentially not the maximum).Noting L the likelihood at y = µ k , the quantity −2(L − L) is chisquare distributed with d degrees of freedom, where d is the number of observables used.Therefore, every model likelihood can be compared to L using a significance level α.We have implemented such afunctionality in pyLIMASS and set a default of α = 0.05.
To conclude, the design of the algorithm can be found in the Figure 2 and can summarized in three steps as follow: • Collect the observables x and pass them to the algorithm in the form of distributions.It can be on the form of Normal distribution or MCMC chains for example • Fit the observed distribution with a GM model to estimate θ f , via the EM algoritm for example.
• Use the trained GM model to estimate the posterior probabilities of parameters P given the observables and priors via Equation 7 and Equation 103.TESTS ON SIMULATED EXAMPLES

A typical single lens event
In this section, we have simulated a "standard" microlensing event as it could be seen by the Roman microlensing survey.It is anticipated that Roman would observe the microlensing fields with at least two filters (Penny et al. 2019).Most of the images would be taken with the F146 wide filter but complementary observations with a ∼ daily cadence will be obtained with the narrow band F087 filter8 .However, we note that the exact strategy is still not decided at the time of writing9 and multiple scenarii exist, that potentially include extra narrow band filters.We set the source to be a Sun like star with M S = 1M ⊙ , R S = 1R ⊙ , T ef f = 5778 K, F e S = 0 dex, logg S = 4.437 cgs at D S = 8 kpc and with a visual absorption A V = 4 mag, with a (geocentric) proper motion µ S,N = µ S,E = −3 mas/yr.The lens is assumed to be an M-dwarf of cgs at D L = 4 kpc and a visual absorption A V = 2 mag, with a (geocentric) proper motion µ L,N = 2 mas/yr and µ L,E = −2 mas/yr.Using these information, it is straightforward to generate any observables such as t E , π E , the lens and source magnitudes in the F087 and F146 bands and so on.Based on this simulation, several reconstruction fits have been performed using different set of observables.This intends to reproduce real case events where some quantities or effects are measured or not.The uncertainties of all observables have been arbitrary set to 10% and without correlation.Results are summarized in Table 1.The first experiment (Fit 1) used constraints from θ E , D S , π EN and π EE .This scenario could occur if Roman data measures the source parallax π S as well as the astrometric shift induced by the microlensing event (Fardeen et al. 2023), and therefore measuring θ E directly, in addition of the measurement of the microlensing parallax vector π E from the lightcurve modeling.In this case, the lens mass and distance are reconstructed with high fidelity, as one should expect according to Equation 4.However, the mass of the lens cannot be reconstructed to better than 14.6% (i.e.±0.073M ⊙ , via error propagation) using only constraints from θ E , π EN and π EE , and this is exactly what is recovered in Table 1.Similarly, the lens distance is expected to be measure with a precision of 9 % (i.e.±0.36 kpc) because it was assumed that the source distance is (somehow) known with a 10 % precision.Fit 2 used θ S , ρ * , π EN and π EE constraints and the reported precision on the lens mass measurement decreased as expected (the anticipated precision is 17% in this case).But in this case the lens distance is not known because the source distance is only poorly constrained by θ S and ρ * .This scenario is realistic, as it is expected that ρ * should be measured on a significant number of events with lightcurves observed by the Roman space telescope, while θ S should be routinely derived from the study of the color-magnitudes diagrams of Roman fields (Yoo et al. 2004;Nataf et al. 2013).The last scenario (Fit 3) corresponds to a case where few measurements can be made from the lightcurve itself (a single lens single source event for example), but the relative proper motion vector µ rel as well as the magnitudes of the source and lens in two bands have been measured via the modeling of the Roman images at the location of the microlensing event.This technique, often referred as "lens-flux analysis", is expected to provide constraints for a large fraction of events.In this case, constraints from t E , F 087 S , F 146 S , F 087 L , F 146 L , µ rel,N and µ rel,E have been used.In this scenario, while the precision of the lens mass decreases significantly, stellar parameters Table 1.Estimated parameters of the source and the lens by using different set of observables.Fit 1 includes constraints from θE, DS, πEN and πEE.Fit 2 uses ρ * , θS, πEN and πEE while Fit 3 uses tE, F 087S, F 146S, F 087L, F 146L, µ rel,N and µ rel,E .
of the source and the lens, as well as the absorption, are known with better precision.The main exception is the metallicity, which is not constrained for all runs.Indeed, the reconstructed distribution of the metallicity is uniform and match the distribution of the isochrones (e.g.F e = −0.8± 0.7).This is expected, not only because the constraints on this parameter are the narrow F 087 and wide F 146 bands, but also because this parameter is highly degenerate with the age (Worthey 1994).Similarly, the proper motion of the source and lens is always poorly constrained (but not their difference µ rel ).Again, this is due to the lack of information in the provided observables.

Simulation of a sample of lenses seen by Roman
We generate a realistic sample of simulated events that can be expected from the Roman Galactic Bulge Time Domain Survey (Penny et al. 2019), with the Galactic model of Koshimoto et al. (2021) 10 .We generated 1000 single source and single lens events and simulated every lightcurve of the catalog using the pyLIMA software (Bachelet et al. 2017), including parallax (Gould 2004) and assuming Gaussian noise (with a zero point of Z p = 27.4 mag and an exposure time of 50 s).To ensure minimal signal in the lightcurves, u 0 was drawn from an uniform distribution U (−1, 1) while t 0 has been forced to be observed in the ∼ 72 day windows (Penny et al. 2019).We also consider that no other blends were present in the light of sight.Then, the gradient-like method included in pyLIMA was used to fit the lightcurve (and by starting from the true solution).The estimated model parameters and their corresponding uncertainties were ultimately used as a realistic set of observables for pyLIMASS.The parameters t E , π EN , π EE were used as observables, as well as constraints from A V,S and θ S with a precision of 0.1 mag and 10%, respectively, which are typical uncertainties reported in the literature.The apparent magnitudes of the source m S and lens m L in the F087 and F146 bands were also used, with a minimum precision of 0.001 mag.These measurements will be based on the lightcurve modeling as well as the modeling of the Roman images near the event location.A double-PSF fit provides strong constraints on the source and lens magnitudes as well as the relative proper motion µ rel (Bachelet et al. 2022b), which is included as the last observable.The typical astrometric precision can be approximated by ∆X = F W HM πSN R (Lindegren 1978;Sozzetti 2005), where F W HM is the Full-Width Half Maximum of the optical system and SN R is the photometric signal-to-noise ratio.The anticipated Roman FWHM will be on the order of 1 pixel (i.e.110 mas)11 and therefore we assumed the relative precision of the relative proper motion to be equal to the less precise lens/source magnitudes components, i.e. ∆µ rel /µ rel ≈ SN R = ∆ F 146 = max(∆m S , ∆m L ) F 146 .The 1000 events were then analyzed with an MCMC exploratory run of the posterior using the 16 parameters P (after a search for the global maximum) and, for each event, we computed the mean and standard deviation of the chains for each parameters after the burn-in phase for convergence of the Markov chains.
In the Figure 3, the reconstructed mass accuracy ∆M L = M L,T rue − M L,pyLIM ASS is presented versus the accuracy on the lens distance ∆D L .The overall accuracy is excellent with, a median ∆M L = −0.007± 0.0025M ⊙ and ∆D L = 0.08 ± 0.75 kpc.However, a strong bias can be observed in the reconstructed lens mass for events where one of the components brightness (the source or the lens) is poorly known (i.e.very faint), that we quantified with ∆ F 146 .It is clear that lens masses are poorly reconstructed (both in terms of accuracy and precision) for ∆ F 146 ≥ 1 mag.This is expected because not only these magnitudes place constraints on the object masses via Equation 6, but also on µ rel (and ultimately θ E ).Moreover, a clear correlation is visible in the Figure 3.This is due to the fact that θ E can be reconstructed from the injected observables (namely t E and µ rel ) and that the lens mass can be approximated by M L ≈ θ 2 E D L /κ (Equation 4). Figure 4 shows the distribution of the accuracy, normalized accuracy, uncertainties and precision for six reconstructed parameters.The main result is that pyLIMASS can reconstruct the mass of the lenses with a median precision of σ M L M L = 20% and almost no bias.The visual absorption values towards the sources are accurately reconstructed for most of the events, albeit with a small bias.This is to be expected, given the constraint on the absorption in the observables.However, less accurate and precise values are obtained for the distance and mass of the sources.This is because the algorithm can only put weak constraints on the source properties (via θ S and the source magnitudes).And while the absorption is known accurately, the algorithm can adjust the stellar isochrones and the source distances to match the observed magnitudes.For example, if one selects isochrones within a specific color range |F 087 − F 146 − 0.5| ≤ 0.1 mag, the range of corresponding stellar parameters is wide, with M ∈ [0.6, 2] M ⊙ , T ef f ∈ [5000, 6000] K, [M/H] ∈ [−2, 0.7] dex and logg ∈ [4.8, 1.9].Therefore, we argue that the accuracy and precision observed in the reconstructed parameters are fundamental issues rather than intrinsic problems of the algorithm, given the set of observables provided.

OBSERVED EXAMPLES
In this section, we apply pyLIMASS on observed events published in the literature.Because all of these events, with the exception of Gaia16aye (Wyrzykowski et al. 2023), are located towards the Galactic Bulge, we can only make a conservative assumption about the upper limit of the source distance D S < 15 kpc (Penny et al. 2019).A summary of the results can be found in the Table 2 but every events are discussed in more details thereafter.

MOA-2009-BLG-266Lb
Muraki et al. ( 2011) reported the discovery of the anomalous microlensing event MOA-2009-BLG-266Lb.The lightcurve of the event exhibits finite-source effects and an annual parallax signal leading to the measurement of the Einstein ring radius θ E = 0.98 ± 0.04 mas and π E = (π EN , π EE ) = (0.1665 ± 0.0503, 0.1439 ± 0.0035).Ultimately, the authors found the lens system to be composed of a super-Earth mass planet m p = 10.4 ± 1.7M ⊕ orbiting a M L = 0.56 ± 0.09M ⊙ host at D L = 3.04 + 0.33 kpc.They reported the source to be a red giant with V S = 17.677 ± 0.03 mag, I S = 15.856 ± 0.030 mag and H S = 13.780± 0.030 mag.They also obtained high-resolution images with the NACO instrument on the Very Large Telescope (VLT) and the total flux at the target location is H tot = 13.755± 0.05 mag.As displayed in the Figure 5, by combining in the observables the source magnitudes (in V, I and H bands) as well as the baseline magnitude in H band (and assuming a conservative 0.1 mag errors for each measurements), the parallax vector in polar coordinates π E = (π E , ϕ E ) and ρ * = 0.00539 ± 0.00011, pyLIMASS estimates the lens system to be composed of a m p = 10.5 ± 1.9 M ⊕ planet orbiting a M L = 0.56 ± 0.12 M ⊙ host located at D L = 3.63 ± 0.55 kpc from Earth, a solution in excellent agreement with Muraki et al. (2011).However, it is relevant to point small discrepancies.First, pyLIMASS estimates the lens to be slightly more distant.This is mainly due to the fact that the source distance is not well constrained by the given observables with D S = 11.4 ± 2.2 kpc, while Muraki et al. (2011) assumed D S = 8.8 kpc.pyLIMASS also estimates the angular source radius to be θ S = 4.94 ± 0.36 µas, slightly smaller than the value reported by Muraki et al. (2011) with θ S = 5.2 ± 0.2 µas.We note however that using the color-angular radius relation of Adams et al. (2018), the angular source radius is predicted to be θ S = 4.6µas and θ S = 5.1µas (using the (V-I) and (V-H) relations of the "All Stars", respectively).Because ρ * is precisely measured, a small decrease in the angular source radius automatically leads to a smaller Einstein ring radius, that pyLIMASS estimates to be θ E = 0.919 ± 0.069 mas, and ultimately a smaller π rel since M L is precisely measured.Finally, we note that the pyLIMASS estimate of the source visual absorption A V,S = 1.89 ± 0.34 mag is lower than the measurement  of Muraki et al. (2011) A V,S = 2.14 mag, but in agreement with the OGLE Extinction Calculator12 , based on the measurements of Nataf et al. (2013), that reports A V,S = 1.84 mag.

OGLE-2005-BLG-169Lb
Originally described in Gould et al. (2006) (G06), Batista et al. (2015) (Ba15) and Bennett et al. (2015) (Be15) have reanalyzed this event using high-resolution observations from HST and Keck and they were able to measure the relative proper motion as well as the source and lens magnitudes in the B, V, I and H bands.
As observables, we used t E = 41.8 ± 2.9 days, ρ * = 0.00048 ± 0.00004, the source magnitudes B S = 23.382± 0.06 mag, V S = 22.21 ± 0.04 mag, I S = 20.555± 0.05 mag H S = 18.81 ± 0.08 mag and the lens magnitudes B l = 24.72 ± 0.15 mag, V L = 22.783 ± 0.07 mag, I L = 20.493± 0.05 mag and H L = 18.20 ± 0.10 from Ba15 and Be15.We also used µ rel,hel,N = 4.87±0.12mas/yr and µ rel,hel,E = 5.63±0.12mas/yr from Ba15.pyLIMASS estimates are consistent with the three studies, as shown in the Figure 6, with M L = 0.72 ± 0.10M ⊙ and D L = 4.94 ± 0.70 kpc.pyLIMASS estimates of θ S = 0.401 ± 0.023µas and θ E = 0.823 ± 0.040 mas are compatible but smaller compared to G06 measurements of θ S = 0.44 ± 0.04µas and θ E = 1.00 ± 0.22.However, we note, that the angular source radius is precisely known from the measurement of the relative proper motion, since: Using the µ rel,geo = 7.0 ± 0.2 mas/yr estimate from Ba15 with the ρ * and t E values from Be15, we found θ S = 0.387 ± 0.048µas, close to the pyLIMASS estimate, but different from the values of G06 and Ba15 (Be15 does not report θ S values in their Table 1).We note that pyLIMASS estimates of the Einstein ring are very similar to the one obtained by Be15, i.e. θ E = 0.848 ± 0.027 mas, when the constraints of the relative proper motion µ rel are integrated in the modeling.Again, the discrepancy between the lens distance is due to the difference in the source distance estimated by pyLIMASS with D S = 11.7 ± 1.9 kpc, while Ba15 assumed D S = 8.5 ± 1.5 kpc.

OGLE-2017-BLG-1254L
As presented in Zang et al. (2020), this event has been monitored by ground based telescopes as well as the Spitzer Space Telescope.The separation between the Earth and Spitzer at the time of the event allowed a precise measurement of the parallax vector.Despite a bimodal degeneracy, the impact on the lens mass and distance is not substantial (Zang et al. 2020).It is, however, a good test case for pyLIMASS as the two different modes can be injected simultaneously in the observables as the GM can handle multi-modal distributions, as can be seen in the Figure 7.We note that initial guesses on the position needs to be passed on to the GM to obtain good convergence of the model.Using π EN , π EE , ρ * , I S and H S from Zang et al. (2020), we obtain the lens-source distance D LS = 1.21 ± 0.35 kpc and M L = 0.62±0.07M⊙ , in agreement with the published results M L = 0.60±0.03and D LS = 0.53±0.11kpc.Again, the relative low precision on D LS is mainly due to the low precision on the source distance estimate with D S = 12.4 ± 1.8, to be compared with D S = 7.8 ± 0.8 kpc from Zang et al. (2020).Indeed, as pointed by Zang et al. (2020), the distance between the lens and the source D LS = D S (ϵ D − 1) is generally determined more precisely than D S and D L separately.As can be seen in the Figure 8, if one assumes the source distance as an observable, the mass estimate from pyLIMASS stays unchanged with M L = 0.59 ± 0.06M ⊙ but the lens-source distances are almost identical to the published values with D LS = 0.52 ± 0.10 kpc.

OGLE-2019-BLG-0960Lb
Presented by Yee et al. (2021) as the "smallest microlensing planet" (it was rather the smallest mass ratio microlensing planet at the time), the lens system is well characterized because both finite-source effets with ρ ∼ 3 × 10 −5 and the parallax vector are well measured, despite the four-fold degeneracy (Gould 2004).Yee et al. (2021) also secure three bands measurements of the source magnitudes that includes V s = 21.34 ± 0.1 mag and I s = 19.81± 0.05 mag.Coupled with their measurements of the extinction A V,S = 1.64 ± 0.1 mag, they estimated θ S = 0.625 ± 0.028µas and then θ E = 1.9 − 2.1 mas depending on the lightcurve modeling solutions.Ultimately, the four solutions are in agreement at the 1-σ level with M L ∼ 0.45M ⊙ at D L ∼ 0.8 kpc, compatible with the measured blend light I b = 18.65 ± 0.1 mag (Yee et al. 2021).We note that their source color (V − I) 0 = 0.86 mag is compatible with an early K-dwarf star Bessell & Brett (1988); Pecaut & Mamajek (2013).However, this source stellar type is in tension with their adopted source distance of D S = 7.56 kpc.Indeed, at this distance, one can anticipate a significantly smaller source angular radius of θ S ∼ 0.5µas.
We have injected the four different solutions of Yee et al. (2021) in pyLIMASS and used the observables t E , ρ * , π EN , π EE , the source magnitudes in the V S and I S bands as well as the lens magnitude in the I band.For the latter, we assumed two different distributions (and so two runs of pyLIMASS).The first run is assuming that the entire blend light is due to the lens (I L = I blend ) and the second run uses only a lower limit (I L > I blend ).Results are presented in the Figure 9.Both scenario are well compatible and found a slightly smaller Einstein ring radius θ E = 1.6 ± 0.2 mas, due to the smaller angular source radius θ S = 0.51 ± 0.05µas.Both runs agree at the 1-σ level with Yee et al. ( 2021) that the lens is an M-dwarf M L = 0.49 ± 0.09M ⊙ at D L = 1.3 ± 0.3 kpc (for I L = I blend ) and M L = 0.5 ± 0.1M ⊙ at D L = 1.3 ± 0.30 kpc (with I L ≥ I blend ), and a similar source distance of D S = 9.3 ± 2.7 kpc.  2021).(Right) The same distribution but assuming an lower limit on the lens brightness.

MOA-2013-BLG-220Lb
Originally published by Yee et al. ( 2014) (Y14 thereafter), MOA-2013-220Lb shows strong anomalies at its peaks that allows a precise measurement of the lens-planet mass ratio q = 0.003 as well as ρ * and t E .Unfortunately, only a single mass-distance relation can be extracted from the lightcurve (θ E ) and Y14 could not make any firm conclusions on the nature of the host.However, they were able to place upper limits on the lens mass M L < 0.77M ⊙ and distance D L < 6.5 kpc, based on blend flux arguments.Vandorou et al. (2020) (V20 thereafter) published a follow-up study of this event where they observed the event location with high-resolution imagers with the Keck telescope in 2015 and 2019.Their results are in contradiction with Y14 limits as they conclude the lens to be M L = 0.88 ± 0.05M ⊙ at D L = 6.72 ± 0.59 kpc.Using slightly different methods, both studies agree on the measurement of the angular source radius θ S , and finally on the measurement of θ E since ρ * is well constrained by the lightcurve.For instance, V20 found θ S = 0.689 ± 0.052 µas based on the source properties ((V − I), I) 0,S = (0.592 ± 0.038, 18.057 ± 0.054).However, this source color is in tension with their Keck source measurement in the K s band.Indeed, V20 reported K s = 17.20 ± 0.02 mag and therefore, assuming A Ks ∼ 0.2 mag based on their measurement of A V = 2.1 mag, the measured color is (I − K s ) 0 = 1.06 mag, 0.3 mag redder that the expected infrared color derived from their optical color and the tabulations of Bessell & Brett (1988) and Pecaut & Mamajek (2013).Similarly, the spectral type of such main sequence star implies R S ∼ 1.4 R ⊙ (Pecaut & Mamajek 2013), leading to a significantly different angular radius θ S = 0.795µas if the source is located at D S = 8.19 ± 0.76 kpc (V20).In fact, similarly to OGLE-2005-BLG-169Lb, the source and Einstein ring angular radii are known with great precision from the measurements of t E , ρ * and µ rel (i.e.Equation 13), leading to θ E = 0.454 ± 0.005 mas and θ S = 0.70 ± 0.01µas.As can be seen in the Figure 10, the angular radius of the source measured from the proper motion accurately matches the one predicted by V20, strongly raising the confidence in the (V − I) S color and absorption measurements.We also present in the Figure 10 the result of a Monte-Carlo experiment designed to recover the angular radius of the source based on the closest match between (V, I) S mag and the stellar isochrones of Bressan et al. (2012), that is also in good agreement with V20.However, if one assumes D S = 8.19 ± 0.76 kpc from V20, the predicted K s magnitude of the source from this sample is K s,S = 17.6 ± 0.2 mag, in tension with the K s,S = 17.20 ± 0.02 mag from V20.The cause of these tension is not clear and it could be due to incorrect calibrations, a wrong estimation of the absorption and distance to the source, a source/lens misendification or the presence of a companion to the source.The latter hypothesis seems plausible.Indeed, if one assumes the infrared magnitude of the primary to be K s = 17.6 mag, then the potential companion predicted magnitude would be K s = 18.5 mag, too faint to be detected in the V and I bands.The FWHM of 2019 Keck images presented in V20 were ∼ 60 mas, so any binary stars with a projected separation a ⊥ ≤ 60 mas × 8.19 kpc ≤ 500 AU would be included in the FWHM of the images, which is realistic.2020) and the stellar isochrones of Bressan et al. (2012), assuming the source at DS = 8.19 ± 0.76 kpc.This is in good agreement with the estimation of V20 (yellow) as well as the measurement from Equation 13. (Right) The estimated magnitude of the source in the Ks band from the isochrones is in tension in the measurement of the source (yellow) magnitude reported by V20.The lens magnitude is alsor reported in red.
It is clear nevertheless from Figure 10 that the optical colors of the source are in good agreement with the derived relative proper motion from the Keck images.Therefore, we included these information in the pyLIMASS estimation, along with the measurements of t E , ρ, A V,S and K s,L of V20, and results are presented in the Figure 11.The estimated mass agrees with, but is less precise than, V20 with M L = 0.99 ± 0.16M ⊙ , while lens distance is poorly constrained with D L = 8.7 ± 1.2 kpc, mainly because the source distance is almost unconstrained by the observables with D S = 11.3 ± 1.8kpc.The main reason is that, while the effective temperature of the source has been well measured with T ef f,S = 6200 ± 200 K, the absolute magnitudes of the source is not well constrained.Indeed, a wide range of metallicity (F e = −1.00± 0.70) and surface gravity values (logg = 3.95 ± 0.16) is possible at this location of the Hertzsprung-Russel diagram.The decrease in precision for the lens mass is also due to the fact that V20 used isochrones with a fixed age of 6.4 Gyr, while pyLIMASS considers a wider range of ages ([1-10] Gyr), allowing for more flexibility in term of mass and luminosity.Moreover, the estimated absorption towards the lens is poorly measured with A V,L = 1.16 ± 0.65 mag, mainly because only the K s band has been used, while V20 used a precise empirical relation for the lens absorption.It is relevent to note that pyLIMASS estimation is compatible with the blend flux constraint highlighted by Y14.However, as underlined by V20, the blend brightness reported by Y14 is extremely faint with I blend = 21.10 ± 0.05 mag and is potentially biased.It could be affected by an incorrect background estimation when performing crowded field photometry.While this does not affect the estimation of the source magnitudes, it could significantly bias the estimate of such a faint blend.For completness, a second run was conducted by including the additional constraint I L = I blend .From the Keck images presented in V20, it is clear that the field of view around the target is not overly crowded and only two components appear to be present and, therefore, it seems natural to assume I L = I blend .The solution leads to a closer and lighter lens with M L = 0.624 ± 0.058M ⊙ at D L = 3.96 ± 0.35 kpc.We note that in this scenario, the source distance decreases significantly to D S = 4.72 ± 0.47 kpc, and its predicted infrared magnitude is K s,S = 17.507 ± 0.065 mag.A source at this distance is potentially possible, but given the already mentionned cautions about the faint blend, this scenario seems unlikely.

Gaia16aye
Gaia16aye is a long duration event detected by the Gaia telescope (Wyrzykowski et al. 2020) (W20 thereafter).Thanks to a dense monitorig from the ground, this event has been precisely characterized.In particular, the lens orbital motion and the space parallax have been detected (W20), and the lens is a binary stars of M L,1 = 0.57 ± 0.05 and M L,2 = 0.36 ± 0.03 at 780 ± 60 pc.We used information from t E = 111.09± 0.41 d, log 10 (ρ) = −2.519± 0.003, π EN = −0.373± 0.002, π EE = −0.145± 0.001, T ef f,S = 3933 ± 133 K, logg S = 2.2 ± 1.4, F e s = −0.08 ± 0.41 as well as the magnitude of the source (V, R, I) S = (16.61± 0.02, 15.62 ± 0.02, 14.70 ± 0.02) mag.Because the lens is a binary star, we set pyLIMASS to consider the lens as non stellar (so no constraints from the isochrones are imposed).Results can be found in Figure 12 and Figure 13.With a total mass M L = 0.97 ± 0.13M ⊙ at D L = 0.759 ± 0.092 kpc, pyLIMASS results are fully consistent results with W20.Moreover, as can be seen in Figure 13, both θ S = 10 ± 1µas and θ E = 3.2 ± 0.5 mas are in good agreement with W20 that reported θ S = 9.2 ± 0.7µas and θ E = 3.04 ± 0.24 mas.Similarly, the source distance D S = 15.3 ± 3.3 kpc is in excellent agreement with W20 estimation D S = 15.7 ± 3.0 kpc.

CONCLUSIONS
We present a new algorithm to estimate the most probable physical properties of microlensing events.Based on stellar isochrones and a Gaussian Mixture approach to combine the various information collected on events, the algorithm explores the entire parameter space and recovers the most plausible scenarii.With more than thirty possible observables currently available, the user can provide as much as information available to constrain the parameters of any events.The algorithm is solely based on these observables (e.g.measurements) and the stellar isochrones, and does not assume any distribution on the parameters.In particular, it does not use any models of the Milky Way but users can easily adapt the algorithm to include such priors information.
We have demonstrated the robustness and fidelity of the algorithm on a simulated sample of planetary hosts that can be observed by the future Galactic Bulge Time Domain Survey of the Nancy Grace Roman Space Telescope.Using a set of observables that should be routinely obtained from observations, we found the algorithm to accurately and precisely reconstruct the simulated parameters as well as returning realistic uncertainties.We also run the algorithm on several events published in the literature.The algorithm generally found an agreement at a 1-σ level with the published results.The discrepancies are generally due to tensions between the measurements presented in the original studies, such ad for OGLE-2019-0960Lb or MOA-2013-BLG-220Lb, rather than the algorithm.
It is clear from Section 3.2 and Section 4 that pyLIMASS tends to systematically overestimate the source distances.This is mainly due for two reasons.First, the algorithm does not assume any models of the Milky Way.Therefore, there are no informative priors about the source distances provided (by default), which is a major difference with most of analysis presented in the literature.Indeed, most studies used galactic models priors or simply assume a source distance equal to the distance of the Red Giant Clump.A second explanation of this effect is detailed in the Section 3.2.The combination of the absorption flexibility with the diversity of stars inside the isochrones allows a vast panel of models to accurately replicate the source apparent magnitudes injected in the observables.Fortunately, the imprecise measurement of the source distance moderately impact the estimation of the lens parameters, but explains for some part the lower precision measured by pyLIMASS versus the literature seen in Section 4.
We note that the algorithm is fast to run.For example, it takes about 300 seconds to run 20000 chains with 32 walkers (with the emcee package of Foreman-Mackey et al. ( 2013)) on a personal laptop equipped with a Intel(R) Xeon(R) W-10855M CPU @ 2.80GHz and uses about 400 Mo of memory.Therefore, it is an ideal tool to estimate the parameters of microlensing events at large scale.The presented results have been obtained with a python implementation of the algorithm, named pyLIMASS, that will be officially released in the upcoming release 2.0 of the pyLIMA software (Bachelet et al. in prep.) but that can already be found online13 .
In the future, several improvements could be made to the algorithm.In particular, the inclusion of a secondary body, both for the source and the lens, would bring extra constraints.Indeed, not only the observed magnitudes would be better replicated, but also the orbital elements (and therefore the mass) of the systems could be constraints via the observation of second order effects during the events, namely the orbital of the lens (Skowron et al. 2011;Udalski et al. 2018) and the xallarap effect (Rota et al. 2021).

A. IMPLEMENTATION OF THE PARAMETERS PRIORS
As usual with Bayesian approaches, the choice of the parameters priors p(P) is critical (Tak et al. 2018).While the priors distributions can be changed by the user, we described in this section the default options.
As discussed previously, the microlensing system is described by sixteen parameters of interest P .The more trivial, but not innocent (Bailer-Jones 2015; Eadie et al. 2023), choice is to use a uniform prior.This is a valid approach, as long as lower and upper limits are finite and the range not too narrow (Tak et al. 2018), and are often referred as weakly informative prior.But in the present case, such priors on stellar parameters could create non-astrophysical objects and generate non-physical observables.To address this, pyLIMASS incorporates a prior probability distribution on the stellar parameters by comparing the stellar parameters M , T ef f , F e and logg with the stellar isochrones.In practice, the squared Euclidean distance δ 2 is computed between the drawn parameters and all of the isochrones catalog, and the minimum distance is selected.The distance is then turned into a probability p(P ) via: where σ is a relaxing term that is set to σ ∼ 0.05.We note that this needs to be systematically applied to the source, but not necessarily to the lens object.Indeed, the lens could be an object that is not present in the catalog of stellar isochrones, such as a free-floating planet, a brown dwarf or a stellar remnant.The isochrone constraints for the lens are not set by default, but the user can force the lens to be a stellar object by activating this constraint.If not, then uniform prior on the stellar parameters are used.Similarly, uniform priors for the distances D, proper motions µ, visual absorption A v,s are used and a summary can be found in the Table 3.This is a similar approach as Mróz & Wyrzykowski (2021) and Kruszyńska et al. (2022).U(0,1) Table 3. List of default priors used in pyLIMASS.There are tow possibilities for the lens, depending if the lens is forced to be stellar (see text).U(a,b) indicates a uniform distribution of the prior between a and b.

B. LIST OF PYLIMASS DEFAULT OBSERVABLES
The following list shows the default observables that can be included in pyLIMASS.There is no mandatory selection of observables, but at least one observable needs to be provided by the user.We note that in cases where the heliocentric proper motions are provided as observable, the reference time t 0,par as well as the equatorial event coordinates (α and δ) must be passed to the algorithm to calculate the geocentric to heliocentric transformation (Skowron et al. 2011).

• Lens apparent magnitudes
• Baseline apparent magnitudes tutions, in particular the institutions participating in the Gaia Multilateral Agreement.This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France.

Figure 1 .
Figure 1.(Left) Simulated V versus U distribution (top left) and the different reconstructed GM distributions with different number of Gaussian components.As more components are added to the model, the simulated distribution is reconstructed with more fidelity.(Right) BIC (plain blue) and AIC (dashed red) from a validation sample versus the number of Gaussian component used in the GM model.Both metrics reaches a plateau after ∼ 20 components, indicating potential overfitting afterwards.

Figure 2 .
Figure 2. Summary of the proposed algorithm.The first step is to collect the observables and inject them in the GM object.Then, the observables distribution are fitted by the GM model to found the best hyperparameters θ f .Finally, the likelihood function of the GM model is used in conjunction with user defined priors, stellar isochrones and microlensing parameters to explore the poterior probability p(P |X).

Figure 3 .
Figure 3. (Left)Lens mass accuracy distribution versus lens distance accuracy distribution for the sample of the Roman events presented in the text.The various levels of colours and lines indicates the 39%, 86%, 98.9% and 99.9% density contours from pyLIMASS while the 1-d distribution are also presented.The red cross indicated the median uncertainty, in both direction, reported by pyLIMASS.(Right) Distribution of the mass accuracy versus the maximum magnitude uncertainty between the lens and the source in the F146 band ∆F146.

Figure 4 .
Figure 4. Various distributions for six parameters of prime interest of the Roman sample.From top to bottom are displayed the distribution of the parameters accuracy, normalized accuracy, uncertainties and precision.The displayed range of parameters have been limited to the [10%,90%] percentiles for plotting purposes.

Figure 5 .
Figure 5. Lens mass versus lens distance for MOA-2009-266Lb (Muraki et al. 2011).The various levels of colours and lines indicates the 39%, 86%, 98.9% and 99.9% density contours from pyLIMASS while the measurement ofMuraki et al. (2011) is displayed in red.The contours have been smoothed with a Gaussian kernel for plotting purposes.

Figure 6 .
Figure 6.Same as Figure 5 but for the microlensing event OGLE-2005-169Lb.The measurements of G06, Ba15 and Be15 are displayed in black, orange and red respectively.

Figure 7 .
Figure 7. Distribution of the parallax components for the two models of OGLE-2017-1254L described in Zang et al. (2020) (blue stars) and samples generated by the GM model (orange dots).

Figure 8 .
Figure 8. Distribution of the lens mass and the lens-source distance for the microlensing event OGLE-2017-1254L (Zang et al. 2020).(Left)The pyLIMASS distribution with weak constraints on the source distance (i.e. 2 < DS < 15kpc) and in red the measurements from Zang et al. (2020).(Right) The same distribution but assuming the source distance DS = 7.8 ± 0.8 kpc as an observable.

Figure 9 .
Figure 9. Distribution of the lens mass and the lens distance for the microlensing event OGLE-2019-0960Lb (Yee et al. 2021).(Left)The pyLIMASS distribution with strong constraint on the lens flux and in red the measurements from Yee et al. (2021).(Right) The same distribution but assuming an lower limit on the lens brightness.

Figure 10 .
Figure 10.(Left) Reconstruced angular source radius of the source of MOA-2013-220Lb from the (V, I)S measurements of Vandorou et al. (2020) and the stellar isochrones of Bressan et al. (2012), assuming the source at DS = 8.19 ± 0.76 kpc.This is in good agreement with the estimation of V20 (yellow) as well as the measurement from Equation13.(Right) The estimated magnitude of the source in the Ks band from the isochrones is in tension in the measurement of the source (yellow) magnitude reported by V20.The lens magnitude is alsor reported in red.

Figure 11 .
Figure 11.Distribution of the lens mass and the lens distance for the microlensing event MOA-2013-BLG-220Lb (Yee et al. 2014; Vandorou et al. 2020).(Left)The pyLIMASS distribution with the lens flux in Ks bands and in red the measurements from V20. (Right) The same distribution but assuming the extra constraint on the lens flux IL = 21.10 ± 0.1 mag.

Figure 12 .
Figure12.Distribution of the lens mass and the lens distance for the microlensing event Gaia16aye(Wyrzykowski et al. 2020).The measurements of W20 are represented in red.

Figure 13 .
Figure 13.pyLIMASS posterior distribution of the angular source radius of the source (Left), Einstein angular radii (Middle) and source distance (Right) for the event Gaia16aye (Wyrzykowski et al. 2020).The 1σ measurements of W20 are represented in red.

Table 2
. pyLIMASS estimated lenses mass ML, lenses distance DL and sources distance for several events published in the literature.For the event MOA-2009-BLG-266Lb, (VS, IS, HS), H baseline , πE, ϕE, ρ * , AV,S have been used as observables.