human
I am working on the file "Bayes.ipynb".
The first task is:
Now assume that we have data of failure moments of devices given by $D={10, 12, 15}$ and the question would be what the value of $\theta$ might be. Common sense will tell you that $\theta<10$, but probably also that $\theta>5$ or there would be more data close to 5. So what is it?
The details are, as you may have expected by now on a blog by Jake VanderPlas. In the interest of time we will not go into the nitty-gritty, but the result is outlined here.
Because of the low number of measurements, namely 3, the commonly employed "normal approximation", that is worked out in the given link, breaks down. This normal approximation happens to give nonsense results that have our intuitive answer outside the 95% confidence interval, which is [10.2 - 12.5].
For small $N$, the normal approximation will not apply, and we must instead compute the confidence integral from the actual sampling distribution, which is the distribution of the mean of $N$ variables each distributed according to $p(\theta)$. The sum of random variables is distributed according to the convolution of the distributions for individual variables, so we can exploit the convolution theorem and use the method of characteristic functions to find the following sampling distribution for the sum of $N$ variables distributed according to our particular $p(x~|~\theta)$:
$$
f(\theta~|~D) \propto
\left{
\begin{array}{lll}
z^{N - 1}\exp(-z) &,& z > 0\
0 &,& z < 0
\end{array}
\right}
;~ z = N(\hat{\theta} + 1 - \theta)
$$
To compute the 95% confidence interval, we can start by computing the cumulative distribution: we integrate $f(\theta~|~D)$ from $0$ to $\theta$ (note that we are not actually integrating over the parameter $\theta$, but over the estimate of $\theta$. Frequentists cannot integrate over parameters).
This integral is relatively painless if we make use of the expression for the incomplete gamma function:
$$
\Gamma(a, x) = \int_x^\infty t^{a - 1}e^{-t} dt
$$
which looks strikingly similar to our $f(\theta)$.
Using this to perform the integral, we find...
Can you write Python code for it?

=====
What is a Python code to find the following sampling distribution for the sum of $N$ variables distributed according to our particular $p(x~|~\theta)$:?
assistant
#Here is Python code to find the following sampling distribution for the sum of $N$ variables distributed according to our particular $p(x~|~\theta)$:
.

from scipy.special import gammaincc, erfinv
from scipy import optimize

# This approximate CI is the result of the much simpler normal approximation, 
# it serves as a simple helper function here, for the optimization down below.
def approx_CI(D, sig=0.95):
    """Approximate truncated exponential confidence interval"""
    # use erfinv to convert percentage to number of sigma
    Nsigma = np.sqrt(2) * erfinv(sig)
    D = np.asarray(D)
    N = D.size
    theta_hat = np.mean(D) - 1
    return [theta_hat - Nsigma / np.sqrt(N),
            theta_hat + Nsigma / np.sqrt(N)]

def exact_CI(D, frac=0.95):
    """Exact truncated exponential confidence interval"""
    D = np.asarray(D)
    N = D.size
    theta_hat = np.mean(D) - 1

    def f(theta, D):
        z = theta_hat + 1 - theta
        return (z > 0) * z ** (N - 1) * np.exp(-N * z)

    def F(theta, D):
        return gammaincc(N, np.maximum(0, N * (theta_hat + 1 - theta))) - gammaincc(N, N * (theta_hat + 1))
    
    def eqns(CI, D):
        """Equations which should be equal to zero"""
        theta1, theta2 = CI
        return (F(theta2, D) - F(theta1, D) - frac,
                f(theta2, D) - f(theta1, D))
    
    guess = approx_CI(D, 0.68) # use 1-sigma interval as a guess
    result = optimize.root(eqns, guess, args=(D,))
    if not result.success:
        print("warning: CI result did not converge!")
    return result.x

D = [10, 12, 15]
print("CI: ({0:.1f}, {1:.1f})".format(*exact_CI(D)))
print("The approximate CI using the normal approximation would give: ({0:.1f}, {1:.1f})".format(*approx_CI(D)))