Repository: wichmann-lab/python-psignifit

Issue:
Credible intervals incorrect for values below threshold?
# Context
Great work on psignifit. While reproducing some of my results, however, I unfortunately stumbled across a potential issue with the credible intervals. Below a code example (three versions) which shows that credible intervals for performance values below the threshold seem to be incorrect. 

# Minimal example
```
import numpy as np
import psignifit as ps
from psignifit import psigniplot
import matplotlib.pyplot as plt


plotVersion = 3  # select which of the proposed versions to plot (1-3)

# Example from documentation:
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'}

res = ps.psignifit(data, **options)

plt.figure()
psigniplot.plot_psychometric_function(res)


# Plot credible intervals across the whole curve
if plotVersion == 1: # produces an error
    y = np.linspace(0.5, 1, 50)  # get credible intervals between 50-100% performance
    CIlow = []; CIhigh = []
    for i in y:
        threshold, CI = res.threshold(i)  # use scaled sigmoid
        CIlow.append(CI[0][0]); CIhigh.append(CI[0][1])
    plt.plot(CIlow, y, label="lower bound"); plt.plot(CIhigh, y, label="upper bound")

elif plotVersion == 2:
    y = np.linspace(0.02, .5, 50) # get credible intervals according to error message of version 1
    CIlow = []; CIhigh = []
    for i in y:
        threshold, CI = res.threshold(i) # use scaled sigmoid
        CIlow.append(CI[0][0]); CIhigh.append(CI[0][1])
    plt.plot(CIlow, y+0.5, label="lower bound"); plt.plot(CIhigh, y+0.5, label="upper bound")

elif plotVersion == 3:
    y = np.linspace(0., 1, 50)  # get credible intervals over whole range for unscaled sigmoid
    CIlow = []; CIhigh = []
    for i in y:
        threshold, CI = res.threshold(i, unscaled=True) # use unscaled sigmoid instead
        CIlow.append(CI[0][0]); CIhigh.append(CI[0][1])
    
    y = np.linspace(0.5, 1, 50) # redefine to 50-100%
    plt.plot(CIlow, y, label="lower bound"); plt.plot(CIhigh, y, label="upper bound")

plt.legend()
plt.show()
```

# Expected behavior
I would have expected that the psychometric function is always inbetween the credible intervals across the whole performance-range.

# Observed behavior
The credible intervals for performance values below threshold (75% in the example) dont include the fitted function values.

## Version 2 of code
![image](https://github.com/user-attachments/assets/9deecaec-bd3c-4f7b-85cc-012685075de6)

## Version 3 of code
![image](https://github.com/user-attachments/assets/d2e61b14-0f6c-4928-8ac3-ee83a0f9f3dc)

Discussion:
@pberkes: Hi Lynn, thanks for the report. We'll investigate quickly.

Two things:
1) which version of `psignifit` are you using? We've been working on the library; ideally, you are using 4.2 (we will soon have a 4.3 release)
2) the confidence intervals are over all the parameters, not only the threshold; the two curves you are generating are the two sigmoids using the parameters at the extreme of the confidence intervals. In that sense, you won't necessarily find the optimal sigmoid inside the two other sigmoids. For example, when the CI of the threshold is tiny, and the CI of the width is large, then you would get two sigmoids that intersect in the same point but with very different slopes

In any case, we'll have a close look. In the latest version of the code, I'm getting to identical plots for versions 1 and 3, and an error for version 2, which I think it's correct because 2AFC should have proportion correct above 0.5@LynnSchmittwilken: Hi Pietro! Sounds great!
1. I've been using version 4.2; and I just reproduced the exact same plots with a fresh virtual environment when only installing psignifit 4.2
2. I see, that makes sense! However, I dont think that the upper and lower bound should ever reverse order regardless.

Since in your version of the code, the error does not exist for plotVersion==1 but rather plotVersion==2, it sounds like psignifit 4.3 might have already have fixed this issue (in part?). I agree, that in theory plotVersion==2 should throw an error. I'll just wait for the update then!@pberkes: 1. I'll check if there have been fixes to the CIs between 4.2 and 4.3dev
2. There might be a misunderstanding about the CIs here. If you think about having only a CI over the slope of the sigmoid, for example, then of course the curves would intersect.

