diff --git a/psignifit/_result.py b/psignifit/_result.py
index 3655a04..eb1f743 100644
--- a/psignifit/_result.py
+++ b/psignifit/_result.py
@@ -2,12 +2,14 @@
 import json
 from typing import Any, Dict, Tuple, List, Optional, TextIO, Union
 from pathlib import Path
+import warnings
 
 import numpy as np
 from numpy.typing import NDArray
 
 from ._configuration import Configuration
 from ._typing import EstimateType
+from ._utils import cast_np_scalar
 from .tools import psychometric_with_eta
 
 class NumpyEncoder(json.JSONEncoder):
@@ -115,14 +117,29 @@ def threshold(self, proportion_correct: np.ndarray, unscaled: bool = False, retu
             (thresholds, ci): stimulus values along with confidence intervals
 
         """
+        if return_ci:
+            warnings.warn("""The confidence intervals computed by this method are only upper bounds. 
+                          To get a more accurate confidence interval at a level of proportion
+                          correct other than `thresh_PC`, you need to redefine the threshold, that is, redefine at
+                          which level the threshold parameter is set. You can do that by changing the argument
+                          'thresh_PC', and call psignifit again. For an example, see documentation,
+                          page "Advanced options").
+                          
+                          If instead you want to get the range of uncertainty for the psychometric 
+                          function fit, then you need to sample from the posterior and visualize those 
+                          samples. The function `plot_posterior_samples` in psigniplot does exactly that.
+                          You find an example of that visualization in the documentation, page Plotting.""")
+
         proportion_correct = np.asarray(proportion_correct)
         sigmoid = self.configuration.make_sigmoid()
 
         estimate = self.get_parameter_estimate(estimate_type)
         if unscaled:  # set asymptotes to 0 for everything.
             lambd, gamma = 0, 0
+            proportion_correct_unscaled = proportion_correct
         else:
             lambd, gamma = estimate['lambda'], estimate['gamma']
+            proportion_correct_unscaled = (proportion_correct - gamma) / (1- lambd - gamma)
         new_threshold = sigmoid.inverse(proportion_correct, estimate['threshold'],
                                         estimate['width'], gamma, lambd)
         if not return_ci:
@@ -138,9 +155,24 @@ def threshold(self, proportion_correct: np.ndarray, unscaled: bool = False, retu
             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_threshold_ci[coverage_key] = [ci_min, ci_max]
+
+            mask_above = proportion_correct_unscaled > self.configuration.thresh_PC
+
+            ci_min = np.zeros(proportion_correct.shape)
+            ci_max = np.zeros(proportion_correct.shape)
+
+            ci_min[mask_above] = sigmoid.inverse(proportion_correct[mask_above],
+                                                 thres_ci[0], width_ci[0], gamma_ci[0], lambd_ci[0])
+            ci_max[mask_above] = sigmoid.inverse(proportion_correct[mask_above],
+                                                 thres_ci[1], width_ci[1], gamma_ci[1], lambd_ci[1])
+
+            ci_min[~mask_above] = sigmoid.inverse(proportion_correct[~mask_above],
+                                                 thres_ci[0], width_ci[1], gamma_ci[0], lambd_ci[0])
+            ci_max[~mask_above] = sigmoid.inverse(proportion_correct[~mask_above],
+                                                 thres_ci[1], width_ci[0], gamma_ci[1], lambd_ci[1])
+
+            new_threshold_ci[coverage_key] = [cast_np_scalar(ci_min),
+                                              cast_np_scalar(ci_max)]
 
         return new_threshold, new_threshold_ci
 
diff --git a/tests/test_result.py b/tests/test_result.py
index 7ff01bf..4933fe9 100644
--- a/tests/test_result.py
+++ b/tests/test_result.py
@@ -104,8 +104,8 @@ def test_threshold_slope_ci_scaled(result):
     _, 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]],
+        '0.95': [[1.00059, 1.001, 1.001171], [1.004908 , 1.005, 1.005969]],
+        '0.9': [[1.004322, 1.005224, 1.006097], [1.009345, 1.01, 1.012941]],
     }
     assert list(threshold_cis.keys()) == ['0.95', '0.9']
     for coverage_key, cis in threshold_cis.items():
@@ -129,8 +129,8 @@ def test_threshold_slope_ci_unscaled(result):
     _, 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]],
+        '0.95': [[1.000615, 1.001, 1.001159], [1.004923, 1.005, 1.005797]],
+        '0.9': [[1.00423, 1.005, 1.005797], [1.009615, 1.01, 1.011594]],
     }
     assert list(threshold_cis.keys()) == ['0.95', '0.9']
     for coverage_key, cis in threshold_cis.items():