Merge remote-tracking branch 'origin/main' into move-tests-and-demos-out

psignifit/_confidence.py

@@ -14,24 +14,26 @@
 
 
 def confidence_intervals(probability_mass: np.ndarray, grid_values: np.ndarray,
-                         p_values: Sequence[float], mode: str) -> Dict[str, list]:
+                         confidence_levels: Sequence[float], mode: str) -> Dict[str, list]:
     """ Confidence intervals on probability grid.
 
     Supports two methods:
 
         - 'project', projects the confidence region down each axis.
           Implemented using :func:`grid_hdi`.
-        - 'percentiles', finds alpha/2 and 1-alpha/2 percentiles (alpha = 1-p_value).
+        - 'percentiles', finds alpha/2 and 1-alpha/2 percentiles (alpha = 1 - confidence_level).
           Implemented using :func:`percentile_intervals`.
 
     Args:
         probability_mass: Probability mass at each grid point, shape (n_points, n_points, ...)
         grid_values: Parameter values along grid axis in the same order as zerocentered_normal_mass dimensions,
                      shape (n_dims, n_points)
-        p_values: Probabilities of confidence in the intervals.
+        confidence_levels: Fraction of the posterior to be included in the corresponding confidence interval. For
+            example, a confidence level of 0.92 will return a confidence interval containing 92% of the probability
+            mass in the posterior
         mode: Either 'project' or 'percentiles'.
     Returns:
-        A dictionary mapping p_values as a string to a list containing the start and end grid-values for the
+        A dictionary mapping confidence_levels as a string to a list containing the start and end grid-values for the
         confidence interval per dimension, shape (n_dims, 2).
     Raises:
         ValueError for unsupported mode or sum(probability_mass) != 1.
@@ -45,8 +47,8 @@ def confidence_intervals(probability_mass: np.ndarray, grid_values: np.ndarray,
     if not np.isclose(probability_mass.sum(), 1):
         raise ValueError(f'Expects sum(probability_mass) to be 1., got {probability_mass.sum():.4f}')
     intervals = {}
-    for p_value in p_values:
-        intervals[str(p_value)] = calc_ci(probability_mass, grid_values, p_value).tolist()
+    for level in confidence_levels:
+        intervals[str(level)] = calc_ci(probability_mass, grid_values, level).tolist()
 
     return intervals
 
@@ -67,7 +69,7 @@ def grid_hdi(probability_mass: np.ndarray, grid_values: np.ndarray, credible_mas
         grid_values: Parameter values along grid axis, shape (n_dims, n_points)
         credible_mass: Minimal mass in highest density region
     Returns:
-        Grid value at interval per dimension, shape (n_dims, 2)
+        Start and end grid-values for the confidence interval per dimension, shape (n_dims, 2)
 
     .. _stats.stackexchange: https://stats.stackexchange.com/questions/148439/what-is-a-highest-density-region-hdr
     """
