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

Author: pberkes

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/_posterior.py

@@ -146,7 +146,6 @@ def log_posterior(data: np.ndarray, sigmoid: Sigmoid, priors: Dict[str, Prior],
 
     if gamma is None:
         gamma = lambd
-    scale = 1 - gamma - lambd
 
     pbin = 0
     p = 0
@@ -159,7 +158,7 @@ def log_posterior(data: np.ndarray, sigmoid: Sigmoid, priors: Dict[str, Prior],
         if trials == 0:
             continue
 
-        psi = sigmoid(level, thres, width) * scale + gamma
+        psi = sigmoid(level, thres, width, gamma=gamma, lambd=lambd)
 
         # Separate cases to avoid warnings problems with np.log(...)
         # for correct_trials == 0 or (trials - correct_trials) == 0

psignifit/psigniplot.py

@@ -56,8 +56,13 @@ def plot_psychometric_function(result: Result,  # noqa: C901, this function is t
     x = np.linspace(x_data.min(), x_data.max(), num=1000)
     x_low = np.linspace(x[0] - extrapolate_stimulus * (x[-1] - x[0]), x[0], num=100)
     x_high = np.linspace(x[-1], x[-1] + extrapolate_stimulus * (x[-1] - x[0]), num=100)
-    y = sigmoid(np.r_[x_low, x, x_high], params['threshold'], params['width'])
-    y = (1 - params['gamma'] - params['lambda']) * y + params['gamma']
+    y = sigmoid(
+        np.r_[x_low, x, x_high],
+        threshold=params['threshold'],
+        width=params['width'],
+        gamma=params['gamma'],
+        lambd=params['lambda'],
+    )
     ax.plot(x, y[len(x_low):-len(x_high)], c=line_color, lw=line_width, clip_on=False)
     ax.plot(x_low, y[:len(x_low)], '--', c=line_color, lw=line_width, clip_on=False)
     ax.plot(x_high, y[-len(x_high):], '--', c=line_color, lw=line_width, clip_on=False)
@@ -110,8 +115,13 @@ def _plot_residuals(x_values: np.ndarray,
     data = result.data
     sigmoid = result.configuration.make_sigmoid()
 
-    std_model = params['gamma'] + (1 - params['lambda'] - params['gamma']) * sigmoid(
-        data[:, 0], params['threshold'], params['width'])
+    std_model = sigmoid(
+        data[:, 0],
+        threshold=params['threshold'],
+        width=params['width'],
+        gamma=params['gamma'],
+        lambd=params['lambda'],
+    )
     deviance = data[:, 1] / data[:, 2] - std_model
     std_model = np.sqrt(std_model * (1 - std_model))
     deviance = deviance / std_model
@@ -329,8 +339,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],
@@ -348,8 +359,13 @@ def plot_prior(result: Result,
         for param_value, color in zip(x_percentiles, colors):
             this_sigmoid_params = dict(sigmoid_params)
             this_sigmoid_params[param] = param_value
-            y = sigmoid(sigmoid_x, this_sigmoid_params['threshold'], this_sigmoid_params['width'])
-            y = (1 - estimate['gamma'] - this_sigmoid_params['lambda']) * y + estimate['gamma']
+            y = sigmoid(
+                sigmoid_x,
+                threshold=this_sigmoid_params['threshold'],
+                width=this_sigmoid_params['width'],
+                gamma=estimate['gamma'],
+                lambd=this_sigmoid_params['lambda'],
+            )
             plt.plot(sigmoid_x, y, linewidth=line_width, color=color)
 
         plt.scatter(data[:, 0], np.zeros(data[:, 0].shape), s=marker_size*.75, c='k', clip_on=False)
@@ -376,13 +392,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))
@@ -428,7 +472,7 @@ def plot_bias_analysis(data: np.ndarray, compare_data: np.ndarray,
 
     fig = plt.figure(constrained_layout=True, figsize=(5, 15))
     gs = fig.add_gridspec(6, 1)
-    
+
     ax1 = fig.add_subplot(gs[0:2, 0])
     plot_psychometric_function(result_combined, ax=ax1, estimate_type=estimate_type)
     plot_psychometric_function(result_data, ax=ax1, line_color=[1, 0, 0], data_color=[1, 0, 0],
@@ -440,25 +484,25 @@ def plot_bias_analysis(data: np.ndarray, compare_data: np.ndarray,
     ax3 = fig.add_subplot(gs[3, 0])
     ax4 = fig.add_subplot(gs[4, 0])
     ax5 = fig.add_subplot(gs[5, 0])
-    
+
     axesmarginals = [ax2, ax3, ax4, ax5]
-    
+
     for param, ax in zip(['threshold', 'width', 'lambda', 'gamma'], axesmarginals):
 
-        plot_marginal(result_combined, param, ax=ax, plot_prior=False, 
+        plot_marginal(result_combined, param, ax=ax, plot_prior=False,
                       line_color=[0, 0, 0], estimate_type=estimate_type,
                       plot_ci=False)
-        
+
         plot_marginal(result_data, param, ax=ax, plot_prior=False,
-                      line_color=[1, 0, 0], estimate_type=estimate_type, 
+                      line_color=[1, 0, 0], estimate_type=estimate_type,
                       plot_ci=False)
-        
-        
+
+
         plot_marginal(result_compare_data, param, ax=ax, plot_prior=False,
                       line_color=[0, 0, 1], estimate_type=estimate_type,
                       plot_ci=False)
-     
+
     for ax in axesmarginals:
         ax.autoscale()
-     
+
 

psignifit/sigmoids.py

@@ -58,33 +58,38 @@ def __eq__(self, o: object) -> bool:
                 and o.alpha == self.alpha
                 and o.negative == self.negative)
 
-    def __call__(self, stimulus_level: N, threshold: N, width: N) -> N:
+    def __call__(self, stimulus_level: N, threshold: N, width: N, gamma: N = 0, lambd: N = 0) -> N:
         """ Calculate the sigmoid value at specified stimulus levels.
 
+        See Eq 1 in Schuett, Harmeling, Macke and Wichmann (2016).
+
         Args:
-            stimulus_level: Stimulus level value at which to calculate the slope
-            threshold: Parameter value for threshold at PC
-            width: Parameter value for width of the sigmoid
+            stimulus_level: Stimulus level value
+            threshold: Threshold at `PC`
+            width: Width of the sigmoid
+            gamma: Guess rate (lower asymptote of the sigmoid)
+            lambd: Lapse rate (upper asymptote of the sigmoid)
         Returns:
             Proportion correct at the stimulus values.
         """
 
         value = self._value(stimulus_level, threshold, width)
+        value = gamma + (1.0 - lambd - gamma) * value
 
         if self.negative:
-            return 1 - value
-        else:
-            return value
+            value = 1 - value
+
+        return value
 
     def slope(self, stimulus_level: N, threshold: N, width: N, gamma: N = 0, lambd: N = 0) -> N:
         """ Calculate the slope at specified stimulus levels.
 
         Args:
             stimulus_level: Stimulus level value at which to calculate the slope
-            threshold: Parameter value for threshold at PC
-            width: Parameter value for width of the sigmoid
-            gamma: Parameter value for the lower offset of the sigmoid
-            lambd: Parameter value for the upper offset of the sigmoid
+            threshold: Threshold at `PC`
+            width: Width of the sigmoid
+            gamma: Guess rate (lower asymptote of the sigmoid)
+            lambd: Lapse rate (upper asymptote of the sigmoid)
         Returns:
             Slope at the stimulus values.
         """
@@ -105,10 +110,10 @@ def inverse(self, prop_correct: N, threshold: N, width: N,
 
         Args:
             prop_correct: Proportion correct at the threshold to calculate.
-            threshold: Parameter value for threshold at PC
-            width: Parameter value for width of the sigmoid
-            gamma: Parameter value for the lower offset of the sigmoid
-            lambd: Parameter value for the upper offset of the sigmoid
+            threshold: Threshold at `PC`
+            width: Width of the sigmoid
+            gamma: Guess rate (lower asymptote of the sigmoid)
+            lambd: Lapse rate (upper asymptote of the sigmoid)
         Returns:
             Threshold at the proportion correct values.
         """
@@ -151,7 +156,7 @@ def _value(self, stimulus_level: N, threshold: N, width: N) -> N:
 
         Args:
             stimulus_level: Stimulus level values at which to calculate the sigmoid value
-            threshold: Threshold value at PC
+            threshold: Threshold at `PC`
             width: Width of the sigmoid
         Returns:
             Proportion correct at the stimulus values.
@@ -168,7 +173,7 @@ def _slope(self, stimulus_level: N, threshold: N, width: N) -> N:
 
         Args:
             stimulus_level: Stimulus level value at which to calculate the slope
-            threshold: Threshold value at PC
+            threshold: Threshold at `PC`
             width: Width of the sigmoid
         Returns:
             Slope at the stimulus level value
@@ -185,7 +190,7 @@ def _inverse(self, prop_correct: N, threshold: N, width: N) -> N:
 
         Args:
             prop_correct: Proportion correct values at which to calculate the stimulus level values.
-            threshold: Threshold value at PC
+            threshold: Threshold at `PC`
             width: Width of the sigmoid
         Returns:
             Stimulus values corresponding to the proportion correct values.
@@ -206,7 +211,7 @@ def _standard_parameters(self, threshold: N, width: N) -> list:
         For negative slope sigmoids, we return the same parameters as for the positive ones.
 
         Args:
-            threshold: Threshold value at PC
+            threshold: Threshold at `PC`
             width: Width of the sigmoid
         Returns:
             Standard parameters (loc, scale) for the sigmoid subclass.
@@ -376,8 +381,8 @@ def assert_sigmoid_sanity_checks(sigmoid, n_samples: int, threshold: float, widt
 
     Args:
          n_samples: Number of stimulus levels between 0 (exclusive) and 1 for tests
-         threshold: Parameter value for threshold at PC
-         width: Width of the sigmoid
+            threshold: Threshold at `PC`
+            width: Width of the sigmoid
     Raises:
           AssertionError if a sanity check fails.
     """