Skip to content

Commit

Permalink
Merge pull request wichmann-lab#178 from pberkes/fix-threshold
Browse files Browse the repository at this point in the history 
Fix computation of confidence interval `Result.threshold` method
  • Loading branch information
@otizonaizit
otizonaizit authored Dec 4, 2024
2 parents 87eaa0e + 0b0d0eb commit 4c1457f
Show file tree
Hide file tree
Showing 2 changed files with 112 additions and 36 deletions.
21 changes: 12 additions & 9 deletions psignifit/_result.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,18 +126,21 @@ def threshold(self, proportion_correct: np.ndarray, unscaled: bool = False, retu
if not return_ci:
return new_threshold

param_cis = [self.confidence_intervals[param] for param in ('threshold', 'width', 'lambda', 'gamma')]
if unscaled: # set asymptotes to 0
param_cis[2] = np.zeros_like(param_cis[2])
param_cis[3] = np.zeros_like(param_cis[3])

new_ci = []
for (thres_ci, width_ci, lambd_ci, gamma_ci) in zip(*param_cis):
new_threshold_ci = {}
for coverage_key in self.confidence_intervals['threshold'].keys():
thres_ci = self.confidence_intervals['threshold'][coverage_key]
width_ci = self.confidence_intervals['width'][coverage_key]
if unscaled:
gamma_ci = np.array([0.0, 0.0])
lambd_ci = np.array([0.0, 0.0])
else:
gamma_ci = self.confidence_intervals['gamma'][coverage_key]
lambd_ci = self.confidence_intervals['lambda'][coverage_key]
ci_min = sigmoid.inverse(proportion_correct, thres_ci[0], width_ci[0], gamma_ci[0], lambd_ci[0])
ci_max = sigmoid.inverse(proportion_correct, thres_ci[1], width_ci[1], gamma_ci[1], lambd_ci[1])
new_ci.append([ci_min, ci_max])
new_threshold_ci[coverage_key] = [ci_min, ci_max]

return new_threshold, new_ci
return new_threshold, new_threshold_ci

def slope(self, stimulus_level: np.ndarray, estimate_type: Optional[EstimateType]=None) -> np.ndarray:
""" Slope of the psychometric function at a given stimulus levels.
Expand Down
127 changes: 100 additions & 27 deletions psignifit/tests/test_result.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,32 +3,31 @@
import numpy as np
import pytest

from psignifit import Configuration
from psignifit import Result
from psignifit import Configuration, Result, psignifit


@pytest.fixture
def result():
parameter_estimate_MAP = {
'threshold': 0.005,
'threshold': 1.005,
'width': 0.005,
'lambda': 0.0123,
'gamma': 0.021,
'eta': 0.0001
}
parameter_estimate_mean = {
'threshold': 0.002,
'threshold': 1.002,
'width': 0.002,
'lambda': 0.013,
'gamma': 0.024,
'eta': 0.001
}
confidence_intervals = {
'threshold': [[0.001, 0.005], [0.005, 0.01], [0.1, 0.2]],
'width': [[0.001, 0.005], [0.005, 0.01], [0.1, 0.2]],
'lambda': [[0.001, 0.005], [0.005, 0.01], [0.1, 0.2]],
'gamma': [[0.001, 0.005], [0.005, 0.01], [0.1, 0.2]],
'eta': [[0.001, 0.005], [0.005, 0.01], [0.1, 0.2]]
'threshold': {'0.95': [1.001, 1.005], '0.9': [1.005, 1.01]},
'width': {'0.95': [0.001, 0.005], '0.9': [0.005, 0.01]},
'lambda': {'0.95': [0.03, 0.08], '0.9': [0.1, 0.2]},
'gamma': {'0.95': [0.03, 0.08], '0.9': [0.05, 0.2]},
'eta': {'0.95': [0.001, 0.005], '0.9': [0.005, 0.01]},
}
return _build_result(parameter_estimate_MAP, parameter_estimate_mean, confidence_intervals)

Expand Down Expand Up @@ -92,13 +91,63 @@ def test_threshold_raises_error_when_outside_valid_range(result):

def test_threshold_slope(result):
proportion_correct = np.linspace(0.2, 0.5, num=1000)
stimulus_levels, confidence_intervals = result.threshold(proportion_correct)
stimulus_levels = result.threshold(proportion_correct, return_ci=False)
np.testing.assert_allclose(
result.slope(stimulus_levels),
result.slope_at_proportion_correct(proportion_correct)
)


def test_threshold_slope_ci_scaled(result):
proportion_correct = [0.4, 0.5, 0.7]
_, threshold_cis = result.threshold(proportion_correct, return_ci=True, unscaled=False)

expected = {
'0.95': [[1.000918, 1.001, 1.001171], [1.00454 , 1.005, 1.005969]],
'0.9': [[1.004661, 1.005112, 1.006097], [1.008691, 1.01, 1.012941]],
}
assert list(threshold_cis.keys()) == ['0.95', '0.9']
for coverage_key, cis in threshold_cis.items():
# one CI per proportion_correct
assert threshold_cis[coverage_key][0].shape[0] == len(proportion_correct)
assert threshold_cis[coverage_key][1].shape[0] == len(proportion_correct)
np.testing.assert_allclose(
threshold_cis[coverage_key][0],
expected[coverage_key][0],
atol=1e-6
)
np.testing.assert_allclose(
threshold_cis[coverage_key][1],
expected[coverage_key][1],
atol=1e-6
)


def test_threshold_slope_ci_unscaled(result):
proportion_correct = [0.4, 0.5, 0.7]
_, threshold_cis = result.threshold(proportion_correct, return_ci=True, unscaled=True)

expected = {
'0.95': [[1.000923, 1.001, 1.001159], [1.004615, 1.005, 1.005797]],
'0.9': [[1.004615, 1.005, 1.005797], [1.00923, 1.01, 1.011594]],
}
assert list(threshold_cis.keys()) == ['0.95', '0.9']
for coverage_key, cis in threshold_cis.items():
# one CI per proportion_correct
assert threshold_cis[coverage_key][0].shape[0] == len(proportion_correct)
assert threshold_cis[coverage_key][1].shape[0] == len(proportion_correct)
np.testing.assert_allclose(
threshold_cis[coverage_key][0],
expected[coverage_key][0],
atol=1e-6
)
np.testing.assert_allclose(
threshold_cis[coverage_key][1],
expected[coverage_key][1],
atol=1e-6
)


def test_threshold_value():
# This test fails before PR #139
lambda_ = 0.1
Expand All @@ -112,21 +161,28 @@ def test_threshold_value():
'eta': 0.0,
}
confidence_intervals = {
'threshold': [[0.5, 0.5]],
'width': [[width, width]],
'lambda': [[0.05, 0.2]],
'gamma': [[0.1, 0.3]],
'eta': [[0.0, 0.0]]
'threshold': {'0.95': [0.5, 0.5]},
'width': {'0.95': [width, width]},
'lambda': {'0.95': [0.05, 0.2]},
'gamma': {'0.95': [0.1, 0.3]},
'eta': {'0.95': [0.0, 0.0]},
}
result = _build_result(parameter_estimate, parameter_estimate, confidence_intervals)

# The threshold at the middle of the gamma-to-(1-lambda) range must be 0.5 for a Gaussian
thr, thr_ci = result.threshold(
thr = result.threshold(
proportion_correct=np.array([(1 - lambda_ - gamma) / 2.0 + gamma]),
unscaled=False, return_ci=True,
unscaled=False, return_ci=False,
)
expected_thr = np.array(0.5) # by construction
np.testing.assert_allclose(thr, expected_thr)

expected_thr = np.array(0.5)
# ... except when unscaled is True
thr = result.threshold(
proportion_correct=np.array([(1 - lambda_ - gamma) / 2.0 + gamma]),
unscaled=True, return_ci=False,
)
expected_thr = np.array([0.5343785]) # computed by hand
np.testing.assert_allclose(thr, expected_thr)

# Compare to results computed by hand
Expand All @@ -136,8 +192,8 @@ def test_threshold_value():
)
expected_thr = np.array([0.654833]) # Computed by hand
np.testing.assert_allclose(thr, expected_thr, atol=1e-4)
expected_thr_ci = [[[0.648115], [0.730251]]] # Computed by hand
np.testing.assert_allclose(thr_ci, expected_thr_ci, atol=1e-4)
expected_thr_ci = {'0.95': [[0.648115], [0.730251]]} # Computed by hand
np.testing.assert_allclose(thr_ci['0.95'], expected_thr_ci['0.95'], atol=1e-4)


def _close_numpy_dict(first, second):
Expand Down Expand Up @@ -197,13 +253,7 @@ def test_standard_parameters_estimate():
'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]]
}
confidence_intervals = {}
result = _build_result(parameter_estimate, parameter_estimate, confidence_intervals)

# For a Gaussian sigmoid with alpha=0.05, PC=0.5
Expand All @@ -214,3 +264,26 @@ def test_standard_parameters_estimate():
loc, scale = result.standard_parameters_estimate()
np.testing.assert_allclose(loc, expected_loc)
np.testing.assert_allclose(scale, expected_scale)


def test_threshold_bug_172():
# Reproduce bug in issue #172

data = np.array([[0.0010, 45.0000, 90.0000], [0.0015, 50.0000, 90.0000],
[0.0020, 44.0000, 90.0000], [0.0025, 44.0000, 90.0000],
[0.0030, 52.0000, 90.0000], [0.0035, 53.0000, 90.0000],
[0.0040, 62.0000, 90.0000], [0.0045, 64.0000, 90.0000],
[0.0050, 76.0000, 90.0000], [0.0060, 79.0000, 90.0000],
[0.0070, 88.0000, 90.0000], [0.0080, 90.0000, 90.0000],
[0.0100, 90.0000, 90.0000]])

options = {
'sigmoid': 'norm',
'experiment_type': '2AFC'
}

result = psignifit(data, **options)
threshold = result.threshold(0.9, return_ci=False) # which should be 0.0058

expected_threshold = 0.0058
np.testing.assert_allclose(threshold, expected_threshold, atol=1e-4)

0 comments on commit 4c1457f

Please sign in to comment.