ebachelet / pyLIMA

GNU General Public License v3.0
31 stars 8 forks source link

PSPL fit with parallax not returning correct values #46

Closed caseylam closed 1 year ago

caseylam commented 2 years ago

I'm trying to fit a microlensing event with parallax. To validate my procedure, I'm testing on some fake data that I generated (so I know the input parameters). Using pyLIMA, I'm not able to get the correct values back from the fit. For example, the blending parameter g I input is 0.25, but it returns -0.05059 with an error 0.01565.

I am setting up the fitter like this (following the example jupyter notebook).

telescope_1 = telescopes.Telescope(name='OGLE', camera_filter='I',
                                   light_curve_magnitude=data_1)

your_event.telescopes.append(telescope_1)

your_event.check_event()

model_1 = microlmodels.create_model('PSPL', your_event)

your_event.fit(model_1,'LM')

model_2 = microlmodels.create_model('PSPL', your_event, parallax=['Full', 57000])

model_2.parameters_guess = your_event.fits[0].fit_results[:3]+[0,0]

your_event.fit(model_2,'LM')

your_event.fit(model_2,'DE',DE_population_size=5)

your_event.fits[-1].produce_outputs()
print ('Chi2_LM :',your_event.fits[-1].outputs.fit_parameters.chichi)

I've also tried parallax=Annual, which gives the same thing. I've also tried leaving out the DE fit, which results in a similar answer but with even smaller errors.

Might you have any suggestions on what I might be doing incorrectly? Thanks!

ebachelet commented 2 years ago

For parallax fit (and also in general), the time unit is important and need to be JD (or HJD), both in the data injected in the telescopes, but as well in the t0,par argument:

model_2 = microlmodels.create_model('PSPL', your_event, parallax=['Full', 57000+2400000.5]) (assuming your time is MJD) I would give it a try with this to see how it change the results.

caseylam commented 2 years ago

Thank you! Yes, the bulk of the problem was the MJD vs (H)JD issue.

The value it is returning is much closer now to what I put in. However, there's still some small residuals left in the fit and it seems like the exact solution is still not being found (or the uncertainties are too small): g is now being found as 0.22277 with an error 0.00051.

Is this just as well as can be expected to be done? Or do uncertainties tend to be underestimated?

Or is the discrepancy due to another mistake I made somewhere?

ebachelet commented 2 years ago

It is hard to know without going more deeply in the data. However, gradient-like methods sometimes return unreliable covariance matrices. I would recommend you to try the fit with MCMC method to see if g errors are better estimated.

rpoleski commented 2 years ago

I'd suggest to simulate the data with very small uncertainties and fit them. You should get parameters that are very close to the assumed values.