substantial changes to the priors set up, and mcmc and added strong s…

…upport for the unaccounted uncertainty

uravu/optimize.py

@@ -78,15 +78,14 @@ def ln_likelihood(
         (float): ln-likelihood between model and data.
     """
     if unaccounted_uncertainty:
-        variables = variables[:-1]
-        uu_f = variables[-1]
+        var = variables[:-1]
+        log_f = variables[-1]
     else:
-        uu_f = 0
-    model = function(abscissa.m, *variables)
+        var = variables
+        log_f = -np.inf
     y_data = unp.nominal_values(ordinate.m)
     dy_data = unp.std_devs(ordinate.m)
+    model = function(abscissa.m, *var)
 
-    sigma2 = dy_data ** 2 + uu_f ** 2 * model ** 2
-    return -0.5 * np.sum(
-        (model - y_data) ** 2 / sigma2 + np.log(2 * np.pi * sigma2)
-    )
+    sigma2 = dy_data ** 2 + model ** 2 * np.exp(2 * log_f)
+    return -0.5 * np.sum((model - y_data) ** 2 / sigma2 + np.log(sigma2))

uravu/plotting.py

@@ -69,8 +69,8 @@ def plot_relationship(
             )
     if relationship.unaccounted_uncertainty:
         var = relationship.variable_medians
-        additional_uncertainty = var[-1] * relationship.function(
-            relationship.x_n, *var[:-1]
+        additional_uncertainty = np.abs(
+            var[-1] * relationship.function(relationship.x_n, *var[:-1])
         )
         axes.fill_between(
             relationship.x_n,

uravu/relationship.py

@@ -478,29 +478,25 @@ def max_likelihood(self, x0=None, **kwargs):
         """
         self.variables = optimize.max_ln_likelihood(self, x0, **kwargs)
 
-    def prior(self, array):
+    def prior(self):
         """
-        *Standard priors*. Broad, uniform distributions spread evenly
-        around the current values for the variables.
-
-        This is used as the default where no priors are given.
-
-        Args:
-            array (array_like): An array of random uniform numbers (0, 1].
-                The shape of which is M x N, where M is the number of
-                parameters and N is the number of walkers.
+        *Standard priors* for the relationship. These priors are broad,
+        uninformative, for normal variables running the range
+        [x - x * 10, x + x * 10) (where x is the variable). For an unaccounted
+        uncertainty natural log factor the range is [-10, 1).
 
         Returns:
-            (array_like): an array of random uniform numbers broadly
-            distributed in the range [x - x * 5, x + x * 5), where x is the
-            current variable value.
-        """
-        broad = np.copy(array)
-        for i, variable in enumerate(self.variable_medians):
-            loc = variable - variable * 5
-            scale = (variable + variable * 5) - loc
-            broad[i] = uniform.ppf(broad[i], loc=loc, scale=scale)
-        return broad
+            (list of scipy.stats.rv_continuous): scipy.stats functions
+                describing the priors.
+        """
+        priors = []
+        for var in self.variable_medians:
+            loc = var - np.abs(var) * 10
+            scale = (var + np.abs(var) * 10) - loc
+            priors.append(uniform(loc=loc, scale=scale))
+        if self.unaccounted_uncertainty:
+            priors[-1] = uniform(loc=-10, scale=11)
+        return priors
 
     def mcmc(
         self,

uravu/sampling.py

@@ -49,22 +49,25 @@ def mcmc(
         (list of uravu.distribution.Distribution): a
             uravu.distribution.Distribution to describe each of the variables.
     """
-    uniform = np.random.uniform(size=(len(relationship.variables), walkers))
     if prior_function is None:
-        initial_prior = relationship.prior(uniform).T
-    else:
-        initial_prior = prior_function(uniform).T
+        prior_function = relationship.prior
+
+    initial_prior = relationship.variable_medians + 1e-4 * np.random.randn(
+        walkers, len(relationship.variable_medians)
+    )
     ndims = initial_prior.shape[1]
+    k = prior_function()
 
     sampler = emcee.EnsembleSampler(
         walkers,
         ndims,
-        optimize.ln_likelihood,
+        ln_probability,
         args=[
             relationship.function,
             relationship.abscissa,
             relationship.ordinate,
             relationship.unaccounted_uncertainty,
+            k,
         ],
     )
 
@@ -79,6 +82,20 @@ def mcmc(
     return distributions
 
 
+def ln_probability(
+    variables, function, abscissa, ordinate, unaccounted_uncertainty, priors
+):
+    """
+    Determine the natural log probability for a given set of variables.
+    """
+    log_prior = 0
+    for i, var in enumerate(variables):
+        log_prior += priors[i].logpdf(var)
+    return log_prior + optimize.ln_likelihood(
+        variables, function, abscissa, ordinate, unaccounted_uncertainty
+    )
+
+
 def nested_sampling(
     relationship, prior_function=None, progress=True, **kwargs
 ):
@@ -102,17 +119,42 @@ def nested_sampling(
     """
     if prior_function is None:
         prior_function = relationship.prior
+    priors = prior_function()
     sampler = dynesty.NestedSampler(
         optimize.ln_likelihood,
-        prior_function,
+        nested_prior,
         len(relationship.variables),
         logl_args=[
             relationship.function,
             relationship.abscissa,
             relationship.ordinate,
             relationship.unaccounted_uncertainty,
         ],
+        ptform_args=[priors,],
     )
     sampler.run_nested(print_progress=progress, **kwargs)
     results = sampler.results
     return ufloat(results["logz"][-1], results["logzerr"][-1])
+
+
+def nested_prior(array, priors):
+    """
+    *Standard priors*. Broad, uniform distributions spread evenly
+    around the current values for the variables.
+
+    This is used as the default where no priors are given.
+
+    Args:
+        array (array_like): An array of random uniform numbers (0, 1].
+            The shape of which is M x N, where M is the number of
+            parameters and N is the number of walkers.
+
+    Returns:
+        (array_like): an array of random uniform numbers broadly
+        distributed in the range [x - x * 5, x + x * 5), where x is the
+        current variable value.
+    """
+    broad = np.copy(array)
+    for i, prior in enumerate(priors):
+        broad[i] = prior.ppf(broad[i])
+    return broad