diff --git a/psignifit/_result.py b/psignifit/_result.py index 85d04f8..29a8426 100644 --- a/psignifit/_result.py +++ b/psignifit/_result.py @@ -170,3 +170,25 @@ def slope_at_proportion_correct(self, proportion_correct: np.ndarray, unscaled: """ stimulus_levels = self.threshold(proportion_correct, unscaled, return_ci=False, estimate_type=estimate_type) return self.slope(stimulus_levels) + + def standard_parameters_estimate(self, estimate_type: Optional[EstimateType]=None): + """ Get the parameters of the psychometric function in the standard format. + + `psignifit` uses the same intuitive parametrization, threshold and width, for all + sigmoid types. However, each different type of sigmoid has its own standard parametrization. + + The interpretation of the standard parameters, location and scale, depends on the sigmoid class used. + For instance, for a Gaussian sigmoid, the location corresponds to the mean and the scale to the standard + deviation of the distribution. + + For negative slope sigmoids, we return the same parameters as for the positive ones. + + Args: + proportion_correct: proportion correct at the threshold you want to calculate + Returns: + Standard parameters (loc, scale) for the sigmoid subclass. + """ + sigmoid = self.configuration.make_sigmoid() + estimate = self.get_parameters_estimate(estimate_type) + loc, scale = sigmoid.standard_parameters(estimate['threshold'], estimate['width']) + return loc, scale diff --git a/psignifit/sigmoids.py b/psignifit/sigmoids.py index ec67391..8b7d734 100644 --- a/psignifit/sigmoids.py +++ b/psignifit/sigmoids.py @@ -5,24 +5,8 @@ from typing import Optional, TypeVar import numpy as np -import scipy as sp import scipy.stats -# Alias common distribution to be reused all over the place. - -# - Normal distribution: -# - This one is useful when we want mean=0, std=1 -# Percent point function -> inverse of cumulative normal distribution function -# returns percentiles -norminv = scipy.stats.norm(loc=0, scale=1).ppf -# - also instantiate a generic version -norminvg = scipy.stats.norm.ppf -# - Cumulative normal distribution function -normcdf = scipy.stats.norm.cdf - -# T-Student with df=1 -t1cdf = scipy.stats.t(1).cdf -t1icdf = scipy.stats.t(1).ppf # sigmoid can be calculated on single floats, or on numpy arrays of floats N = TypeVar('N', float, np.ndarray) @@ -31,11 +15,16 @@ class Sigmoid: """ Base class for sigmoid implementation. - Handles logarithmic input and negative output - for the specific sigmoid implementations. + Handles negative output for the specific sigmoid implementations. - Sigmoid classes should derive from this class and implement - the methods '_value', '_slope', and '_threshold'. + Sigmoid classes should derive from this class and implement the methods '_value', '_slope', '_threshold', + and `_standard_parameters`. + + In many cases, a sigmoid corresponds to the CDF of a probability distribution. This is the case + for all the sigmoids included in `psignifit`. If the probability distribution + is one of the distributions defined in `scipy.stats`, it is possible to define a `Sigmoid` subclass + by implementing only the methods `_scipy_distr` and `_standard_parameters`. The methods + '_value', '_slope', '_threshold' will be automatically defined in such a case. The stimulus levels, threshold and width are parameters of method calls. They correspond to the object attributes PC and alpha in the following way: @@ -73,14 +62,13 @@ def __call__(self, stimulus_level: N, threshold: N, width: N) -> N: """ Calculate the sigmoid value at specified stimulus levels. Args: - prop_correct: Proportion correct at the threshold to calculate. + 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 Returns: Proportion correct at the stimulus values. """ + value = self._value(stimulus_level, threshold, width) if self.negative: @@ -92,7 +80,7 @@ def slope(self, stimulus_level: N, threshold: N, width: N, gamma: N = 0, lambd: """ Calculate the slope at specified stimulus levels. Args: - prop_correct: Proportion correct at the threshold to calculate. + 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 @@ -101,7 +89,8 @@ def slope(self, stimulus_level: N, threshold: N, width: N, gamma: N = 0, lambd: Slope at the stimulus values. """ - slope = (1 - gamma - lambd) * self._slope(stimulus_level, threshold, width) + raw_slope = self._slope(stimulus_level, threshold, width) + slope = (1 - gamma - lambd) * raw_slope if self.negative: return -slope @@ -132,97 +121,178 @@ def inverse(self, prop_correct: N, threshold: N, width: N, prop_correct = 1 - prop_correct result = self._inverse(prop_correct, threshold, width) + return result - def _value(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - raise NotImplementedError("This should be overwritten by an implementation.") + def standard_parameters(self, threshold: N, width: N) -> tuple: + """ Transforms parameters threshold and width to a standard parametrization. + + The interpretation of the standard parameters, location and scale, depends on the sigmoid class used. + For instance, for a Gaussian sigmoid, the location corresponds to the mean and the scale to the standard + deviation of the distribution. + + For negative slope sigmoids, we return the same parameters as for the positive ones. + + Args: + threshold: Threshold value at PC + width: Width of the sigmoid + Returns: + Standard parameters (loc, scale) for the sigmoid subclass. + """ + return self._standard_parameters(threshold, width) + + # --- Private interface + + def _value(self, stimulus_level: N, threshold: N, width: N) -> N: + """ Compute the sigmoid value at specified stimulus levels. + + This method can be overwritten by subclasses to define new sigmoids. The base implementation uses the + CDF of a continuous distribution function defined in `scipy.stats` and returned by `_scipy_distr`. + + Args: + stimulus_level: Stimulus level values at which to calculate the sigmoid value + threshold: Threshold value at PC + width: Width of the sigmoid + Returns: + Proportion correct at the stimulus values. + """ + loc, scale = self._standard_parameters(threshold=threshold, width=width) + value = self._cdf(stimulus_level, loc=loc, scale=scale) + return value + + def _slope(self, stimulus_level: N, threshold: N, width: N) -> N: + """ Compute the slope of the sigmoid at a given stimulus level. + + This method can be overwritten by subclasses to define new sigmoids. The base implementation uses the + PDF of a continuous distribution function defined in `scipy.stats` and returned by `_scipy_distr`. - def _slope(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: + Args: + stimulus_level: Stimulus level value at which to calculate the slope + threshold: Threshold value at PC + width: Width of the sigmoid + Returns: + Slope at the stimulus level value + """ + loc, scale = self._standard_parameters(threshold=threshold, width=width) + raw_slope = self._pdf(stimulus_level, loc=loc, scale=scale) + return raw_slope + + def _inverse(self, prop_correct: N, threshold: N, width: N) -> N: + """ Compute the stimulus value at different proportion correct values. + + This method can be overwritten by subclasses to define new sigmoids. The base implementation uses the + PPF of a continuous distribution function defined in `scipy.stats` and returned by `_scipy_distr`. + + Args: + prop_correct: Proportion correct values at which to calculate the stimulus level values. + threshold: Threshold value at PC + width: Width of the sigmoid + Returns: + Stimulus values corresponding to the proportion correct values. + """ + loc, scale = self._standard_parameters(threshold=threshold, width=width) + result = self._ppf(prop_correct, loc=loc, scale=scale) + return result + + def _standard_parameters(self, threshold: N, width: N) -> list: + """ Transforms parameters threshold and width to a standard parametrization. + + This method must be overwritten by subclasses to define new sigmoids. + + The interpretation of the standard parameters, location and scale, depends on the sigmoid class used. + For instance, for a Gaussian sigmoid, the location corresponds to the mean and the scale to the standard + deviation of the distribution. + + For negative slope sigmoids, we return the same parameters as for the positive ones. + + Args: + threshold: Threshold value at PC + width: Width of the sigmoid + Returns: + Standard parameters (loc, scale) for the sigmoid subclass. + """ raise NotImplementedError("This should be overwritten by an implementation.") - def _inverse(self, prop_correct: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: + def _scipy_distr(self): + """ Get the scipy.stats implementation of the underlying cdf. + + Returns + scipy_distr: scipy.stats implementation of the distribution + distr_kwarg: A dictionary of extra keyword arguments needed to initialize the scipy.stats distribution + """ raise NotImplementedError("This should be overwritten by an implementation.") + def _ppf(self, x, loc=0, scale=1): + distr, distr_kwargs = self._scipy_distr() + return distr.ppf(x, loc=loc, scale=scale, **distr_kwargs) + + def _pdf(self, x, loc=0, scale=1): + distr, distr_kwargs = self._scipy_distr() + return distr.pdf(x, loc=loc, scale=scale, **distr_kwargs) + + def _cdf(self, x, loc=0, scale=1): + distr, distr_kwargs = self._scipy_distr() + return distr.cdf(x, loc=loc, scale=scale, **distr_kwargs) + class Gaussian(Sigmoid): """ Sigmoid based on the Gaussian distribution's CDF. """ - def _value(self, stimulus_level, threshold, width): - C =(norminv(1 - self.alpha) - norminv(self.alpha)) - return normcdf(stimulus_level, (threshold - norminvg(self._PC, 0, width / C)), width / C) - def _slope(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = norminv(1 - self.alpha) - norminv(self.alpha) - m = (threshold - norminvg(self._PC, 0, width / C)) - normalized_stimulus_level = (stimulus_level - m) / width * C - normalized_slope = sp.stats.norm.pdf(normalized_stimulus_level) - return normalized_slope * C / width + def _scipy_distr(self): + return scipy.stats.norm, {} + + def _standard_parameters(self, threshold: N, width: N) -> list: + scale = width / (self._ppf(1 - self.alpha) - self._ppf(self.alpha)) + loc = threshold - self._ppf(self._PC, 0, scale) + return [loc, scale] - def _inverse(self, prop_correct: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = norminv(1 - self.alpha) - norminv(self.alpha) - return norminvg(prop_correct, threshold - norminvg(self._PC, 0, width / C), width / C) class Logistic(Sigmoid): """ Sigmoid based on the Logistic distribution's CDF. """ - def _value(self, stimulus_level, threshold, width): - return 1 / (1 + np.exp(-2 * np.log(1 / self.alpha - 1) / width * (stimulus_level - threshold) - + np.log(1 / self._PC - 1))) - def _slope(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = 2 * np.log(1 / self.alpha - 1) / width - unscaled_slope = np.exp(-C * (stimulus_level - threshold) + np.log(1 / self._PC - 1)) - return C * unscaled_slope / (1 + unscaled_slope)**2 + def _scipy_distr(self): + return scipy.stats.logistic, {} - def _inverse(self, prop_correct: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - return (threshold - width * (np.log(1 / prop_correct - 1) - np.log(1 / self._PC - 1)) / 2 - / np.log(1 / self.alpha - 1)) + def _standard_parameters(self, threshold: N, width: N) -> list: + scale = width / (2 * np.log(1 / self.alpha - 1)) + loc = threshold - self._ppf(self._PC, loc=0, scale=scale) + return [loc, scale] class Gumbel(Sigmoid): """ Sigmoid based on the Gumbel distribution's CDF. """ - def _value(self, stimulus_level, threshold, width): - C = np.log(-np.log(self.alpha)) - np.log(-np.log(1 - self.alpha)) - return 1 - np.exp(-np.exp(C / width * (stimulus_level - threshold) + np.log(-np.log(1 - self._PC)))) - def _slope(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = np.log(-np.log(self.alpha)) - np.log(-np.log(1 - self.alpha)) - unscaled_stimulus_level = np.exp(C / width * (stimulus_level - threshold) + np.log(-np.log(1 - self._PC))) - return C / width * np.exp(-unscaled_stimulus_level) * unscaled_stimulus_level + def _scipy_distr(self): + return scipy.stats.gumbel_l, {} - def _inverse(self, prop_correct: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = np.log(-np.log(self.alpha)) - np.log(-np.log(1 - self.alpha)) - return threshold + (np.log(-np.log(1 - prop_correct)) - np.log(-np.log(1 - self._PC))) * width / C + def _standard_parameters(self, threshold: N, width: N) -> list: + scale = width / (np.log(-np.log(self.alpha)) - np.log(-np.log(1-self.alpha))) + loc = threshold - self._ppf(self._PC, loc=0, scale=scale) + return [loc, scale] class ReverseGumbel(Sigmoid): """ Sigmoid based on the reversed Gumbel distribution's CDF. """ - def _value(self, stimulus_level, threshold, width): - C = np.log(-np.log(1 - self.alpha)) - np.log(-np.log(self.alpha)) - return np.exp(-np.exp(C / width * (stimulus_level - threshold) + np.log(-np.log(self._PC)))) - def _slope(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = np.log(-np.log(1 - self.alpha)) - np.log(-np.log(self.alpha)) - unscaled_stimulus_level = np.exp(C / width * (stimulus_level - threshold) + np.log(-np.log(self._PC))) - return -C / width * np.exp(-unscaled_stimulus_level) * unscaled_stimulus_level + def _scipy_distr(self): + return scipy.stats.gumbel_r, {} - def _inverse(self, prop_correct: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = np.log(-np.log(1 - self.alpha)) - np.log(-np.log(self.alpha)) - return threshold + (np.log(-np.log(prop_correct)) - np.log(-np.log(self._PC))) * width / C + def _standard_parameters(self, threshold: N, width: N) -> list: + scale = width / (np.log(-np.log(self.alpha)) - np.log(-np.log(1-self.alpha))) + loc = threshold - self._ppf(self._PC, loc=0, scale=scale) + return [loc, scale] class Student(Sigmoid): """ Sigmoid based on the Student-t distribution's CDF. """ - def _value(self, stimulus_level, threshold, width): - C = (t1icdf(1 - self.alpha) - t1icdf(self.alpha)) - return t1cdf(C * (stimulus_level - threshold) / width + t1icdf(self._PC)) - - def _slope(self, stimulus_level: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = (t1icdf(1 - self.alpha) - t1icdf(self.alpha)) - stimLevel = (stimulus_level - threshold) * C / width + t1icdf(self._PC) - return C / width * sp.stats.t.pdf(stimLevel, df=1) - - def _inverse(self, prop_correct: np.ndarray, threshold: np.ndarray, width: np.ndarray) -> np.ndarray: - C = (t1icdf(1 - self.alpha) - t1icdf(self.alpha)) - return (t1icdf(prop_correct) - t1icdf(self._PC)) * width / C + threshold + + def _scipy_distr(self): + return scipy.stats.t, {'df': 1} + + def _standard_parameters(self, threshold: N, width: N) -> list: + scale = width / (self._ppf(1 - self.alpha) - self._ppf(self.alpha)) + loc = threshold - self._ppf(self._PC, loc=0, scale=scale) + return [loc, scale] class Weibull(Gumbel): diff --git a/psignifit/tests/test_result.py b/psignifit/tests/test_result.py index e53bc9a..f5deb71 100644 --- a/psignifit/tests/test_result.py +++ b/psignifit/tests/test_result.py @@ -185,3 +185,32 @@ def test_estimate_type_default(result): result.estimate_type = 'mean' estimate = result.get_parameters_estimate() assert _close_numpy_dict(estimate, result.parameters_estimate_mean) + + +def test_standard_parameters_estimate(): + width = 2.1 + threshold = 0.87 + parameter_estimate = { + 'threshold': threshold, + 'width': width, + 'lambda': 0.0, + 'gamma': 0.0, + 'eta': 0.0, + } + confidence_intervals = { + 'threshold': [[threshold, threshold]], + 'width': [[width, width]], + 'lambda': [[0.05, 0.2]], + 'gamma': [[0.1, 0.3]], + 'eta': [[0.0, 0.0]] + } + result = _build_result(parameter_estimate, parameter_estimate, confidence_intervals) + + # For a Gaussian sigmoid with alpha=0.05, PC=0.5 + expected_loc = threshold + # 1.644853626951472 is the normal PPF at alpha=0.95 + expected_scale = width / (2 * 1.644853626951472) + + loc, scale = result.standard_parameters_estimate() + np.testing.assert_allclose(loc, expected_loc) + np.testing.assert_allclose(scale, expected_scale) diff --git a/psignifit/tests/test_sigmoids.py b/psignifit/tests/test_sigmoids.py index aac4cd3..ca4fc69 100644 --- a/psignifit/tests/test_sigmoids.py +++ b/psignifit/tests/test_sigmoids.py @@ -1,4 +1,5 @@ import numpy as np +from scipy import stats import pytest from psignifit import sigmoids @@ -101,3 +102,36 @@ def test_sigmoid_sanity_check(sigmoid_name): sigmoids.assert_sigmoid_sanity_checks( sigmoid, n_samples=10000, threshold=threshold, width=0.7, ) + + +@pytest.mark.parametrize( + 'subclass, distr, distr_kwargs', + [ + (sigmoids.Gaussian, stats.norm, {}), + (sigmoids.Logistic, stats.logistic, {}), + (sigmoids.Gumbel, stats.gumbel_l, {}), + (sigmoids.ReverseGumbel, stats.gumbel_r, {}), + (sigmoids.Student, stats.t, {'df': 1}), + ] +) +def test_standard_parameters(subclass, distr, distr_kwargs): + PC = 0.4 + alpha = 0.083 + + # Create a sigmoid from standard parameters + original_loc = 3.1 + original_scale = 1.34 + original_distr = distr(loc=original_loc, scale=original_scale, **distr_kwargs) + + # Measure width and threshold + x = np.linspace(0, 6, 10000) + psi = original_distr.cdf(x) + threshold = x[np.argwhere(psi > PC)][0, 0] + width = original_distr.ppf(1 - alpha) - original_distr.ppf(alpha) + + # Check that the sigmoid method returns the original parameters + sigmoid = subclass(PC=PC, alpha=alpha) + loc, scale = sigmoid.standard_parameters(threshold=threshold, width=width) + + np.testing.assert_allclose(loc, original_loc, atol=1e-3) + np.testing.assert_allclose(scale, original_scale, atol=1e-3)