Notes on fits

Below are some notes, details and advices about the fits. Most of the fits methods relies on scipy. While several algorithms are implemented, see here, pyLIMA relies mainly of three main fitting methods.

Differential Evolution (DE)

DE is a population-based algorithm for global optimization of complex functions by Storn & Price. There are plenty of litterature discussing the remarkable performance about DE and its derivatives.

DE is the workhorse of pyLIMA fits and it has proven its reliability and effectiveness on many fits.

Gradient-Like methods

There are two methods to performs gradient-like fits in pyLIMA, the Trust-Reflective Function (TRF) and Levenberg-Marquardt (LM). They are almost identical, but the former accounts of parameters boundaries (which is desirable when the event is not very well constrained). They are very efficient to find the best models as soon as a minima is found. Jacobian are implemented for simplest models (i.e. PSPL and FSPL without second-order effects).

MCMC

MCMC is implemeented via the awesome emcee package. The number of walkers and links can be adjusted. Uniform priors at the parameters boundaries are set by default.

Fancy parameters

Modeling performance can be significantly improved by exploring the parameter space in log-space. This is especially recommended for the modeling of binary lenses. While users can pass their own fancy parameters, pyLIMA includes an option to use \(log_{10}(t_E)\), \(log_{10}(\rho)\), \(log_{10}(s)\) and \(log_{10}(q)\):

from pyLIMA.models import USBL_model, pyLIMA_fancy_parameters
from pyLIMA.fits import DE_fit

fancy_parameters = pyLIMA_fancy_parameters.standard_fancy_parameters

usbl = USBL_model.USBLmodel(current_event, fancy_parameters=fancy_parameters)

de = DE_fit.DEfit(usbl)
de.fit()

The boundaries are automatically adjusted to the new parameters. However, the parameters guess (for gradient-like and MCMC methods) needs to be given in the new scale:

from pyLIMA.models import USBL_model, pyLIMA_fancy_parameters
from pyLIMA.fits import TRF_fit

fancy_parameters = pyLIMA_fancy_parameters.standard_fancy_parameters.copy()

usbl = USBL_model.USBLmodel(current_event, fancy_parameters=fancy_parameters)

trf = TRF_fit.TRFfit(usbl)
# t0, u0, log10(tE),log10(rho),log10(s),log10(q),alpha
trf.model_parameters_guess = [ 2.45986449e+06,  -1.81080674e-01,  2.63277601e+00, -1.45312165e+00, -1.85934608e-01,-3.12504456e+00,  5.42587262e+00]
trf.fit()

Parallelization

DE and MCMC methods can significantly be speed-up by implementing a pool of workers via multiprocessing:

import multiprocessing as mul
pool = mul.Pool(processes = 4)

my_fit.fit(computational_pool = pool)

Priors

pyLIMA now includes the possibility to add user-defined priors. While they are no priors by default, uniform and gaussian priors are available. Users can also define their own functions as long as they return a pdf for a given parameters as well as a rvs method, for example with a Cauchy distribution:

class CauchyDistribution(object):

    def __init__(self, mean, gamma):

        self.mean = mean
        self.gamma = gamma

    def pdf(self, x):

        denominator = np.pi*self.ggam*(1+(x-self.mean)**2/self.gamma**2)
        probability = 1 / denominator

        return probability

    def rvs(self, size):

        sample = np.random.standard_cauchy(size)

        samples = self.mean+self.gamma*sample

        return samples

from pyLIMA.models import PSPLmodel
from pyLIMA.fits import DEfit

model = PSPLmodel(event)
thefit = DEfit(model)

t0prior =  CauchyDistribution(2459856,0.5)
u0prior =  CauchyDistribution(0.1,0.5)
tEprior =  CauchyDistribution(22,0.5)

thefit.priors = [t0prior,u0prior,tEprior]

Loss functions

By default, pyLIMA implements three loss functions:

\(\chi^2\) : the sum of the normed residuals
\(\log \cal L\) : the ln-likelihood, that includes priors
soft_l1 : the soft_l1 function is close to the Huber loss function and it is very robust against outliers

from pyLIMA.fits import DEfit
thefit = DEfit(model,loss_function='soft_l1')

Fitting algorithms have default loss functions described in pyLIMA Modules. The sign of the loss function will depends if the fitting algorithms maximize or minimize the objective function.

Advices on fitting binary lightcurves

For fitting binary models, DE has proven to be reliable to locate global minima. However, we recommand to explore \(s\le1\) and \(s\ge1\) separetely, to explore carefully the close/wide degeneracy (see). One the minimas are found, each of them should be explored using MCMC.

We note that some wide binary systems can be hard, if not impossible, to model with the default pyLIMA settings. OGLE-2015-BLG-0060 is a good example. In this case, it is recomanded to change the origin of the system, for example to the primary body:

from pyLIMA.models import USBLmodel

usbl = USBLmodel(current_event,origin=['primary',[0,0]])