human
I am working on the file "Sessions/Session04/Day1/LSSTC-DSFP4-Juric-FrequentistAndBayes-03-Credibility.ipynb".
The first task is:
We immediately see an issue. By our simple common sense argument, we've determined that it is impossible for $\theta$ to be greater than 10, yet the entirety of the 95% confidence interval is above this range!  Perhaps this issue is due to the small sample size: the above computation is based on a large-$N$ approximation, and we have a relatively paltry $N = 3$.
Maybe this will be improved if we do the more computationally intensive exact approach?
The answer is no. If we compute the confidence interval without relying on large-$N$ Gaussian eapproximation, the result is $(10.2, 12.2)$.
Note: you can verify yourself by evaluating the code in the sub-slides.
3. Frequentist approach #2: Exact Sampling Distribution
Computing the confidence interval from the exact sampling distribution takes a bit more work.
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...
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
from scipy import optimize

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