Merge pull request #184 from pberkes/lambda_gamma_in_sigmoid

Move rescaling of psychometric function related to lambda and gamma in `Sigmoid` object

Author: pberkes

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
@@ -349,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)
@@ -457,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],
@@ -469,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.
     """

psignifit/tests/test_sigmoids.py

@@ -45,7 +45,9 @@ def test_sigmoid_by_name(sigmoid_name):
 def test_sigmoid_values(subclass, expected_y):
     sigmoid = subclass(PC=0.6, alpha=0.1)
     x = np.array([9.5, 10.0, 11.5])
-    y = sigmoid(x, threshold=10, width=3)
+    y = sigmoid(x, threshold=10, width=3, gamma=0.08, lambd=0.12)
+    # Rescale expected_y to take into account gamma and lambda
+    expected_y = 0.08 + expected_y * 0.8
     np.testing.assert_allclose(y, expected_y, atol=1e-6)
 
 

psignifit/tools.py

@@ -54,9 +54,9 @@ def pool_blocks(data: np.ndarray, max_tol=0, max_gap=np.inf, max_length=np.inf):
 def psychometric(stimulus_level, threshold, width, gamma, lambda_, sigmoid_name):
     """ Psychometric function aka proportion correct function.
 
-    Generates proportion correct values for a range of stimulus levels given a
-    sigmoid.
-    Implementation of Eq 1 in Schuett, Harmeling, Macke and Wichmann (2016)
+    This is a convenience function used mostly for testing and demos, to make the intent of creating a psychometric
+    function more explicit. The implementation simply combines creating a `Sigmoid` object using the sigmoid name,
+    and calling the object to generate the proportion correct values corresponding to `stimulus_level`.
 
     Parameters:
         stimulus_level: array
@@ -69,25 +69,59 @@ def psychometric(stimulus_level, threshold, width, gamma, lambda_, sigmoid_name)
             Guess rate
         lambda_: float
             Lapse rate
-        sigmoid: callable
-            Sigmoid function to use. Default is Gaussian
+        sigmoid_name: callable
+            Name of sigmoid function to use. See `psignifit.sigmoids.ALL_SIGMOID_NAMES` for the list of available
+            sigmoids.
 
     Returns:
         psi: array
-            proportion correct values for each stimulus level
+            Proportion correct values for each stimulus level
 
     """
-    # we use the defaults for pc and alpha in the sigmoids:
-    # pc = 0.5
-    # alpha = 0.05
     sigmoid = sigmoid_by_name(sigmoid_name)
-    sigmoid_values = sigmoid(stimulus_level, threshold=threshold, width=width)
-    psi = gamma + (1.0 - lambda_ - gamma) * sigmoid_values
+    psi = sigmoid(stimulus_level, threshold=threshold, width=width, gamma=gamma, lambd=lambda_)
     return psi
 
 
 def psychometric_with_eta(stimulus_level, threshold, width, gamma, lambda_,
-                 sigmoid_name, eta, random_state=np.random.RandomState(42)):
+                          sigmoid_name, eta, random_state=None):
+    """ Psychometric function with overdispersion.
+
+    This is a convenience function used mostly for testing and demos. Just like the function `psychometric`, it
+    computes proportion correct values for a given psychometric function type, specified by name. In addition,
+    it adds some additional noise to the data, so that its variance is compatible with the overdispersion
+    parameter `eta`.
+
+    See Section 2.2 in Schuett, Harmeling, Macke and Wichmann (2016).
+
+    Parameters:
+        stimulus_level: array
+          Values of the stimulus value
+        threshold: float
+            Threshold of the psychometric function
+        width: float
+            Width of the psychometric function
+        gamma: float
+            Guess rate
+        lambda_: float
+            Lapse rate
+        sigmoid_name: callable
+            Name of sigmoid function to use. See `psignifit.sigmoids.ALL_SIGMOID_NAMES` for the list of available
+            sigmoids.
+        eta: float
+            Overdispersion parameter
+        random_state: np.RandomState
+            Random state used to generate the additional variance in the data.
+            If None, NumPy's default random number generator is used.
+
+    Returns:
+        psi: array
+            Proportion correct values for each stimulus level
+
+    """
+
+    if random_state is None:
+        random_state = np.random.default_rng()
 
     psi = psychometric(stimulus_level, threshold, width, gamma, lambda_, sigmoid_name)
     new_psi = []