Coverage for bilby/core/utils/samples.py: 59%

1import numpy as np 

2from scipy.special import logsumexp 

3 

4 

5class SamplesSummary(object): 

6 """ Object to store a set of samples and calculate summary statistics 

7 

8 Parameters 

9 ========== 

10 samples: array_like 

11 Array of samples 

12 average: str {'median', 'mean'} 

13 Use either a median average or mean average when calculating relative 

14 uncertainties 

15 level: float 

16 The default confidence interval level, defaults t0 0.9 

17 

18 """ 

19 def __init__(self, samples, average='median', confidence_level=.9): 

20 self.samples = samples 

21 self.average = average 

22 self.confidence_level = confidence_level 

23 

24 @property 

25 def samples(self): 

26 return self._samples 

27 

28 @samples.setter 

29 def samples(self, samples): 

30 self._samples = samples 

31 

32 @property 

33 def confidence_level(self): 

34 return self._confidence_level 

35 

36 @confidence_level.setter 

37 def confidence_level(self, confidence_level): 

38 if 0 < confidence_level and confidence_level < 1: 

39 self._confidence_level = confidence_level 

40 else: 

41 raise ValueError("Confidence level must be between 0 and 1") 

42 

43 @property 

44 def average(self): 

45 if self._average == 'mean': 

46 return self.mean 

47 elif self._average == 'median': 

48 return self.median 

49 

50 @average.setter 

51 def average(self, average): 

52 allowed_averages = ['mean', 'median'] 

53 if average in allowed_averages: 

54 self._average = average 

55 else: 

56 raise ValueError("Average {} not in allowed averages".format(average)) 

57 

58 @property 

59 def median(self): 

60 return np.median(self.samples, axis=0) 

61 

62 @property 

63 def mean(self): 

64 return np.mean(self.samples, axis=0) 

65 

66 @property 

67 def _lower_level(self): 

68 """ The credible interval lower quantile value """ 

69 return (1 - self.confidence_level) / 2. 

70 

71 @property 

72 def _upper_level(self): 

73 """ The credible interval upper quantile value """ 

74 return (1 + self.confidence_level) / 2. 

75 

76 @property 

77 def lower_absolute_credible_interval(self): 

78 """ Absolute lower value of the credible interval """ 

79 return np.quantile(self.samples, self._lower_level, axis=0) 

80 

81 @property 

82 def upper_absolute_credible_interval(self): 

83 """ Absolute upper value of the credible interval """ 

84 return np.quantile(self.samples, self._upper_level, axis=0) 

85 

86 @property 

87 def lower_relative_credible_interval(self): 

88 """ Relative (to average) lower value of the credible interval """ 

89 return self.lower_absolute_credible_interval - self.average 

90 

91 @property 

92 def upper_relative_credible_interval(self): 

93 """ Relative (to average) upper value of the credible interval """ 

94 return self.upper_absolute_credible_interval - self.average 

95 

96 

97def kish_log_effective_sample_size(ln_weights): 

98 """ Calculate the Kish effective sample size from the natural-log weights 

99 

100 See https://en.wikipedia.org/wiki/Effective_sample_size for details 

101 

102 Parameters 

103 ========== 

104 ln_weights: array 

105 An array of the ln-weights 

106 

107 Returns 

108 ======= 

109 ln_n_eff: 

110 The natural-log of the effective sample size 

111 

112 """ 

113 return 2 * logsumexp(ln_weights) - logsumexp(2 * ln_weights) 

114 

115 

116def reflect(u): 

117 """ 

118 Iteratively reflect a number until it is contained in [0, 1]. 

119 

120 This is for priors with a reflective boundary condition, all numbers in the 

121 set `u = 2n +/- x` should be mapped to x. 

122 

123 For the `+` case we just take `u % 1`. 

124 For the `-` case we take `1 - (u % 1)`. 

125 

126 E.g., -0.9, 1.1, and 2.9 should all map to 0.9. 

127 

128 Parameters 

129 ========== 

130 u: array-like 

131 The array of points to map to the unit cube 

132 

133 Returns 

134 ======= 

135 u: array-like 

136 The input array, modified in place. 

137 """ 

138 idxs_even = np.mod(u, 2) < 1 

139 u[idxs_even] = np.mod(u[idxs_even], 1) 

140 u[~idxs_even] = 1 - np.mod(u[~idxs_even], 1) 

141 return u