# Interval estimates

#### Estimation framework, reminder

The framework we consider is the following. We have $N$ data points modelled as a vector $y \in \mathbb{R}^N$. We have a model for the data, that is the data is assumed to have a distribution $ p(y|\theta)$ for a certain parameter $\theta$ that we wish to estimate. The function $\theta \rightarrow p(y|\theta)$ is called the likelihood. 

For instance, we have $N$ data points, independent statistically, each data point is assumed to follow a Gaussian distribution of mean $\mu$ and variance $\sigma_k^2$, denoted by $Y_k \sim G(\mu, \sigma_k^2)$. Let us suppose that the $\sigma_k^2$ are known, so that the only unknown parameter is $\mu$ (it plays the role of $\theta$ in the definition of the likelihood). The likelihood of the data point $Y_k$ is then 
$$p(y_k|\mu)  = \frac{1}{\sqrt{2\pi} \sigma_k} \mathrm{e}^{-\frac{1}{2\sigma_k^2} (y_k - \mu)^2}$$
Since all the $Y_k$ are independent, the likelihood of the data $Y$ is the product of the likelihoods of the $Y_k$,
$$p(y|\mu)  = \prod\limits_{k=1}^N  p(y_k|\mu) =\prod\limits_{k=1}^N \frac{1}{\sqrt{2\pi} \sigma_k} \mathrm{e}^{-\frac{1}{2\sigma_k^2} (y_k - \mu)^2} = \frac{1}{\sqrt{2\pi}^N \prod\limits_{k=1}^N \sigma_k} \mathrm{e}^{-\frac{1}{2} \sum\limits_{k=1}^N \frac{(y_k - \mu)^2}{\sigma_k^2} } $$

Furthermore, we might assume that the parameter itself has a prior distrbution. That is, the probability of $\theta$ before seeing the data is $p(\theta)$. In the example above, we could assume that $p(\mu)$ is a Gaussian distribution of mean $0$ and variance $\sigma_\mu$ 

#### Point estimates, reminder

In the previous lesson we have studied the point estimates. In the point estimates view, we have an estimator $\hat{\theta}:y \rightarrow \hat{\theta}(y)$ that takes as argument the data $y$ and outputs a value wanted to be close to  $\theta$. The error bar is then given as the variance or square root mean squared error ($\sqrt{\mathrm{MSE}}$) of $\hat{\theta}$.

Some point estimates ignore the prior distributions, while some take it into account. The most common estimators that do not involve the prior are the maximum likelihood and least square estimates. When the Likelihood of the data is Gaussian and the covariance is known, they are equivalent. In the example above, the maximum likelihood estimate is 
$$\hat{\mu}_{ML} = \arg \max_\mu  p(y|\mu)  = \frac{\sum\limits_{k=1}^N \frac{y_k}{\sigma_k^2} }{\sum\limits_{k=1}^N \frac{1}{\sigma_k^2}} $$

If we assume a prior on $\mu$, $p(\mu)$, the common estimators are the mean, median and a posteriori, that are
$$ \hat{\theta}_{\mathrm{mean}} = \int_{-\infty}^\infty \mu p(\mu|y) \mathrm{d} \mu  =\int_{-\infty}^\infty \mu \frac{p(y|\mu) p(\mu)   }{p(y)}  \mathrm{d} \mu  $$
$$ \hat{\theta}_{\mathrm{median}} = \mathrm{median}(p(\mu|y)) $$
$$ \hat{\theta}_{\mathrm{mode}} = \mathrm{mode}(p(\mu|y)) $$
where the mode is the argument that maximizes the a function, $\mathrm{mode}(p(\mu|y)) = \arg \max_\mu p(\mu|y)$.

In the example above, $$\hat{\theta}_{\mathrm{mean}} = \hat{\theta}_{\mathrm{median}} = \hat{\theta}_{\mathrm{mode}} = \frac{\sum\limits_{k=1}^N \frac{y_k}{\sigma_k^2} }{\frac{1}{\sigma_\mu^2} +\sum\limits_{k=1}^N \frac{1}{\sigma_k^2}} $$.

If the model is correct, the posterior mean and median have respectively minimal mean squared error and mean absolute error. 

#### Interval estimates

In this spreadsheet, we change the viewpoint of the estimation. Instead of aiming at finding an estimator that is optimal in a certain sense, we consider the question: how likely is it that the true value of the parameters lie in a certain interval ? 

''Likely'' is a loose term that needs clarifications. There are two main ways of constructing interval estimates: the confidence intervals and the credible intervals, which have different properties.



## Confidence interval

