Hmmm... the only reason all the reshaping crap at the top is even necessary is that we index the output arrays on L378-379. What if we just do the test on line 378:
```python
if lowerbound.ndim == 0 # lowerbound and upperbound should have same ndim
    if k == 0:
        lowerbound = 0.
    elif k == n:
        upperbound = 1.
else:
    lowerbound[k == 0] = 0.
    upperbound[k == n] = 1.
```
Then I think we could entirely remove all the reshaping stuff at the top and just return `conf_interval` at the end.