From 06672c31e08fd26c81da28edc2219e4266cbae44 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 21 May 2025 21:23:32 +0200 Subject: [PATCH 1/5] Added calc_chisq method to fix issue with calculated chi2 --- backtracks/backtracks.py | 41 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index ed59471..9ab534c 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -475,6 +475,43 @@ def fmodel(self, param): """ return radecdists(self, self.epochs_tt, param) + def calc_chisq(self, param): + """ + Function to calculate the chi2 for correlated (e.g., GRAVITY) and uncorrelated datapoints given certain model parameters. + + Args: + param (np.array of float): set of parameters describing 5D position of background star and 6D position of host star. + + Returns: + Total chi2 given the observed datapoints and a set of model parameters. + """ + + if not np.isfinite(np.sum(param)): + # if nan or inf, avoid calling fmodel + return -np.inf + + xs, ys = self.fmodel(param) + + # separate terms where there is a correlation + corr_terms = ~np.isnan(self.rho) + + chi_sq = 0.0 + + if np.sum(~corr_terms) > 0: + # calculate chi2 for uncorrelated terms + chi_sq += np.sum((self.ras[~corr_terms] - xs[~corr_terms])**2 / self.raserr[~corr_terms]**2 + chi_sq += np.sum((self.decs[~corr_terms] - ys[~corr_terms])**2 / self.decserr[~corr_terms]**2 + + if np.sum(corr_terms) > 0: + # calculate the chi2 for correlated terms + # following equations with correlation terms taken from orbitize! routine + # _chi2_2x2cov (https://orbitize.readthedocs.io/en/latest/modules/orbitize/lnlike.html) + det_C = (self.raserr[corr_terms]**2)*(self.decserr[corr_terms]**2)*(1-self.rho[corr_terms]**2) + covs = self.rho[corr_terms]*self.raserr[corr_terms]*self.decserr[corr_terms] + chi_sq += ((self.ras[corr_terms]-xs[corr_terms])**2*self.decserr[corr_terms]**2+(self.decs[corr_terms]-ys[corr_terms])**2*self.raserr[corr_terms]**2-2*(self.ras[corr_terms]-xs[corr_terms])*(self.decs[corr_terms]-ys[corr_terms])*covs)/det_C + + return chi_sq + def loglike(self, param): """ Function to calculate Log likelihood for correlated (e.g., GRAVITY) and uncorrelated datapoints given certain model parameters. @@ -725,7 +762,7 @@ def fit(self, dlogz=0.5, npool=4, dynamic=False, nlive=200, mpi_pool=False, resu # Compute useful chi2 value self.median_loglike = self.loglike(self.run_median) - self.median_chi2_red = -2.*self.median_loglike/((2*len(self.epochs))-self.ndim) + self.median_chi2_red = self.calc_chisq(self.run_median)/((2*len(self.epochs))-self.ndim) return self.results @@ -762,7 +799,7 @@ def load_results(self, fileprefix: str = './'): # recompute useful chi2 value self.median_loglike = self.loglike(self.run_median) - self.median_chi2_red = -2.*self.median_loglike/((2*len(self.epochs))-self.ndim) + self.median_chi2_red = self.calc_chi2(self.run_median)/((2*len(self.epochs))-self.ndim) def generate_plots( self, From d8a71146290a8c57aef3e9711acf40c215723d86 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 21 May 2025 21:39:10 +0200 Subject: [PATCH 2/5] Fixed typo --- backtracks/backtracks.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index 9ab534c..13001e4 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -499,8 +499,8 @@ def calc_chisq(self, param): if np.sum(~corr_terms) > 0: # calculate chi2 for uncorrelated terms - chi_sq += np.sum((self.ras[~corr_terms] - xs[~corr_terms])**2 / self.raserr[~corr_terms]**2 - chi_sq += np.sum((self.decs[~corr_terms] - ys[~corr_terms])**2 / self.decserr[~corr_terms]**2 + chi_sq += np.sum((self.ras[~corr_terms] - xs[~corr_terms])**2 / self.raserr[~corr_terms]**2) + chi_sq += np.sum((self.decs[~corr_terms] - ys[~corr_terms])**2 / self.decserr[~corr_terms]**2) if np.sum(corr_terms) > 0: # calculate the chi2 for correlated terms From bd765751fd0d999afa82ac944ef80b69bd89262c Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 21 May 2025 21:53:18 +0200 Subject: [PATCH 3/5] Another minor fix --- backtracks/backtracks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index 13001e4..fddb8af 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -799,7 +799,7 @@ def load_results(self, fileprefix: str = './'): # recompute useful chi2 value self.median_loglike = self.loglike(self.run_median) - self.median_chi2_red = self.calc_chi2(self.run_median)/((2*len(self.epochs))-self.ndim) + self.median_chi2_red = self.calc_chisq(self.run_median)/((2*len(self.epochs))-self.ndim) def generate_plots( self, From 838372a03f17490c2d12e21be4e8ea180667f582 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 21 May 2025 22:18:57 +0200 Subject: [PATCH 4/5] Another fix --- backtracks/backtracks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index fddb8af..2d0eaad 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -508,7 +508,7 @@ def calc_chisq(self, param): # _chi2_2x2cov (https://orbitize.readthedocs.io/en/latest/modules/orbitize/lnlike.html) det_C = (self.raserr[corr_terms]**2)*(self.decserr[corr_terms]**2)*(1-self.rho[corr_terms]**2) covs = self.rho[corr_terms]*self.raserr[corr_terms]*self.decserr[corr_terms] - chi_sq += ((self.ras[corr_terms]-xs[corr_terms])**2*self.decserr[corr_terms]**2+(self.decs[corr_terms]-ys[corr_terms])**2*self.raserr[corr_terms]**2-2*(self.ras[corr_terms]-xs[corr_terms])*(self.decs[corr_terms]-ys[corr_terms])*covs)/det_C + chi_sq += np.sum((self.ras[corr_terms]-xs[corr_terms])**2*self.decserr[corr_terms]**2+(self.decs[corr_terms]-ys[corr_terms])**2*self.raserr[corr_terms]**2-2*(self.ras[corr_terms]-xs[corr_terms])*(self.decs[corr_terms]-ys[corr_terms])*covs)/det_C) return chi_sq From 7606b73a09b2c7bd40527b0d312a39ee6c578b48 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 21 May 2025 22:23:44 +0200 Subject: [PATCH 5/5] More fixing --- backtracks/backtracks.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/backtracks/backtracks.py b/backtracks/backtracks.py index 2d0eaad..098f0ac 100644 --- a/backtracks/backtracks.py +++ b/backtracks/backtracks.py @@ -508,7 +508,7 @@ def calc_chisq(self, param): # _chi2_2x2cov (https://orbitize.readthedocs.io/en/latest/modules/orbitize/lnlike.html) det_C = (self.raserr[corr_terms]**2)*(self.decserr[corr_terms]**2)*(1-self.rho[corr_terms]**2) covs = self.rho[corr_terms]*self.raserr[corr_terms]*self.decserr[corr_terms] - chi_sq += np.sum((self.ras[corr_terms]-xs[corr_terms])**2*self.decserr[corr_terms]**2+(self.decs[corr_terms]-ys[corr_terms])**2*self.raserr[corr_terms]**2-2*(self.ras[corr_terms]-xs[corr_terms])*(self.decs[corr_terms]-ys[corr_terms])*covs)/det_C) + chi_sq += np.sum(((self.ras[corr_terms]-xs[corr_terms])**2*self.decserr[corr_terms]**2+(self.decs[corr_terms]-ys[corr_terms])**2*self.raserr[corr_terms]**2-2*(self.ras[corr_terms]-xs[corr_terms])*(self.decs[corr_terms]-ys[corr_terms])*covs)/det_C) return chi_sq