For comparison, this is what happens if I fix all parameters except for the threshold. In that case, the optimal sigmoid lies between the two curves, just as you expect.

<img width="607" alt="image" src="https://github.com/user-attachments/assets/fa9377c9-ae74-4b0e-855f-bc9f2df65677" />

Code (v4.3dev):
```
import numpy as np
import psignifit as ps
from psignifit import psigniplot
import matplotlib.pyplot as plt

# Example from documentation:
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',
    'fixed_parameters': {
        'gamma': np.float64(0.5),
        'eta': np.float64(0.0),
        'lambda': np.float64(0.0),
        'width': np.float64(0.005)
    }
}

res = ps.psignifit(data, **options)

plt.figure()
psigniplot.plot_psychometric_function(res)

y = np.linspace(0.5, 0.98, 50)  # get credible intervals between 50-100% performance
threshold, CI = res.threshold(y)  # use scaled sigmoid
plt.plot(CI['0.95'][0], y, label="lower bound")
plt.plot(CI['0.95'][1], y, label="upper bound")
plt.legend()
```@pberkes: and if I fix everything but the width

<img width="632" alt="image" src="https://github.com/user-attachments/assets/1c968467-958d-4d91-8aa1-0a0164634133" />
@LynnSchmittwilken: Great, this behavior looks good to me!

So, to me, it sounds like 4.3 is fixing a (potentially minor) issue in the credible intervals related to why in 4.2 plotVersion==1 in my example code throws an error.
And then, it seems like there might just be a difference in default behavior as compared to the pre-pip version of psignifit. The old psignifit might have kept all parameters except for the threshold constant by default when using the "threshold"-method to compute the credible intervals, but this I am not sure about without checking the code more thoroughly.

Happy to give this another look, once 4.3 is out! @pberkes: My investigative report:
- In 4.2 there is a bug where `Result.threshold` calls `inverse` with the gamma and lambda parameters inverted. This was fixed in #139 . This is why plot version 2 works even though plot version 1 is the correct one
- I believe the output of the method is correct, even though I have been unable to explain clearly why ... :-)  the blue and orange curves are not bounds on the optimal sigmoids, they are the two most extreme sigmoids that you could get when you set the parameters at the extreme of their confidence intervals. Maybe @guillermoaguilar can jump in and help me out? :-)
- The big difference in the curve between version 2 and 3 is due to the same #139 bug, otherwise they should be identical if you rescale the input variable `y`
@guillermoaguilar: I'm happy to see engagement here, @LynnSchmittwilken and @pberkes ! :D

> * I believe the output of the method is correct, even though I have been unable to explain clearly why ... :-)  the blue and orange curves are not bounds on the optimal sigmoids, they are the two most extreme sigmoids that you could get when you set the parameters at the extreme of their confidence intervals. Maybe @guillermoaguilar can jump in and help me out? :-)

Indeed I think the method is correct, which is maybe counter-intuitive. Maybe another way of understanding is by looking at the posterior of parameters `threshold` and `width`.  You run `psignifit` with argument `debug=True`, and then do the plot with `psigniplot.plot_2D_margin(res, 'threshold', 'width')`.

For the completely unbounded case, no fixed parameters, the 2D posterior has a mass that is *correlated*, meaning a blob that is not just circular but tilted