@@ -87,21 +89,21 @@ def grid_hdi(probability_mass: np.ndarray, grid_values: np.ndarray, credible_mas
     return intervals
 
 
-def percentile_intervals(probability_mass: np.ndarray, grid_values: np.ndarray, p_value: float):
+def percentile_intervals(probability_mass: np.ndarray, grid_values: np.ndarray, confidence_level: float):
     """ Percentile intervals on a grid.
 
-    Find alpha/2 and 1-alpha/2 percentiles on marginal probability mass (alpha = 1 - p_value).
+    Find alpha/2 and 1-alpha/2 percentiles on marginal probability mass (alpha = 1 - confidence_level).
 
     Args:
         probability_mass: Probability mass at each grid point, shape (n_points, n_points, ...)
         grid_values: Parameter values along grid axis, shape (n_dims, n_points)
-        p_value: Probability mass within the returned confidence bounds.
+        confidence_level: Probability mass within the returned confidence bounds.
     Returns:
-        Grid value at interval per dimension, shape (n_dims, 2)
+        Start and end grid-values for the confidence interval per dimension, shape (n_dims, 2)
     """
     mass_margins = stats.contingency.margins(probability_mass)
     intervals = np.empty((len(mass_margins), 2))
-    alpha = 1 - p_value
+    alpha = 1 - confidence_level
     for d, mass in enumerate(mass_margins):
         intervals[d, :] = np.interp([alpha / 2, 1 - alpha / 2], mass.cumsum(), grid_values[d])
     return intervals

psignifit/psigniplot.py

@@ -329,8 +329,9 @@ def plot_prior(result: Result,
         prior_cdf = np.cumsum(prior_val * prior_w)
         q25_index = np.argmax(prior_cdf > 0.25)
         q75_index = np.argmax(prior_cdf > 0.75)
+        prior_mean = np.sum(prior_x * prior_val)/np.sum(prior_val)
 
-        x_percentiles = [estimate[param],
+        x_percentiles = [prior_mean,
                          min(prior_x),
                          prior_x[q25_index],
                          prior_x[q75_index],
@@ -376,13 +377,41 @@ def plot_2D_margin(result: Result,
     other_param_ix = tuple(i for param, i in parameter_indices.items()
                            if param != first_param and param != second_param)
     marginal_2d = np.sum(result.debug['posteriors'], axis=other_param_ix)
+    extent = [result.parameter_values[second_param][0], result.parameter_values[second_param][-1],
+              result.parameter_values[first_param][-1], result.parameter_values[first_param][0]]
+
     if len(np.squeeze(marginal_2d).shape) != 2 or np.any(np.array(marginal_2d.shape) == 1):
-        raise ValueError('The marginal is not two-dimensional. Were the parameters fixed during optimization? If so, then change the arguments to parametes that were unfixed, or use plot_marginal() to obtain a 1D marginal for a parameter.')
+        len_first = len(result.parameter_values[first_param])
+        len_second = len(result.parameter_values[second_param])
+
+        # if first_param is singleton, we copy the marginal into a matrix
+        if len_first == 1 and len_second != 1:
+            marginal_2d = np.broadcast_to(marginal_2d,
+                                          (len(result.parameter_values[second_param]),
+                                           len(result.parameter_values[second_param]))
+                                          )
+            extent[2] = 1  # replace range for a mockup range between 0 and 1
+            extent[3] = 0
+
+        # if second param is singleton
+        elif len_first != 1 and len_second == 1:
+            marginal_2d = np.broadcast_to(marginal_2d,
+                                          (len(result.parameter_values[first_param]),
+                                           len(result.parameter_values[first_param]))
+                                          )
+            extent[0] = 0
+            extent[1] = 1
+
+        # if both params are singletons, we return a matrix full of ones
+        elif len_first == 1 and len_second == 1:
+            marginal_2d = np.ones((len(result.parameter_values[first_param]),
+                                   len(result.parameter_values[second_param]))
+                                  )
+            extent = [0, 1, 1, 0]
 
     if parameter_indices[first_param] > parameter_indices[second_param]:
         marginal_2d = np.transpose(marginal_2d)
-    extent = [result.parameter_values[second_param][0], result.parameter_values[second_param][-1],
-              result.parameter_values[first_param][-1], result.parameter_values[first_param][0]]
+
     ax.imshow(marginal_2d, extent=extent, cmap='Reds_r',  aspect='auto')
     ax.set_xlabel(_parameter_label(second_param))
     ax.set_ylabel(_parameter_label(first_param))

tests/test_confidence.py

@@ -44,20 +44,20 @@ def grid_values():
 
 
 def test_confidence_intervals(zerocentered_normal_mass, grid_values):
-    p_values = [0.05, 0.5, 0.95]
+    confidence_levels = [0.05, 0.5, 0.95]
 
-    intervals = confidence_intervals(zerocentered_normal_mass, grid_values, p_values, mode='project')
-    assert list(intervals.keys()) == [str(p) for p in p_values]
+    intervals = confidence_intervals(zerocentered_normal_mass, grid_values, confidence_levels, mode='project')
+    assert list(intervals.keys()) == [str(p) for p in confidence_levels]
     assert all(len(interval) == len(grid_values) for interval in intervals.values())
 
-    intervals = confidence_intervals(zerocentered_normal_mass, grid_values, p_values, mode='percentiles')
-    assert list(intervals.keys()) == [str(p) for p in p_values]
+    intervals = confidence_intervals(zerocentered_normal_mass, grid_values, confidence_levels, mode='percentiles')
+    assert list(intervals.keys()) == [str(p) for p in confidence_levels]
     assert all(len(interval) == len(grid_values) for interval in intervals.values())
 
     with pytest.raises(ValueError):
-        confidence_intervals(zerocentered_normal_mass, grid_values, p_values, mode='foobar')
+        confidence_intervals(zerocentered_normal_mass, grid_values, confidence_levels, mode='foobar')
     with pytest.raises(ValueError):
-        confidence_intervals(zerocentered_normal_mass * 3, grid_values, p_values, mode='project')
+        confidence_intervals(zerocentered_normal_mass * 3, grid_values, confidence_levels, mode='project')
 
 
 def test_grid_hdi(zerocentered_normal_mass, grid_values):