A confidence interval is constructed in the following way. Given a likelihood $p(y|\theta)$ and data $y$, a confidence interval is constructed by choosing a probability $\alpha$, and two functions of the data $l_\alpha(y)$ and $u_\alpha(y)$ such that 
$$ \mathrm{Pr}\left\{  \theta \in [l_\alpha(y),   u_\alpha(y) ] \;  | \; \theta \right\} = \alpha $$

We first consider an example where we construct a confidence interval for the weighted mean of independent Gaussian variables. 
$$\hat{\mu}_{ML} = \arg \max_\mu  p(y|\mu)  = \frac{\sum\limits_{k=1}^N \frac{y_k}{\sigma_k^2} }{\sum\limits_{k=1}^N \frac{1}{\sigma_k^2}} $$

$\hat{\mu}_{ML}$ has a Gaussian distribution of variance $\sigma_{\hat{\mu}}^2 = \frac{1}{\sum\limits_{k=1}^N \frac{1}{\sigma_k^2}}$. 



```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import scipy.special as sp
```


```python
# Number of simulations
Nsim = 100000
N = 10 # Number of data points 
mean_error = 1 #Mean value of the error bars
mu = 1
#alpha = 4
errors_sigma = mean_error*np.random.chisquare(4,size=N)/4 #Generage values of the error bars
k = 3
conditions = np.ones(Nsim, dtype=bool)

for i in range(Nsim):
    
    y = mu + np.random.randn(N)*errors_sigma
    
    mu_estim = np.sum(y/errors_sigma**2) / np.sum(1/errors_sigma**2)
    sigma_estim = 1/ np.sqrt(np.sum(1/errors_sigma**2))
    
    u = mu_estim  + k*sigma_estim
    l = mu_estim  - k*sigma_estim
    condition = (mu <= u) * (mu >= l)
    
    conditions[i] = condition 
    
print('The true value of the parameter is in the interval in', np.sum(conditions) /Nsim*100, '% of the trials')


```

    The true value of the parameter is in the interval in 99.724 % of the trials


### Question 1

Compute analytically the probability that the true parameter is in the confidence interval if $l = \hat{\mu} - k \sigma_\hat{\mu}$ and $u = \hat{\mu} + k \sigma_\hat{\mu}$ for $k = 1, 2, 3$. 

Does the value of the confidence interval depend on the number of data points ? 

Compute the centered confidence interval that gives an inclusion probability  $ \alpha = 50\%$ 

Given that the distribution of data point is Gaussian with mean $\mu$ and error $\sigma$, we know that, $\alpha$ can be calculated using the following definition,

$$\alpha = \int_{\mu-k\sigma}^{\mu + k\sigma} \frac{1}{\sqrt{2\pi} \sigma} e^{-\frac{(y-\mu)^2}{2\sigma^2}} dy$$

Making a change of variables, $\left(\frac{y-\mu}{\sigma}\right)^2 = u^2$, we can write above equation as,

\begin{equation*}
    \begin{split}
        \alpha &= \frac{1}{\sqrt{2\pi}\sigma} \int_{-k}^{k} e^{-\frac{1}{2}u^2} \sigma du \\
        &= \sqrt{\frac{2}{\pi}} \int_{0}^{k} e^{-\frac{1}{2}u^2} du
    \end{split}
\end{equation*}

In the last line we used the fact that the integral is the even function around the limit. To solve this integration we can make another change of variable $\frac{u^2}{2} = x$,

