Source code for threeML.random_variates

import numpy as np

from import uncertainty_formatter

[docs]class RandomVariates(np.ndarray): """ A subclass of np.array which is meant to contain samples for one parameter. This class contains methods to easily compute properties for the parameter (errors and so on) """ def __new__(cls, input_array, value=None): # Input array is an already formed ndarray instance # We first cast to be our class type obj = np.asarray(input_array).view(cls) # Add the value obj._orig_value = value # Finally, we must return the newly created object: return obj def __array_finalize__(self, obj): # see InfoArray.__array_finalize__ for comments if obj is None: return # Add the value self._orig_value = getattr(obj, "_orig_value", None) def __array_wrap__(self, out_arr, context=None): # This gets called at the end of any operation, where out_arr is the result of the operation # We need to update _orig_value so that the final results will have it out_arr._orig_value = out_arr.median # then just call the parent return super(RandomVariates, self).__array_wrap__(out_arr, context) # def __array_ufunc__(self, ufunc, method, *inputs, **kwargs): # # TODO: must make this return single numbers is needed # args = [] # in_no = [] # for i, input_ in enumerate(inputs): # if isinstance(input_, RandomVariates): # in_no.append(i) # args.append(input_.view(np.ndarray)) # else: # args.append(input_) # outputs = kwargs.pop('out', None) # out_no = [] # if outputs: # out_args = [] # for j, output in enumerate(outputs): # if isinstance(output, RandomVariates): # out_no.append(j) # out_args.append(output.view(np.ndarray)) # else: # out_args.append(output) # kwargs['out'] = tuple(out_args) # else: # outputs = (None,) * ufunc.nout # results = super(RandomVariates, self).__array_ufunc__(ufunc, method, # *args, **kwargs) # if results is NotImplemented: # return NotImplemented # if method == 'at': # return # if ufunc.nout == 1: # results = (results,) # results = tuple((np.asarray(result).view(RandomVariates) # if output is None else output) # for result, output in zip(results, outputs)) # return results[0] if len(results) == 1 else results @property def median(self): """Returns median value""" # the np.asarray casting avoids the calls to __new__ and __array_finalize_ of this class return float(np.median(np.asarray(self))) # @property # def mean(self): # """Returns average value""" # return float(np.asarray(self).mean()) @property def std(self): """Returns sample std value""" return float(np.asarray(self).std()) @property def var(self): """Returns sample variance value""" return float(np.asarray(self).var()) @property def average(self): """Returns average value""" return float(np.asarray(self).mean()) @property def value(self): return float(self._orig_value) @property def samples(self): return np.asarray(self)
[docs] def highest_posterior_density_interval(self, cl=0.68): """ Returns the Highest Posterior Density interval (HPD) for the parameter, for the given credibility level. NOTE: the returned interval is the HPD only if the posterior is not multimodal. If it is multimodal, you should probably report the full posterior, not only an interval. :param cl: credibility level (0 < cl < 1) :return: (low_bound, hi_bound) """ assert 0 < cl < 1, "The credibility level should be 0 < cl < 1" # NOTE: we cannot sort the array, because we would destroy the covariance with other physical quantities, # so we get a copy instead. This copy will live only for the duration of this method (but of course will be # collected only whenevery the garbage collector decides to). ordered = np.sort(np.array(self)) n = ordered.size # This is the probability that the interval should span interval_integral = cl # If all values have the same probability, then the hpd is degenerate, but its length is from 0 to # the value corresponding to the (interval_integral * n)-th sample. # This is the index of the rightermost element which can be part of the interval index_of_rightmost_possibility = int(np.floor(interval_integral * n)) # Compute the index of the last element that is eligible to be the left bound of the interval index_of_leftmost_possibility = n - index_of_rightmost_possibility # Now compute the width of all intervals that might be the one we are looking for interval_width = ( ordered[index_of_rightmost_possibility:] - ordered[:index_of_leftmost_possibility] ) # This might happen if there are too few values if len(interval_width) == 0: raise RuntimeError("Too few elements for interval calculation") # Find the index of the shortest interval idx_of_minimum = np.argmin(interval_width) # Find the extremes of the shortest interval hpd_left_bound = ordered[idx_of_minimum] hpd_right_bound = ordered[idx_of_minimum + index_of_rightmost_possibility] return hpd_left_bound, hpd_right_bound
[docs] def equal_tail_interval(self, cl=0.68): """ Returns the equal tail interval, i.e., an interval centered on the median of the distribution with the same probability on the right and on the left of the mean. If the distribution of the parameter is Gaussian and cl=0.68, this is equivalent to the 1 sigma confidence interval. :param cl: confidence level (0 < cl < 1) :return: (low_bound, hi_bound) """ assert 0 < cl < 1, "Confidence level must be 0 < cl < 1" half_cl = cl / 2.0 * 100.0 low_bound, hi_bound = np.percentile( np.asarray(self), [50.0 - half_cl, 50.0 + half_cl] ) return float(low_bound), float(hi_bound)
# np.ndarray already has a mean() and a std() methods def __repr__(self): # Get representation for the HPD min_bound, max_bound = self.highest_posterior_density_interval(0.68) hpd_string = uncertainty_formatter(self.median, min_bound, max_bound) # Get representation for the equal-tail interval min_bound, max_bound = self.equal_tail_interval(0.68) eqt_string = uncertainty_formatter(self.median, min_bound, max_bound) # Put them together representation = "equal-tail: %s, hpd: %s" % (eqt_string, hpd_string) return representation def __str__(self): return self.__repr__()