![issue_lynn_2dmarginal1](https://github.com/user-attachments/assets/9f6c5315-0aaf-444a-a010-80ffa61706e6)


The CIs are computed from the projection from the high dimensional posterior (in this case 4-D) into each parameter axes. Here we visualized only 2 of those axes, there are two more. The observed behavior (the crossings in the functions) can happen when the mass is like this, because there will be values at the boundary of the blob which is inside the interval for one parameter but outside for the other.

_In other words the range in CI from the posterior is telling you that *with your data all those possible shapes of sigmoids are possible within confidence*, all that family of sigmoids shifted horizontally but also with higher or lower slope._  
 
As a counter-example  when all parameters are fixed but the width (the case of @pberkes), you have the following 2D marginals:

![issue_lynn_2dmarginal_widthfixed](https://github.com/user-attachments/assets/b06f3bb3-839f-41a4-8fc6-57bbc098a25b)

which is trivial as the width is fixed, so the mass is aligned to the axes. And here the CI follows the intuition for the range of possible sigmoids.

More generally, CIs are thought to give you the interval *for the fitted parameter, not the whole function.*

@guillermoaguilar: So I'd say this is not an issue per se.@guillermoaguilar: On second thought, maybe I'm wrong and it doesn't have to do with the correlation in the posterior but just from the fact that here you are estimating threshold and width, both are free, and of course a sigmoid on the lower 'bound' for the width parameter will be shallower than the fitted, and then the function will be above the fitted for values less than the threshold. And the other way around, like here:

<img width="632" alt="397814648-1c968467-958d-4d91-8aa1-0a0164634133" src="https://github.com/user-attachments/assets/8d937ac0-21b6-4b35-b144-51dc46f2cb2e" />

So maybe my explanation above is not the reason why this happens, it's just what happens when you have both threshold and width free to vary.
In any case the method seem to be correct.


@LynnSchmittwilken : if you want to visualize the range of functions that are 'compatible' with your data, I'd generate ones corresponding to varying both width and threshold simultaneously.
@pberkes: Reporting @FelixWichmann 's reply by email:
"""
Well … I certainly sympathise with Lynn here ;-)

The plot of the psychometric function should really show the CIs, not two samples from the posterior -- thus I'd expect as a plot:
a. PF and CIs around the PF. Just a single PF and "error bar" (i.e. without blue & orange PFs)
b. PF and the CI as "sausage" around the PF. This can be obtained by sampling a few hundred PFs from the posterior: The fitted PF will *never* be outside the "sausage" of the random samples from the posterior (once you have enough samples).

In Lynn's plot the light blue and orange are not the lower and upper bound but two specific samples from the posterior with the threshold at the lower/upper and the slope at the upper/lower bound. To me this is misleading---as it was for Lynn---because I'd initially interpreted this as the bounds of the "sausage" (see b. above) ... and then it would be wrong!
If you really want to show a "quick and dirty" visualisation of the "sausage" you should draw four (!) half-PFs:
1. lower bound for threshold and lower bound for slope for proportion correct < 0.5 (nominal)
2. lower bound for threshold and upper bound for slope for proportion correct > 0.5 (nominal)
3. upper bound for threshold and lower bound for slope for proportion correct > 0.5 (nominal)
4. upper bound for threshold and upper bound for slope for proportion correct < 0.5 (nominal)

In the figure only the half-PFs for p > 0.5 are correct, and then you have to swap the slopes below the threshold ("fuse" the orange lower half to the upper blue half and vice versa).

BUT: I'd stick to either a.) or b.) above.
"""@pberkes: @FelixWichmann I can't really disagree... but the function works as advertized, don't you agree  ?

This code implements Felix' proposal b. :


```
options = {
    'sigmoid': 'norm',
    'experiment_type' : '2AFC',
    'fixed_parameters': {
        'eta': np.float64(0.0),
        'lambda': np.float64(0.0),
        'threshold': np.float64(0.005)
    }
}

res = ps.psignifit(data, debug=True, **options)


params_grid = np.meshgrid(*(res.parameter_values.values()), indexing='ij')
params_stack = np.dstack([p.flatten() for p in params_grid]).squeeze()

# Sample from the posterior
samples_idx = np.random.choice(np.arange(params_stack.shape[0]), size=(100,), replace=True, p=res.debug['posteriors'].flatten())
params_samples = params_stack[samples_idx]

# --- Plot
plt.figure()

# Plot the samples from the posterior
sigmoid = res.configuration.make_sigmoid()
x = np.linspace(0.001, 0.01, num=1000)
for idx in range(params_samples.shape[0]):
    y = sigmoid(x, threshold=params_samples[idx, 3], width=params_samples[idx, 4], gamma=params_samples[idx, 1], lambd=params_samples[idx, 2])
    plt.plot(x, y, alpha=0.1, color='gray')

# Plot the MAP estimate
psigniplot.plot_psychometric_function(res, line_width=3)
```

Which gives this result (100 samples from the posterior):
<img width="598" alt="image" src="https://github.com/user-attachments/assets/f4d02af0-2ebb-491f-9505-0c56d5f94b7d" />
@pberkes: For the case with no fixed parameters, the plot looks like this
<img width="597" alt="image" src="https://github.com/user-attachments/assets/7e4d1e87-cbfc-4fda-8fd4-f7e6082e240b" />
@guillermoaguilar: thanks for clarifying Felix,  @FelixWichmann , and thanks @pberkes for implementing the plots, they're really good and beautiful as well.

I think we should definitely include this type of plot to the package, in `psigniplot.py`, and add it to the plotting demos in the documentation. It is really clear to show what is the range of functions for the given confidence level.@LynnSchmittwilken: At this point, I am really not sure whether my original issue "just" stems from the bug in version 4.2. or whether there is also a more conceptual misunderstanding / disagreement.
I'd be happy to wait until 4.3. and give this another look, maybe with my reproduction attempts of what the old "threshold-method" did (pre-pip psignifit) in comparison to the new version. If there is actually a conceptual difference between the two, maybe having a real example helps to make a more informed decision about what the best default behavior of that method should be, and how to get from one to the other (ie which variables need to be changed)

Thanks already for all the effort and quick reactions! @pberkes: There is a consensus that the method is working correctly (in v4.3), with a doubt if the two extreme samples returned by the current method are useful (the current method being picking the parameter values at the extreme of the confidence interval, for all parameters independently).

My suggestion:
1. Remove the `return_ci` option from `Result.threshold`, it's not very useful
2. Add a `Result.posterior_samples` method to sample parameters from the posterior
3. Add a method to plot posterior samples

@pberkes: The only argument I can think of for keeping the `return_ci` is that it does not require the sigmoid to be fit with `debug=True`@guillermoaguilar: > There is a consensus that the method is working correctly (in v4.3), with a doubt if the two extreme samples returned by the current method are useful (the current method being picking the parameter values at the extreme of the confidence interval, for all parameters independently).
> 
> My suggestion:
> 
> 1. Remove the `return_ci` option from `Result.threshold`, it's not very useful
> 2. Add a `Result.posterior_samples` method to sample parameters from the posterior
> 3. Add a method to plot posterior samples

Agree with 2 and 3 for sure! And this is done in PR #196 

But I wouldn't remove the option `return_ci`, point 1. 
After fitting the psychometric function some users just want to get the threshold and its CI.  They will do this for different experimental conditions and compare them. So I def would keep the functionality that it returns the CI.@pberkes: @guillermoaguilar agreed, let's keep the CI option

Closing this issue as "fixed", partly in #139 and partly in #196 

@LynnSchmittwilken please do give 4.3 a run when we release it, and please reopen this issue if you see any remaining issues@guillermoaguilar: I'm reopening this issue as I detected an inconsistency between the new refactored version and the old python version (and MATLAB version).

The inconsistency is in the threshold method of the result object, specifically in how the confidence intervals are calculated for values away from the default threshold.  In the MATLAB code as well as in the python version 4.0, the method uses different bounds depending on whether the requested percentage correct is above or below the percentage at which the threshold is defined (variable thresh_PC). See [here for MATLAB code](https://github.com/wichmann-lab/psignifit/blob/5f7241f6288c7dd3e1b34353118ca919ad28e51b/getThreshold.m#L113) and [here for python code](https://github.com/wichmann-lab/python-psignifit/blob/68d54aa35951930de73f5644b7bdf4c97e0f8953/psignifit/psignifit.py#L725).

In the current implementation there is no distinction made. So the lowest CI is always drawn from the low bound for all parameters... that includes the width. For the width parameter, the low bound is obvs. the steepest function, thus the bound lines will always intercept somewhere. That does not reflect what we want to show, the worse-case confidence interval. 

The correct way should be: 

for stimulus values above thresh_PC:
- the low bound returned should be drawn from the *low* bound for the threshold, the *low* bound for the width, the *low* bound for gamma, and the *high* bound for lambda. 
- the high bound returned should be drawn from the *high* bound for the threshold, the *high* bound for the width, the *high* bound for gamma, and the *low* bound for lambda. 

for stimulus values below  thresh_PC:

- the low bound returned should be drawn from the *low* bound for the threshold, the *high* bound for the width, the *low* bound for gamma, and the *high* bound for lambda. 
- the high bound returned should be drawn from the *high* bound for the threshold, the *low* bound for the width, the *high* bound for gamma, and the *low* bound for lambda. 

This logic is also what Felix mentioned in his response, [here above](https://github.com/wichmann-lab/python-psignifit/issues/191#issuecomment-2557726367)

Note: to calculate CIs like this is still the "worse" , upper-bound confidence interval. The more accurate way to get a threshold at another level is to refit the model and changing the thresh_PC parameter. And to visualize the uncertainty, we have the [posterior sample plot](https://psignifit.readthedocs.io/en/latest/examples/plotting.html#plot-samples-from-the-posterior-over-psychometric-functions). I think we should  include this information as warning, MATLAB version has a warning about it.
@pberkes: ![Image](https://github.com/user-attachments/assets/739819ea-f544-4017-8c7b-0d62b1210702)

@guillermoaguilar if I understand what you want to do, in the image above you'd like the top two lines to both be orange, and the bottom ones to be blue.

I remain unconvinced by this -- we're plotting upper and lower bounds **in the space of psychometric functions**. Your proposal would suggests that the upper and lower bounds are in the space of proportion correct. This interpretation would suggest that most observed data points should lie between the upper and lower bound, but that's not the case: the upper and lower bounds on the psychometric functions have an additional noise term around them. (but I think you agree with this)

If your and Felix' proposal is what is expected in general in the field, I'm not going to insist.@pberkes: To illustrate a bit more my argument, let's say that the psychometric function is parametrized only by threshold and width.

In the space of this parameters, we can compute CIs for the parameters that characterize the most extreme functions, as in the plot below: 

<img width="894" alt="Image" src="https://github.com/user-attachments/assets/9559a6d5-52f9-4279-8903-45f691dc52b7" />

These bounds correspond directly to the two "extreme" functions that we are plotting at the moment, when we plot them in the space of stimulus level vs proportion correct.

If I understand correctly, you'd like to plot this

<img width="1050" alt="Image" src="https://github.com/user-attachments/assets/6b5af043-762e-4adc-a8fa-9b17b249a65e" />

I don't think you're suggesting that the orange and blue lines are confidence intervals over observed data, but then they don't correspond to much else either, except to some intuitive notion that "upper" should be up and "lower" should be down.

My apologies if I'm missing the point and being thick


@guillermoaguilar: @pberkes  no worries, this requires also some mind bending for me too.  

You're right on the heatmap plot, no question about that, those are the bounds in parameter space. The issue is that users want to know what are the bounds of uncertainty in "psychometric function space" too. For that we have the posterior samples plot, and that is a nice visualization. But say someone wants to get the CI for a percentage correct at 0.6, or 0.85. The question is then how do we return those CIs. 

What Felix +Heiko did in their implementation was to return an coarse approximation for that CI, by using the logic described above:
For the low bound  x-value returned by method threshold(): 
- Below the inflection point (thresh_PC) you need to use the high bound for the width and the low bound for the threshold parameter, because that gives you the shallowest function, and thus the lowest possible x-values (the x-value more to the left). 
- Above the inflection point you use the low bound for the width and the low bound for the threshold parameter, because that gives you the steepest function more to the left, and thus also the leftmost possible x-values. 

The other way around for the high bound x-value  returned by method threshold()

These are for sure not accurate and too broad, but at least it gives you something. [ Heiko+Felix in the WIKI for the MATLAB version write](https://github.com/wichmann-lab/psignifit/wiki/How-to-Get-Thresholds-and-Slopes#obtaining-threshold-values):

> Important note: The credible intervals psignifit 4 calculates are based on the proportion correct level specified in your options struct (default is half-way up the sigmoid; information regarding the change of this value can be found [here](https://github.com/wichmann-lab/psignifit/wiki/How-to-Change-the-Threshold-Percent-Correct)). The credible intervals calculated by this function post Bayesian analysis are conservative: The further you move away from the original threshold, the more conservative they become as we simply and independently assume the worst case for all other parameters. For an accurate analysis first [change the threshold](https://github.com/wichmann-lab/psignifit/wiki/How-to-Change-the-Threshold-Percent-Correct) to the desired value and then re-run psignifit 4 using this value.


A similar warning I've included in PR #241. 

Now visually, what we have now is this

![Image](https://github.com/user-attachments/assets/c7df5e9f-71d9-4446-86e3-d11a92f1999e)

which is OK above the inflection point but not below. It should give something like this:

![Image](https://github.com/user-attachments/assets/1e83dcf0-0e07-4190-b607-61613f0daff2)


The first plot was done with this code using branch main, the second one using the branch from my (incomplete) PR #241 

```python
import matplotlib.pyplot as plt
import numpy as np
import psignifit as ps
import psignifit.psigniplot as psp

# Example from documentation:
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]])

fixed_parameters = {'eta': 0.0,
                    'lambda': 0.0,
                    }

result = ps.psignifit(data,
                      experiment_type='2AFC',
                      fixed_parameters=fixed_parameters)

# create x and y vector of psy. function
x = np.linspace(data[:, 0].min(), data[:, 0].max(), num=1000)
y = result.proportion_correct(x)

# %%  gets CIs at the different levels of percentage correct
CIu, CIl = [], []
for i, yi in enumerate(y):
    _, CI = result.threshold(yi)
    CIl.append(CI['0.95'][0])
    CIu.append(CI['0.95'][1])

# %%
fig, axes = plt.subplots(1, 1, figsize=(5, 4))
psp.plot_psychometric_function(result,
                               ax=axes,
                               data_color='k',
                               line_width=1,
                               data_size=0.1,
                               )
axes.fill_betweenx(y, np.array(CIl), np.array(CIu), color='g', alpha=0.2)
fig.savefig('envelope.png', dpi=200)
plt.show()`
```


Again these are CIs that are expected to be broader than necessary (upper-bound for the interval). Still people ask what is the CI for the x-value at a PC different than the default. 

Also because as it is, it doesnt make sense the CI values returned  when you ask for, say, pc =0.55
```
print(result.threshold(0.55)[1]['0.95'])
> [0.0028785460332117234, 0.0026255611358453864]
```
it doesnt make sense that the lower bound > higher bound ... 


I hope it's clearer now.... :D and thank you for reading! @guillermoaguilar: It just occurred to me an alternative....  to use the posterior samples here to return the requested CI. 

So maybe inside the method threshold(), if CIs are requested for a pc not the default, then the method runs posterior samples (a high number say 1000), evaluates the psy. function on those samples, and returns the percentile CI  calculated over the samples for the x-value at the requested percentage correct level.  That would convey a more adequate representation of the underlying uncertainty, no? @pberkes: @guillermoaguilar thanks for the careful explanation, I understand the argument now. We are not computing a real CI over the threshold, but this is what the warning is about.

The alternative way to compute the CI using the samples would marginalize over all the other parameters and indeed it could give an accurate (if non-deterministic) estimate of the CI over the threshold.

I don't have objections with going forward with the implementation in #241. Should we complete that PR and then discuss how to implement a sampling-based CI?@guillermoaguilar: yes let's do that! 

Solution patch:
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 dataclasses
 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 @@ class Result:
             (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 @@ class Result:
             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
 


Test patch:
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():