\begin{equation*}
    \begin{split}
        \alpha &= \sqrt{\frac{2}{\pi}} \int_{0}^{k^2/2} e^{-x} \frac{1}{\sqrt{2x}} dx \\
        &= \frac{1}{\sqrt{\pi}} \int_0^{k^2/2} x^{-1/2} e^{-x} dx \\
        &= \frac{1}{\sqrt{\pi}} \left[-\sqrt{\pi} (1 - erf(\sqrt{x}) \right]^{k^2/2}_0
    \end{split}
\end{equation*}

In the last line we used the integral tables to find the value of given integration. Solving last equation one would get that,

$$\alpha = erf\left(\frac{k}{\sqrt{2}}\right)$$

We can check this formula using scipy as follows,


```python
k1 = np.array([1,2,3])
alpha = sp.erf(k1/np.sqrt(2))

print('For k=1, the probability that the true value would be in interval is ', alpha[0])
print('For k=2, the probability that the true value would be in interval is ', alpha[1])
print('For k=3, the probability that the true value would be in interval is ', alpha[2])
```

    For k=1, the probability that the true value would be in interval is  0.6826894921370859
    For k=2, the probability that the true value would be in interval is  0.9544997361036416
    For k=3, the probability that the true value would be in interval is  0.9973002039367398


This probability would not depend on the number of data points. (In the last calculation, we didn't use number of data points anywhere, we just used the PDF of the data).

Now, we want to calculate the centered confidence interval that gives the inclusion probability $\alpha=0.5$. That means we want to compute $k$ for given alpha which can be done using the inverse error function.

$$k = \sqrt{2} \cdot erf^{-1}(\alpha)$$

We can calculate this using the scipy.


```python
alpha1 = 0.5
kk = np.sqrt(2)*sp.erfinv(alpha1)
print('The confidence interval that gives the inclusion probability of 50% would be at about ' 
      + str(np.around(kk,2)) + 
      '-sigma from the center')
```

    The confidence interval that gives the inclusion probability of 50% would be at about 0.67-sigma from the center


### Question 2

We now consider another example. Suppose we observe
$$Y = \theta + \epsilon$$ where $\epsilon$ follows an exponential distribution
$f(\epsilon) = \frac{1}{\lambda} \exp(-\frac{\epsilon}{\lambda})$ and $\theta$ is the parameter to estimate.

Given the data $y$, construct confidence intervals for  68.27, 95.45 and 99.73 $\%$ for $\theta$ of the form $[y - x_\alpha ,y]$. In other words, find $x_\alpha$ such that $\theta \in [y - x_\alpha ,y]$ with a probability $\alpha$.

Check your calculations with a simulation as above. 



The likelihood function for the given exponential distribution would be,

$$p(y|\theta) = \frac{1}{\lambda}\exp{\left(-\sum_k \frac{y_k}{\lambda}\right)}$$

We can calculate the Maximum Likelihood estimate of $\hat{\lambda}$ as follows,

\begin{equation}
    \begin{split}
        \log p(y|\theta) &= - \log \lambda - \sum_k \frac{y_k}{\lambda} \\
        \Rightarrow \frac{d \log p(y|\theta)}{d\lambda} &= -\frac{1}{\lambda} + \sum_k \frac{y_k}{\lambda^2} \\
        \Rightarrow 0 &= -1 + \sum_k \frac{y_k}{\lambda} \\
        \Rightarrow \hat{\lambda} &= \sum_k y_k
    \end{split}
\end{equation}

Now, we want to find confidence interval for this distribution. We can do so as we did in the previous case. Let, $\alpha$ be probability with which true value of $\lambda$ lies in the given interval $(0,k)$. Then,

\begin{equation}
    \begin{split}
        \alpha &= \int_0^k \frac{1}{\lambda} e^{-x/\lambda} dx \\
        &= \frac{1}{\lambda} \left( \frac{e^{-x/\lambda}}{-1/\lambda} \right)_0^k \\
        &= 1 - e^{-k/\lambda}
    \end{split}
\end{equation}

Here, $\lambda$ would be the ML estimate of the parameter. Using above formula, we can find the interval $(0,k)$ for which the $\hat{\lambda}$ would be in the interval would be $\alpha$.

$$k = \lambda \ln{(1-\alpha)}$$


```python
# Number of simulations
Nsim = 100000
N = 10 # Number of data points 

mu = 2
alpha = 0.95

conditions1 = np.ones(Nsim, dtype=bool)

for i in range(Nsim):
    
    y = np.random.exponential(mu,N)
    
    mu_estim = np.sum(y)
    
    l = 0
    u = mu_estim*np.log(1-alpha)
    condition = (mu <= u) * (mu >= l)
    
    conditions1[i] = condition 
    
print('The true value of the parameter is in the interval in', np.sum(conditions1) /Nsim*100, '% of the trials')
```

    The true value of the parameter is in the interval in 0.0 % of the trials


### Question 3

We now consider another example. Suppose we observe
$$Y = \theta + \epsilon$$ where $\epsilon$ follows a gamma distribution of parameters $\alpha, \beta$.

Given the data $y$, construct confidence intervals for  68.27, 95.45 and 99.73 $\%$ for $\theta$ of the form $[y - x_\alpha ,y]$. In other words, find $x_\alpha$ such that $\mu \in [y - x_\alpha ,y]$ with a probability $\alpha$.

Check your calculations with a simulation as above. 


```python

```

4c63b71f0ad476355c929c9315296e76914998a0
ipynb
Jupyter Notebook
Lec10_Confidence_Intervals.ipynb
Jayshil/Astro-data_science
8f83643197cf05981e09490352caeec3f0cde4ae
MIT
Lec10_Confidence_Intervals.ipynb
Jayshil/Astro-data_science
8f83643197cf05981e09490352caeec3f0cde4ae
MIT
Lec10_Confidence_Intervals.ipynb
Jayshil/Astro-data_science
8f83643197cf05981e09490352caeec3f0cde4ae
MIT
Qwen/Qwen-72B
1. YES
2. YES