# Microlensing Parallax Fitting Tutorial¶

Here you will learn how to fit the observed light-curve with a model that has a point lens, point source, but their relative motion is not rectilinear and includes annual microlensing parallax effect. If you haven’t yet looked at the other tutorials, then you may want to start there and come back here later.

We need some data to fit. Let’s look at some long-timescale event from Wyrzykowski et al. (2015). As an example we can take event with ID 3291 (star BLG234.6.I.218982; also named OGLE-2005-BLG-086 on OGLE EWS website). We can look at the plot of model without parallax: here. We see that model does not fit data very well - there are trends in the residuals around 3400 and 3800.

## Imports and settings¶

We start by importing python modules:

```
import os
import numpy as np
import emcee
import matplotlib.pyplot as plt
import MulensModel as mm
```

We’re using EMCEE package for fitting. You can download it from http://dfm.io/emcee/current/ in case you don’t have it yet. Then we import the data and set the event coordinates:

```
file_name = os.path.join(mm.DATA_PATH,
"photometry_files", "OB05086", "starBLG234.6.I.218982.dat")
my_data = mm.MulensData(file_name=file_name, add_2450000=True)
coords = "18:04:45.71 -26:59:15.2"
```

Note that *add_2450000=True* is very important. The file has time vector
of HJD-2450000 for convenience and it’s fine as long as we don’t fit
the annual parallax. He have to set *MulensData* time vector to
full HJD because this is used to calculate the position of Earth
relative to Sun and we don’t want this to be off by around 6700 years.

We set the starting values of the parameters:

```
params = dict()
params['t_0'] = 2453628.3
params['t_0_par'] = 2453628.
params['u_0'] = 0.37
params['t_E'] = 100.
params['pi_E_N'] = 0.
params['pi_E_E'] = 0.
my_model = mm.Model(params, coords=coords)
my_event = mm.Event(datasets=my_data, model=my_model)
```

We set the parameter reference time (*t_0_par*) for rounded value of *t_0*.
This is common approach. If you don’t set *t_0_par*, then fitting will be
slower, because Earth positions will be re-calculated for every model.

Further we need to specifies which parameters we want to fit and also specify dispersions in starting points. We choose the latter to be some small values:

```
parameters_to_fit = ["t_0", "u_0", "t_E", "pi_E_N", "pi_E_E"]
sigmas = [0.01, 0.001, 0.1, 0.01, 0.01]
```

Some more EMCEE settings - number of walkers, steps, and burn-in steps. Also the list of starting points:

```
n_dim = len(parameters_to_fit)
n_walkers = 40
n_steps = 500
n_burn = 150
start_1 = [params[p] for p in parameters_to_fit]
start = [start_1 + np.random.randn(n_dim) * sigmas
for i in range(n_walkers)]
```

We need one more important piece of information - the function that computes the logarithm of (unnormalized) probability. We split it into three separate functions for clarity. In the first function we also remember the smallest chi2 nad corresponding parameters:

```
def ln_like(theta, event, parameters_to_fit):
""" likelihood function """
for (parameter, value) in zip(parameters_to_fit, theta):
setattr(event.model.parameters, parameter, value)
chi2 = event.get_chi2()
if chi2 < ln_like.chi2_min:
ln_like.chi2_min = chi2
ln_like.chi2_min_theta = theta
return -0.5 * chi2
ln_like.chi2_min = np.inf
```

```
def ln_prior(theta, parameters_to_fit):
"""priors - we only reject obviously wrong models"""
if theta[parameters_to_fit.index("t_E")] < 0.:
return -np.inf
return 0.0
```

```
def ln_prob(theta, event, parameters_to_fit):
""" combines likelihood and priors"""
ln_prior_ = ln_prior(theta, parameters_to_fit)
if not np.isfinite(ln_prior_):
return -np.inf
ln_like_ = ln_like(theta, event, parameters_to_fit)
if np.isnan(ln_like_):
return -np.inf
return ln_prior_ + ln_like_
```

## Running the sampler¶

Ok, we’re ready to run EMCEE:

```
sampler = emcee.EnsembleSampler(
n_walkers, n_dim, ln_prob, args=(my_event, parameters_to_fit))
sampler.run_mcmc(start, n_steps)
samples = sampler.chain[:, n_burn:, :].reshape((-1, n_dim))
```

And now we’re ready to look at the results and best-fitted model:

```
results = np.percentile(samples, [16, 50, 84], axis=0)
print("Fitted parameters:")
form = "{:.5f} {:.5f} {:.5f}"
for i in range(n_dim):
r = results[1, i]
print(form.format(r, results[2, i]-r, r-results[0, i]))
print("\nBest model:")
print(*list(ln_like.chi2_min_theta))
print(ln_like.chi2_min)
```

I hope you got (u_0, t_E, pi_E_N, pi_E_E) of around (0.44, 95, 0.21, 0.10) and chi^2 of 949.5.

At this point you may want to say that the fit is done at this point. But it’s not! We have to check for degenerate solution. We’re fitting single lens model, hence, the search for degenerate solution is easy and it’s enough to start with negative u_0.

Now you have time to do the second fit…

Ok, I hope you got (u_0, t_E, pi_E_N, pi_E_E) of (-0.41, 110, -0.30, 0.11) and chi^2 of 947.0. The difference between the two solutions is small in chi^2 - they are degenerate. And u_0<0 fits data slightly better. It turned out that the second fit was very important!

## Plotting¶

Let’s make a nice plot!

I provide model parameters below. Here is how it goes:

```
plt.figure()
model_0 = mm.Model({'t_0': 2453628.29062, 'u_0': 0.37263,
't_E': 102.387105})
model_1 = mm.Model({'t_0': 2453630.35507, 'u_0': 0.488817,
't_E': 93.611301, 'pi_E_N': 0.2719, 'pi_E_E': 0.1025,
't_0_par': params['t_0_par']}, coords=coords)
model_2 = mm.Model({'t_0': 2453630.67778, 'u_0': -0.415677,
't_E': 110.120755, 'pi_E_N': -0.2972, 'pi_E_E': 0.1103,
't_0_par': params['t_0_par']}, coords=coords)
event_0 = mm.Event(datasets=my_data, model=model_0)
event_1 = mm.Event(datasets=my_data, model=model_1)
event_2 = mm.Event(datasets=my_data, model=model_2)
t_1 = 2453200.
t_2 = 2453950.
plot_params = {'lw': 2.5, 'alpha': 0.3, 'subtract_2450000': True,
't_start': t_1, 't_stop': t_2}
my_event.plot_data(subtract_2450000=True)
event_0.plot_model(label='no pi_E', **plot_params)
event_1.plot_model(label='pi_E, u_0>0', **plot_params)
event_2.plot_model(label='pi_E, u_0<0', color='black', ls='dashed',
**plot_params)
plt.xlim(t_1-2450000., t_2-2450000.)
plt.legend(loc='best')
plt.title('Data and 3 fitted models')
plt.show()
```

I hope you see that parallax models are better than the non-parallax model. If not, then zoom-in around epoch 3800. The non-parallax model has chi^2 higher by about 400.

Slightly modified source code from this tutorial is example 6. Additionally, example 7 shows how to fit parallax model using MultiNest instead of EMCEE algorithm. Note that a single run of MultiNest finds two degenerate modes and reports properties of both of them.

## Exercise¶

As an exercise you may try to fit other events from Wyrzykowski et al. (2015). It’s best to start with long events, that have bright sources, and small impact parameters.