From d45bd669de5749603017975221ba4e053ef2c7e5 Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 12:14:55 -0500 Subject: [PATCH 1/8] Do not include negative values in the Gaussian prior This was causing unphysical (negative) stellar masses to be returned. Fixed. --- orbitize/priors.py | 21 +++++++++++++++++++-- orbitize/results.py | 1 + 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/orbitize/priors.py b/orbitize/priors.py index 3c8c6fa2..754fee93 100644 --- a/orbitize/priors.py +++ b/orbitize/priors.py @@ -50,7 +50,8 @@ def __repr__(self): def draw_samples(self, num_samples): """ - Draw samples from a Gaussian distribution. + Draw positive samples from a Gaussian distribution. + Negative samples will not be returned. Args: num_samples (float): the number of samples to generate @@ -59,14 +60,27 @@ def draw_samples(self, num_samples): numpy array of float: samples drawn from the appropriate Gaussian distribution. Array has length `num_samples`. """ + samples = np.random.normal( loc=self.mu, scale=self.sigma, size=num_samples - ) + ) + bad = np.inf + + while bad != 0: + + bad_samples = np.where(samples <= 0)[0] + bad = len(bad_samples) + + samples[bad_samples] = np.random.normal( + loc=self.mu, scale=self.sigma, size=bad + ) + return samples def compute_lnprob(self, element_array): """ Compute log(probability) of an array of numbers wrt a Gaussian distibution. + Negative numbers return a probability of -inf. Args: element_array (float or np.array of float): array of numbers. We want the @@ -80,6 +94,9 @@ def compute_lnprob(self, element_array): """ lnprob = -0.5*np.log(2.*np.pi*self.sigma) - 0.5*((element_array - self.mu) / self.sigma)**2 + bad_samples = np.where(samples <= 0)[0] + lnprob[bad_samples] = -np.inf + return lnprob class JeffreysPrior(Prior): diff --git a/orbitize/results.py b/orbitize/results.py index 3304fee1..36669d76 100644 --- a/orbitize/results.py +++ b/orbitize/results.py @@ -360,6 +360,7 @@ def plot_orbits(self, parallax=None, total_mass=None, object_mass=0, # Create a linearly increasing colormap for our range of epochs norm = mpl.colors.Normalize(vmin=np.min(epochs), vmax=np.max(epochs[-1,:])) + norm_yr = mpl.colors.Normalize( vmin=np.min(Time(epochs,format='mjd').decimalyear), vmax=np.max(Time(epochs,format='mjd').decimalyear) From 147497f0f237ca027fe11a8fb3673572c7502a02 Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 12:22:23 -0500 Subject: [PATCH 2/8] Fixed a variable name type-o --- orbitize/priors.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbitize/priors.py b/orbitize/priors.py index 754fee93..09f04ca5 100644 --- a/orbitize/priors.py +++ b/orbitize/priors.py @@ -94,10 +94,10 @@ def compute_lnprob(self, element_array): """ lnprob = -0.5*np.log(2.*np.pi*self.sigma) - 0.5*((element_array - self.mu) / self.sigma)**2 - bad_samples = np.where(samples <= 0)[0] + bad_samples = np.where(element_array <= 0)[0] lnprob[bad_samples] = -np.inf - return lnprob + return lnprob class JeffreysPrior(Prior): """ From 853483ff72e98bed2e8e723b36f4f2f0402ad9c9 Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 12:42:17 -0500 Subject: [PATCH 3/8] Fixed unit tests for Gaussian prior --- tests/test_priors.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/test_priors.py b/tests/test_priors.py index 5939af15..9807f8ab 100644 --- a/tests/test_priors.py +++ b/tests/test_priors.py @@ -7,7 +7,7 @@ threshold = 1e-1 initialization_inputs = { - priors.GaussianPrior : [0., 1.], + priors.GaussianPrior : [100., 1.], priors.JeffreysPrior : [1., 2.], priors.UniformPrior : [0., 1.], priors.SinPrior : [], @@ -15,7 +15,7 @@ } expected_means_mins_maxes = { - priors.GaussianPrior : (0.,-np.inf,np.inf), + priors.GaussianPrior : (100.,0.,np.inf), priors.JeffreysPrior : (1/np.log(2),1., 2.), priors.UniformPrior : (0.5, 0., 1.), priors.SinPrior : (np.pi/2., 0., np.pi), @@ -23,7 +23,7 @@ } lnprob_inputs = { - priors.GaussianPrior : np.array([0., -np.inf, np.inf, 3., 3.]), + priors.GaussianPrior : np.array([-3.0, np.inf, 100., 99.]), priors.JeffreysPrior : np.array([-1., 0., 1., 1.5, 2., 2.5]), priors.UniformPrior : np.array([0., 0.5, 1., -1., 2.]), priors.SinPrior : np.array([0., np.pi/2., np.pi, 10., -1.]), @@ -31,7 +31,7 @@ } expected_probs = { - priors.GaussianPrior : np.array([nm(0,1).pdf(0), 0., 0., nm(0,1).pdf(3), nm(0,1).pdf(-3)]), + priors.GaussianPrior : np.array([0., 0., nm(100.,1.).pdf(100.), nm(100.,1.).pdf(99.)]), priors.JeffreysPrior : np.array([0., 0., 1., 2./3., 0.5, 0.])/np.log(2), priors.UniformPrior : np.array([1., 1., 1., 0., 0.]), priors.SinPrior : np.array([0., 0.5, 0., 0., 0.]), @@ -65,6 +65,7 @@ def test_compute_lnprob(): values2test = lnprob_inputs[Prior] lnprobs = TestPrior.compute_lnprob(values2test) + assert np.log(expected_probs[Prior]) == pytest.approx(lnprobs, abs=threshold) From 5cbf0af927c30154c810dbfb3c1ff97345450c94 Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 18:34:56 -0500 Subject: [PATCH 4/8] Add keyword to control whether `GaussianPrior` returns negative samples --- orbitize/priors.py | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/orbitize/priors.py b/orbitize/priors.py index 09f04ca5..453a4ffe 100644 --- a/orbitize/priors.py +++ b/orbitize/priors.py @@ -48,13 +48,15 @@ def __init__(self, mu, sigma): def __repr__(self): return "Gaussian" - def draw_samples(self, num_samples): + def draw_samples(self, num_samples, no_negatives=True): """ Draw positive samples from a Gaussian distribution. Negative samples will not be returned. Args: num_samples (float): the number of samples to generate + no_negatives (bool): if True, only positive samples will + be drawn (default: False). Returns: numpy array of float: samples drawn from the appropriate @@ -66,18 +68,20 @@ def draw_samples(self, num_samples): ) bad = np.inf - while bad != 0: + if no_negatives: - bad_samples = np.where(samples <= 0)[0] - bad = len(bad_samples) + while bad != 0: - samples[bad_samples] = np.random.normal( - loc=self.mu, scale=self.sigma, size=bad - ) + bad_samples = np.where(samples <= 0)[0] + bad = len(bad_samples) + + samples[bad_samples] = np.random.normal( + loc=self.mu, scale=self.sigma, size=bad + ) return samples - def compute_lnprob(self, element_array): + def compute_lnprob(self, element_array, no_negatives=True): """ Compute log(probability) of an array of numbers wrt a Gaussian distibution. Negative numbers return a probability of -inf. @@ -86,6 +90,8 @@ def compute_lnprob(self, element_array): element_array (float or np.array of float): array of numbers. We want the probability of drawing each of these from the appopriate Gaussian distribution + no_negatives (bool): if True, the log(probability) of negative numbers + will be -inf (default:True). Returns: numpy array of float: array of log(probability) values, @@ -94,10 +100,12 @@ def compute_lnprob(self, element_array): """ lnprob = -0.5*np.log(2.*np.pi*self.sigma) - 0.5*((element_array - self.mu) / self.sigma)**2 - bad_samples = np.where(element_array <= 0)[0] - lnprob[bad_samples] = -np.inf + if no_negatives: - return lnprob + bad_samples = np.where(samples <= 0)[0] + lnprob[bad_samples] = -np.inf + + return lnprob class JeffreysPrior(Prior): """ From d9ba726928bc01770940859f03bdb699d792c2c7 Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 18:38:11 -0500 Subject: [PATCH 5/8] moved the `no_negatives` keyword to __init__, as per Jason's suggestion --- orbitize/priors.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/orbitize/priors.py b/orbitize/priors.py index 453a4ffe..a3b0f551 100644 --- a/orbitize/priors.py +++ b/orbitize/priors.py @@ -38,25 +38,26 @@ class GaussianPrior(Prior): Args: mu (float): mean of the distribution sigma (float): standard deviation of the distribution + no_negatives (bool): if True, only positive values will be drawn from + this prior, and the probability of negative values will be 0 (default:True). (written) Sarah Blunt, 2018 """ - def __init__(self, mu, sigma): + def __init__(self, mu, sigma, no_negatives=True): self.mu = mu self.sigma = sigma + self.no_negatives = no_negatives def __repr__(self): return "Gaussian" - def draw_samples(self, num_samples, no_negatives=True): + def draw_samples(self, num_samples): """ Draw positive samples from a Gaussian distribution. Negative samples will not be returned. Args: num_samples (float): the number of samples to generate - no_negatives (bool): if True, only positive samples will - be drawn (default: False). Returns: numpy array of float: samples drawn from the appropriate @@ -68,7 +69,7 @@ def draw_samples(self, num_samples, no_negatives=True): ) bad = np.inf - if no_negatives: + if self.no_negatives: while bad != 0: @@ -81,7 +82,7 @@ def draw_samples(self, num_samples, no_negatives=True): return samples - def compute_lnprob(self, element_array, no_negatives=True): + def compute_lnprob(self, element_array): """ Compute log(probability) of an array of numbers wrt a Gaussian distibution. Negative numbers return a probability of -inf. @@ -90,8 +91,6 @@ def compute_lnprob(self, element_array, no_negatives=True): element_array (float or np.array of float): array of numbers. We want the probability of drawing each of these from the appopriate Gaussian distribution - no_negatives (bool): if True, the log(probability) of negative numbers - will be -inf (default:True). Returns: numpy array of float: array of log(probability) values, @@ -100,7 +99,7 @@ def compute_lnprob(self, element_array, no_negatives=True): """ lnprob = -0.5*np.log(2.*np.pi*self.sigma) - 0.5*((element_array - self.mu) / self.sigma)**2 - if no_negatives: + if self.no_negatives: bad_samples = np.where(samples <= 0)[0] lnprob[bad_samples] = -np.inf From 965babd35cade0e7a1225b4925bbef76c4aa4151 Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 18:45:58 -0500 Subject: [PATCH 6/8] variable name bugfix --- orbitize/priors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/orbitize/priors.py b/orbitize/priors.py index a3b0f551..648fc9aa 100644 --- a/orbitize/priors.py +++ b/orbitize/priors.py @@ -101,7 +101,7 @@ def compute_lnprob(self, element_array): if self.no_negatives: - bad_samples = np.where(samples <= 0)[0] + bad_samples = np.where(element_array <= 0)[0] lnprob[bad_samples] = -np.inf return lnprob From 20c668893ff2f0e2899ee2c90a3ff7b4839b8b8f Mon Sep 17 00:00:00 2001 From: sblunt Date: Fri, 28 Dec 2018 18:53:24 -0500 Subject: [PATCH 7/8] Make numerical prior tests a little more precise --- tests/test_priors.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_priors.py b/tests/test_priors.py index 9807f8ab..3a0a13d6 100644 --- a/tests/test_priors.py +++ b/tests/test_priors.py @@ -7,7 +7,7 @@ threshold = 1e-1 initialization_inputs = { - priors.GaussianPrior : [100., 1.], + priors.GaussianPrior : [1000., 1.], priors.JeffreysPrior : [1., 2.], priors.UniformPrior : [0., 1.], priors.SinPrior : [], @@ -15,7 +15,7 @@ } expected_means_mins_maxes = { - priors.GaussianPrior : (100.,0.,np.inf), + priors.GaussianPrior : (1000.,0.,np.inf), priors.JeffreysPrior : (1/np.log(2),1., 2.), priors.UniformPrior : (0.5, 0., 1.), priors.SinPrior : (np.pi/2., 0., np.pi), @@ -23,7 +23,7 @@ } lnprob_inputs = { - priors.GaussianPrior : np.array([-3.0, np.inf, 100., 99.]), + priors.GaussianPrior : np.array([-3.0, np.inf, 1000., 999.]), priors.JeffreysPrior : np.array([-1., 0., 1., 1.5, 2., 2.5]), priors.UniformPrior : np.array([0., 0.5, 1., -1., 2.]), priors.SinPrior : np.array([0., np.pi/2., np.pi, 10., -1.]), @@ -31,7 +31,7 @@ } expected_probs = { - priors.GaussianPrior : np.array([0., 0., nm(100.,1.).pdf(100.), nm(100.,1.).pdf(99.)]), + priors.GaussianPrior : np.array([0., 0., nm(1000.,1.).pdf(1000.), nm(1000.,1.).pdf(999.)]), priors.JeffreysPrior : np.array([0., 0., 1., 2./3., 0.5, 0.])/np.log(2), priors.UniformPrior : np.array([1., 1., 1., 0., 0.]), priors.SinPrior : np.array([0., 0.5, 0., 0., 0.]), From 7ae867d71a88cc8748c7b977b75e11eaed9bfa23 Mon Sep 17 00:00:00 2001 From: sblunt Date: Sat, 29 Dec 2018 13:56:46 -0500 Subject: [PATCH 8/8] Update equality statement in checking for negative numbers <= --> <, as per @henry-ngo 's suggestion --- orbitize/priors.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbitize/priors.py b/orbitize/priors.py index 648fc9aa..6288ca90 100644 --- a/orbitize/priors.py +++ b/orbitize/priors.py @@ -73,7 +73,7 @@ def draw_samples(self, num_samples): while bad != 0: - bad_samples = np.where(samples <= 0)[0] + bad_samples = np.where(samples < 0)[0] bad = len(bad_samples) samples[bad_samples] = np.random.normal( @@ -101,7 +101,7 @@ def compute_lnprob(self, element_array): if self.no_negatives: - bad_samples = np.where(element_array <= 0)[0] + bad_samples = np.where(element_array < 0)[0] lnprob[bad_samples] = -np.inf return lnprob