\documentclass[a4paper]{article}

\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{fullpage}
\usepackage{mathtools}
\usepackage{float}
\usepackage{enumerate}
\usepackage{hyperref}
\usepackage{comment}

\setlength{\parindent}{0pt}

\usepackage{amsthm}
\newtheorem{theorem}{Theorem}[section]
\newtheorem*{corollary}{Corollary}
\newtheorem{lemma}{Lemma}[section]
\newtheorem{claim}{Claim}[section]
\newtheorem{axiom}{Axiom}[section]

\theoremstyle{definition}
\newtheorem{definition}{Definition}[section]

\theoremstyle{remark}
\newtheorem{remark}{Remark}[section]
\newtheorem{example}{Example}[section]
\newtheorem{observation}{Observation}[section]

\let\originalleft\left
\let\originalright\right
\renewcommand{\left}{\mathopen{}\mathclose\bgroup\originalleft}
\renewcommand{\right}{\aftergroup\egroup\originalright}
\newcommand{\laplace}[1]{\mathcal{L}\left\{#1\right\}}
\newcommand{\pset}[1]{\mathcal{P}\left(#1\right)}
\newcommand{\pmat}[1]{\begin{pmatrix}#1\end{pmatrix}}
\newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\C}{\mathbb{C}}
\renewcommand{\P}{\mathbb{P}}
\newcommand{\inprod}[1]{\left\langle#1\right\rangle}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\val}{val}
\DeclareMathOperator{\ord}{ord}
\DeclareMathOperator{\charac}{char}
\DeclareMathOperator{\im}{im}
\DeclareMathOperator{\id}{id}
\DeclareMathOperator{\intemb}{int}
\DeclareMathOperator{\com}{com}
\DeclareMathOperator{\vecspan}{span}
\DeclareMathOperator{\defequiv}{\mkern4mu :\mkern-8mu \iff}

\newcommand{\ind}[1]{\mathbf{1}{\left\{#1\right\}}}
\newcommand{\norm}[1]{\left\lVert#1\right\rVert}
\newcommand{\hateq}{\mathrel{\widehat{=}}}

\makeatletter
\newcommand{\axcirc}[1]{\mathbin{\mathpalette\make@circled#1}}
\newcommand{\make@circled}[2]{%
	\ooalign{$\m@th#1\largecirc{#1}$\cr\hidewidth$\m@th#1#2$\hidewidth\cr}%
}
\newcommand{\largecirc}[1]{%
	\vcenter{\hbox{\scalebox{1.5}{$\m@th#1\bigcirc$}}}%
}
\makeatother

\newcommand{\letprob}{Let $(S, \Sigma, P)$ be a probability space. {}}
\let\letprobalways\letprob

% (un)comment to remove useless fluff
\renewcommand{\letprob}{}

\let\oldhat\hat
\let\oldtilde\tilde
\let\tilde\widetilde
\let\hat\widehat

\title{Mathematical Statistics 2}
\author{Koen Oostveen}

\begin{document}
	\maketitle
	\tableofcontents
	\section{Preliminaries}
	Let $V$ be a finite-dimensional $\R$-vector space that is also integrable. A \textbf{random variable} $X$ on $V$ and a probability space $(S, \Sigma, P)$ is a map $X : S \to V$. If $X$ has a density $f : V \to \R$, we define the \textbf{expectation} of $X$ through
	\[ EX := \int_V x f(x) \, dx \]
	which is well-defined due to all densities being equal up to $L^2$ quotient.
	\begin{lemma}
		We have that (in any basis) $(EX)^i = E(X^i)$, where $X^i : S \to \R$ is the $i$-th component of $X$ (pointwisely).
	\end{lemma}
	\begin{lemma}
		For any two random variables $X, Y$ on the same measure and vector spaces we have $E(X + Y) = EX + EY$. That is, $E$ is a linear operator from the space of random variables to $V$. The space of random variables is also a vector space (actually a module), where the field (ring) is the random variables $S \to \R$.
	\end{lemma}
	\begin{proof}
		Go into a basis and establish the equality componentwise through elementary results of probability theory.
	\end{proof}
	\begin{lemma}
		Let $A : V \to W$ be a linear map on vector spaces $V, W$. Then it holds that for any random variable $X$ we have $E[AX] = AE[X]$.
	\end{lemma}
	\begin{proof}
		Trivial due to linearity. We can pick a basis $e_1, \ldots, e_n$ for $V$ for example to have
		\[ E[AX] = A(e_i) E[X^i] = A(e_i) (EX)^i = A(EX) \]
	\end{proof}
	\begin{definition}
		The \textbf{covariance operator} $\operatorname{Cov}(X)$ for a random variable $X$ is given in a basis through
		\[ \operatorname{Cov}(X)^i_j := \operatorname{Cov}(X^i, X^j) \]
	\end{definition}
	\begin{lemma}
		On an inner product space with an orthonormal basis, we also have that
		\[ \operatorname{Cov}(X) = E[\inprod{X - EX, X - EX}] \]
		where the inner product is understood componentwise I suppose.
	\end{lemma}
	\begin{lemma}
		Let $A : V \to W$ be a linear map and let $X$ be a random variable on $V$. We have that
		\[ \operatorname{Cov}(AX) = A\operatorname{Cov}(X)A^* \]
		where $A^*$ is the Hermitian adjoint. Remember that this is the unique linear map such that for any two vectors $x, y \in V$ we have $\inprod{A(x), y} = \inprod{x, A^*(y)}$.
	\end{lemma}
	\begin{proof}
		We pick an orthonormal basis and we work it out. Remember that on an orthonormal basis we have the conjugate transpose blabla. Observe that
		\begin{align*}
			\operatorname{Cov}(AX)^i_j &= \operatorname{Cov}(AX^i, AX^j) \\
			&= \operatorname{Cov}(A^i_k X^k, A^j_k X^k) \\
			&= A^i_k A^j_k \operatorname{Cov}(X^k, X^k) \\
			&= A^i_l (A^*)^k_j \operatorname{Cov}(X^l, X^k) \\
			&= A^i_l \operatorname{Cov}(X)^l_k (A^*)^k_j \\
			&= (A\operatorname{Cov}(X))^i_k (A^*)^k_j \\
			&= (A\operatorname{Cov}(X)A^*)^i_j
		\end{align*}
	\end{proof}
	\begin{definition}
		A linear map $A : V \to W$ on inner product spaces is said to be \textbf{positive semi-definite} is for all $v \in V$ we have $\inprod{v, A(v)} \geq 0$.
	\end{definition}
	\begin{lemma}
		The covariance operator is positive semi-definite.
	\end{lemma}
	\begin{definition}
		A random variable $X$ is said to be \textbf{multivariate normally distributed} with mean $\mu \in V$ and invertible covariance matrix $\Sigma$ if	its density $f_X : V \to \R$ satisfies
		\[ f_X(x) = \frac{1}{\sqrt{\det(2 \pi \Sigma)}} \exp\left( -\frac{1}2 \inprod{x - \mu, \Sigma^{-1}(x - \mu)} \right) \]
		and we write $X \sim N(\mu, \Sigma)$.
	\end{definition}
	\begin{lemma}
		If $X \sim N(\mu, \Sigma)$ then $EX = \mu$ and $\operatorname{Var}(X) = \Sigma$.
	\end{lemma}
	\begin{remark}
		You don't necessarily have that if $X$ and $Y$ are normal, then $(X, Y)$ is multivariate normal. Example, let $X \sim N(0, 1)$. Let $\varepsilon$ be independent noise, any r.v, i.e. $\operatorname{Cov}(X, \varepsilon) = 0$, such that also $P(\varepsilon = 1) = P(\varepsilon = -1) = 1/2$. Let $Y = \varepsilon X$. We can show that $X$ and $Y$ are distributionally equal, i.e. $P(Y \leq y) = P(X \leq y)$. Then the support of the standard normal is not all of $V$, namely, only the subset where $Y = X$ or $Y = -X$.
	\end{remark}
	\begin{lemma}
		If $X$ and $Y$ are normally distributed and independent, then together they are bivariate normal.
	\end{lemma}
	\begin{proof}
		The joint density can be computed. Let $X \sim N(\mu, \sigma^2)$ and $Y \sim N(\tilde \mu, \tilde \sigma^2)$. Let $\mu' := (\mu, \tilde \mu)$ and
		\[ \Sigma := \begin{bmatrix}
			\sigma^2 & 0\\
			0 & \tilde \sigma^2
		\end{bmatrix}, \qquad \Sigma^{-1} = \begin{bmatrix}
		1/\sigma^2 & 0\\
		0 & 1/\tilde \sigma^2
		\end{bmatrix} \]
		Then we have $(X, Y) \sim N(\mu', \Sigma)$. We can see this by
		\begin{align*}
			f_{X, Y}(x, y) &= f_X(x) f_Y(y) \\
			&= \frac{1}{\sqrt{2 \pi \sigma^2} \sqrt{2 \pi \tilde \sigma^2}} \exp\left(- \frac{1}{2} \left( \left(\frac{x - \mu}{\sigma}\right)^2 + \left(\frac{y - \tilde \mu}{\tilde \sigma}\right)^2 \right) \right) \\
			&= \frac{1}{\sqrt{\det(2 \pi \Sigma)}} \exp\left( -\frac{1}2 \inprod{(x, y) - \mu', \Sigma^{-1}((x, y) - \mu')} \right) \qedhere
		\end{align*}
	\end{proof}
	\begin{corollary}
		Let $X_1, \ldots, X_n \sim N(\mu, \sigma^2)$ be i.i.d. Then $X := (X_1, \ldots, X_n)$ is multivariate normal with $\mu' := (\mu, \ldots, \mu)$ and $\Sigma := \sigma^2 \id$ so that $X \sim N(\mu', \Sigma)$.
	\end{corollary}
	\begin{lemma}
		Let $V$ and $W$ be appropriate vector spaces and let $X$ be a random variable on $V$. Let $A : V \to W$ be linear and invertible. Let $Y := AX$. We have that
		\[ f_Y(y) = f_X(A^{-1}(y))|\det(A^{-1})| \]
	\end{lemma}
	\begin{proof}
		The density is the unique function that satisfies that for all measurable $U \subseteq W$ we have
		\[ P(Y \in U) = \int_U f_Y(y) \, dy \]
		So we compute the above:
		\begin{align*}
			P(Y \in U) &= P(AX \in U) \\
			&= P(X \in A^{-1}(U)) \\
			&= \int_{A^{-1}(U)} f_X(x) \, dx \\
			&= \int_U f_x(A^{-1}(y)) |\det(A^{-1})| \, dy \qedhere
		\end{align*}
	\end{proof}
	\begin{lemma}
		Let $V$ and $W$ be appropriate vector spaces and let $X \sim N(\mu, \Sigma)$ on $V$. Let $A : V \to W$ be linear and invertible. We have that
		\[ AX \sim N(A\mu, A \circ \Sigma \circ A^*) \]
		where $A^*$ is the Hermitian adjoint. Recall that this is the unique map that satisfies
		\[ \inprod{A(x), y} = \inprod{x, A^{*}(y)} \]
	\end{lemma}
	\begin{proof}
		We go into components in an orthonormal basis. We then compute the density. Let $Y := AX$.
		\begin{align*}
			f_{Y}(y) &= f_X(A^{-1}(y)) |\det(A^{-1})| \\
			&= \frac{1}{\sqrt{\det(2 \pi \Sigma)}} \exp\left( -\frac{1}2 \inprod{A^{-1}(y) - \mu, \Sigma^{-1}(A^{-1}(y) - \mu)} \right) |\det(A^{-1})| \\
			&= \frac{1}{\sqrt{\det(2 \pi \Sigma)}} \exp\left( -\frac{1}2 \inprod{A^{-1}(y - A\mu), (A^* \circ (A^*)^{-1} \circ \Sigma^{-1} \circ A^{-1})(y - A\mu)} \right) |\det(A)|^{-1} \\
			&= \frac{1}{\sqrt{\det(2 \pi A^2 \Sigma)}} \exp\left( -\frac{1}2 \inprod{y - A\mu, (A \circ \Sigma \circ A^*)^{-1}(y - A\mu)} \right) \\
			&= \frac{1}{\sqrt{\det(2 \pi A^2 \Sigma)}} \exp\left( -\frac{1}2 \inprod{y - A\mu, (A \circ \Sigma \circ A^*)^{-1}(y - A\mu)} \right)
		\end{align*}
		We are done if we can show that $\det(A^2 \circ \Sigma) = \det(A \circ \Sigma \circ A^*)$. Well, we have that $\det(A^*) = \det(A)$ due to previous linear algebra results, and because everything is invertible this is clear.
	\end{proof}
	\begin{lemma}
		Let $X_1, \ldots, X_n \sim N(\mu, \sigma^2)$ be independent. Then
		\[ \frac{(n-1) \hat \sigma^2 }{\sigma^2} = \frac{1}{\sigma^2} \sum_{i=1}^n (X_i - \bar X)^2 \sim \chi^2_{n-1} \]
		where $\bar X := \frac{1}{n} \sum_{i=1}^n X_i$ and $\hat \sigma^2 := \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar X)^2$.
	\end{lemma}
	\begin{proof}{\textbf{TODO}}
		Let $X := (X_1, \ldots, X_n)$.
		\[ (n-1)\hat \sigma^2 = \inprod{X - \bar X, X - \bar X} \qedhere \]
	\end{proof}
	\begin{theorem}[Multivariate CLT]
		Let $(X_k)_{k \in \N^*}$ be a sequence of random variables in $\R^d$, independent and identically distributed (iid). Suppose for all $k \in \N^*$ we have $E X_k = \mu$ for some $\mu \in \R$ and $\operatorname{var}(X_k) = \sigma^2$ for $\sigma > 0$. Then
		\[ \sqrt{n}(\overline{X}_n - \mu) \to N(0_{\R^d}, \Sigma) \]
	\end{theorem}
	\newpage
	\section{Regression}
	Let $\beta_0 \in \R$ and $\beta \in \R^n$, $\varepsilon_1, \ldots, \varepsilon_m$ be independent noise with mean $0$, let additionally for $i=1, \ldots, m$ we have $X_i$ be random variables in $\R^n$. We say that $Y := (Y_1, \ldots, Y_m)$ follows a \textbf{linear model} if
	\[ Y_i = \beta_0 + \beta_j X_i^j + \varepsilon_i \]
	Let
	\[ X := \begin{bmatrix}
		1 & X_1^1 & \cdots & X_1^n \\
		\vdots & \vdots & \ddots & \vdots \\
		1 & X_m^1 & \cdots & X_m^n \\
	\end{bmatrix}, \beta := \bmat{\beta_0 \\ \vdots \\ \beta_n} \]
	so that we also have $Y = X\beta + \varepsilon$.
	\begin{lemma}
		If $\varepsilon \sim N(0, \sigma^2 \id)$ then $Y \sim N(X \beta, \sigma^2 \id)$ (where $X_i$ are not assumed to be random variables but concrete data).
	\end{lemma}
	\begin{lemma}
		Let $A : \R^n \to \R^m$ be linear (where we equip $\R^d$ with the standard inner product). Let $b \in \R^m$. If $A^T A$ is invertible, then
		\[ \arg \min_{x \in \R^n} \norm{Ax - b}_2^2 = (A^T A)^{-1}(A^T b) \]
	\end{lemma}
	\begin{proof}
		Let $y := (A^T A)^{-1}(A^T b)$. Let $r := b - Ay$. We claim that for all $p \in \im A$ we have $\inprod{r, p} = 0$. Let $\tilde p$ be a preimage of $p$, then
		\[ \inprod{r, p} = \inprod{b - A(A^T A)^{-1}(A^T b), A\tilde p} = \inprod{A^T b - A^T b, \tilde p} = 0 \]
		Consider that
		\begin{align*}
			\norm{Ax - b}^2_2 &= \norm{A(x - y) - (b - Ay)} \\
			&= \norm{A(x - y)}^2_2 + \norm{b - Ay}^2_2 \\
			&\geq \norm{A(x - y)}^2_2
		\end{align*}
		where we have equality if and only if $A(x - y) = 0 \implies x = y$, so $y$ is the unique minimizer.
	\end{proof}
	\begin{theorem}
		If $X^* \circ X$ is invertible, then the MLE and the least-squares estimator $\hat \beta$ for $\beta$ coincide and
		\[ \hat \beta = (X^* \circ X)^{-1}(X^*(Y)) \]
		with
		\[ \hat \beta \sim N(\beta, \sigma^2 (X^* \circ X)^{-1}) \]
	\end{theorem}
	\begin{proof}
		The likelihood of $Y$ is
		\[ f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left( -\frac{1}2 \left(\frac{\norm{x - X \beta}}{\sigma}\right)^2 \right) \]
		This is maximized by
		\[ \hat \beta := \operatorname*{argmin}_{\beta} \norm{Y - X \beta} \]
		which is exactly the least-squares estimator given by solving the normal equation. It is unbiased due to
		\[ E\hat \beta = \beta + E(X^* \circ X)^{-1}(X^*(\varepsilon)) = \beta \]
		and normality and the other properties can be derived from the earlier lemma's.
	\end{proof}
	\begin{remark}
		If $X^* X$ is not invertible, we can interpret that as the minimum being attained by infinitely many $\beta$, or alternatively that there are simply multiple linear models possible given the data. Such a model is called \textbf{non-identifiable}.
	\end{remark}
	\newpage
	\section{Logistic regression}
	Behold, a majestic smooth interpolator and actually bijection between $(0, 1)$ and $\R$: $f : \R \to (0, 1)$ defined by $x \mapsto f(x) := \left(1 + \exp(-x)\right)^{-1}$. So if we want to do a regression on a categorical variable $X \in \{A, B\}$, we can simply model the probability of category $A$, by taking a linear model (something we understand well) and stuffing it into the logistic function, to get a probability out. I.e. we use the model
	\[ P(X = A) = \left(1 + \exp\left(-\left(\beta_0 + \beta_1 x_1 + \ldots + \beta_n x_n\right)\right)\right)^{-1} \]
	Now let's invert the logistic function, which is by the way invertible as we saw. For $x \in \R$ and $y \in (0, 1)$ we have
	\[ y = (1 + \exp(-x))^{-1} \iff x = -\log(y^{-1} - 1) = \log\left(\frac{y}{1-y}\right) =: \operatorname{logit}(y) \]
	i.e. we define the \textbf{logit} function to be the inverse of the logistic function. We are left with
	\[ \operatorname{logit}(P(X = A)) = \beta_0 + \beta_1 x_1 + \ldots + \beta_n x_n \]
	If we want to now make an estimator, we got to have a statistical model. We still observe data like before and assume that the outputs $Y_i$ take either $0$ or $1$ with probability $p_i$, i.e. Bernoulli distributed.
	Let us compute the likelihood
	\[ L = \prod_{i=1}^n (1-p_i)^{1 - Y_i} p_i^{Y_i} \]
	so that the log-likelihood is
	\[ \log L = \sum_{i=1}^n (1 - Y_i) \log(1-p_i) + Y_i \log(p_i) \]
	where
	\[ p_i = 1/(1 + \exp\left(-\left(\beta_0 + \beta_1 x_{i,1} + \ldots + \beta_n x_{i,n})\right)\right) \]
	well, good luck maximizing that\textellipsis
	
	\subsection{Extra: neural networks}
	Compare: a linear model can be seen as just a vector $\beta \in \R^n$ and intercept $\beta_0 \in \R$ such that the model can be defined by a function $f : \R^n \to \R$ such that $x \mapsto \inprod{x, \beta} + \beta_0$, whereas a (shallow) neural network is defined as a similar function such that $\sigma : \R^m \to \R^k$ is a (continuous) activation function, $W : \R^n \to \R^m$ a linear map, and $a \in \R^k$ some sort of coefficients, then the model $f$ is defined by
	\[ f(x) := \inprod{a, \sigma(Wx)} \]
	and of course if you choose $\sigma$, $W$ and $a$ appropriately, you recover the case of linear regression.
	Of course in general you cannot say anything about $\sigma$ apart from maybe convexity or continuity conditions, but in generality you cannot guarantee much, so you get a nonconvex optimization problem out, you require some sort of gradient descent algorithm (you assume the whole thing is partially differentiable on some basis).
	\newpage
	\section{The $\chi^2$ goodness-of-fit test}
	In lower generality, we have considered before tests of proportions, namely, a model with a parameter $p \in [0, 1]$, and we consider for some $p_0 \in [0, 1]$ the hypothesis testing problem
	\[ H_0 : p = p_0 \qquad \text{vs.} \quad H_1 : p \neq p_0 \]
	
	But what if we want to consider multiple proportions, namely, what if we measure proportions $p_1, \ldots, p_n \in [0, 1]$ and we want to test whether at least one $p_i \neq p_j$ for $i \neq j$. The setup is as follows.\\
	
	We measure categorical random variables $X_1, \ldots, X_n \in \{1, \ldots, k\}$, we have $k \in \N$ categories. Let $p_1, \ldots, p_k \in [0, 1]$ be the proportions of occurrences of categories $1, \ldots, k$ respectively, which we measure of course by (for $i = 1, \ldots, k$)
	\[ p_i := \frac{1}{n} \sum_{j=1}^n \ind{X_j = i} \]
	If $p_1, \ldots, p_n$ are the true proportions ($H_0$) then obviously
	\[ p_1 + \ldots + p_n = 1 \]
	or in other words $P(X_i = j) = p_j$ from which the above obviously follows. We first define the random variables
	\[ N_j := \sum_{j=1}^n \ind{X_j = i} =: \#\{i : X_i = j\} \]
	and then the random variable $N := (N_1, \ldots, N_k) \in \R^k$ which then by definition follows a \textbf{multinomial distribution}, we write $N \sim MN(n, (p_1, \ldots, p_k))$, where $n$ is obviously the amount of trials. Using combinatorical arguments we can show that 
	\[ P(N = (n_1, \ldots, n_k)) = \binom{n}{n1, \ldots, n_k} p_1^{n_1} \cdot \ldots \cdot p_k^{n_k} \ind{n_1 + \ldots +  n_k = 1} \]
	where
	\[ \binom{n}{n_1, \ldots, n_k} := \frac{n!}{n_1! \cdot \ldots \cdot n_k!} \]
	We can see that this is consistent with the binomial coefficient, since for $N \sim B(n, p)$ we have
	\[ P(N = k) = \binom{n}{k} p^k (1-p)^{n-k} = \binom{n}{k, n-k} p_1^{n_1} p_2^{n_2} \]
	where $n_1 := k$, $n_2 := n - k$, $p_1 := p$, $p_2 := 1 - p$.
	We can show that additionally for $N \sim MN(n, (n_1, \ldots, n_k))$
	\[ EN = (np_1, \ldots, np_k) \tag{easily consistent with binomial dist} \]
	\[ \operatorname{Cov}(N)^i_j = \operatorname{Cov}(N_i, N_j) = \begin{cases*}
		np_i (1 - p_i) & if $i = j$ \\ 
		-np_ip_j & if $i \neq j$
	\end{cases*} \]
	which is once again consistent with the binomial distribution, so we have found a reasonable generalization. And just like the binomial distribution, we have asymptotic normality, namely for `large $n$' we have
	\[ N \sim N(EN, \operatorname{Cov}(N)) \]
	With rescaling we can find a standard multivariate normal, namely
	\[ Z := \left(\frac{N_1 - np_1}{\sqrt{np_1}}, \ldots, \frac{N_k - np_k}{\sqrt{np_k}}\right) \]
	with $EZ = 0$ and
	\[ \operatorname{Cov}(Z) = \id(1 - p_1 - \ldots - p_k) \]
	and from that it (apparently) follows that (approximately)
	\[ \sum_{j=1}^k \frac{(N_j - np_j)^2}{np_j} \sim \chi^2_{k-1} \]
	which we can compute quantiles on and that way do hypothesis testing. Namely, we test
	\[ H_0 : P(X_i = j) = p_j, \quad i = 1, \ldots, k \]
	where we reject if
	\[ \underbrace{\sum_{j=1}^k \frac{(N_j - np_j)^2}{np_j}}_{=: X} > c_\alpha \]
	with $\alpha$ the level and $P(X > c_\alpha) = \alpha$.\\
	\section{Contingency tables / The $\chi^2$ independence test}
	Let's figure out if two samples of categorical variables are related in some sense. We observe a sample of $n$ pairs of categorical random variables $(X_i, Y_i)$ for every $i = 1, \ldots, n$, where $X_i \in \{1, \ldots, s\}$ and $Y_i \in \{1, \ldots, t\}$, with proportions $p_1, \ldots, p_s$ and $q_1, \ldots, q_t$ respectively. We define the symbol
	\[ r_{k,l} := P(X_i = k, Y_i = l) \]
	Now if we have independence of the $X$ and $Y$ variables, we have
	\[ r_{k,l} = P(X_i = k, Y_i = l) = P(X_i = k) P(Y_i = l) = p_k q_l \]
	for all $k, l$, which is then our null hypothesis $H_0$, whereas the alternative states that independence does not hold for at least one pair $k, l$. Let's find an estimator.
	\[ \hat r_{k, l} = \frac{1}{n} \sum_{i=1}^n \ind{x_i = k, y_i = l} \]
	and likewise we can find estimators $\hat p_k$ and $\hat q_l$. We can from there show that the statistic
	\[ T := n \sum_{k, l} \frac{(\hat r_{k, l} - \hat p_k \hat q_l)^2}{\hat p_k \hat q_l} \]
	is such that $T \sim \chi^2_{(s-1)(t-1)}$ approximately / asymptotically. Yay.
	\newpage
	
	\section{Nonparameteric tests}
	Hypothesis tests that do not involve parameters of distributions. Given i.i.d. $X_1, \ldots, X_n$ random variables. We pose the testing problem: $H_0 : X_1 \text{ generated from normal distribution}$ v.s. not of course. We do not know a priori anything about the distribution of the $X_i$'s. Advantage: we do not need many assumptions, disadvantage: it will be not very powerful.\\
	\subsection{Test for normality}
	Well, we can use the fact that $X_1, \ldots, X_n$ are normal under $H_0$. So, we can measure how much the data `deviates' from the normal distribution. We use the skewness and kurtosis for this.
	\begin{definition}
		Let $X$ be a random variable with finite 2nd and 3rd moment, mean $\mu$, and standard deviation $\sigma > 0$. We define the \textbf{skewness}, $\gamma_1(X)$, of $X$ by
		\[ \gamma_1(X) := E[(X-\mu)^3/\sigma^3] \]
		recall that $\sigma := \sqrt{\operatorname{Var}(X)}$.
	\end{definition}
	\begin{lemma}
		If the p.d.f. of $X$ is symmetric around $\mu$, then the skewness is zero.
	\end{lemma}
	\begin{proof}
		Let $f : \R \to [0, \infty)$ be the p.d.f. of $X$. Then
		\begin{align*}
			\gamma_1(X) &= E[(X-\mu)^3/\sigma^3] \\
			&= \int_\R \frac{(x-\mu)^3}{\sigma^3} f(x) \, \mathrm dx \\
			&= \frac{1}{\sigma^3} \int_\R (x-\mu)^3 f(x) \, \mathrm dx \\
			&= \frac{1}{\sigma^3} \underbrace{\int_\R x^3 f(x + \mu) \, \mathrm dx}_0 = 0
		\end{align*}
		where the last step is justified, because the integrand $x \mapsto x^3 f(x + \mu)$ is odd,
		\[ (-x)^3 f(-x + \mu) = -x^3 f(x + \mu) \qedhere \]
	\end{proof}
	\begin{corollary}
		The skewness of a normally distributed $X$ is $0$.
	\end{corollary}
	\begin{definition}
		Let $X$ be a random variable with finite 2nd and 4th moment, mean $\mu$, and standard deviation $\sigma > 0$. We define the \textbf{excess kurtosis}, $\gamma_2(X)$, of $X$ by
		\[ \gamma_2(X) := E[(X-\mu)^4/\sigma^4] - 3 \]
	\end{definition}
	Why the $-3$? This is exactly because the normal distribution that way also has excess kurtosis $0$.
	\subsubsection{Jarque-Bera test}
	We measure the sample skewness and kurtosis by replacing the expectations and moments by estimators. Let $\bar X$ be the sample mean, let
	\[ s_X := \sqrt{\frac{1}{n} \sum_{i=1}^n (X_i - \bar X)^2 } \]
	Then we use these estimators as follows:
	\[ S := \frac{1}{n} \sum_{i=1}^n (X_i - \bar X)^3/s_X^3, \qquad K' := \frac{1}{n} \sum_{i=1}^n (X_i - \bar X)^4/s_X^4 - 3 \]
	One can show that the statistic
	\[ JB := \frac{n}{6}\left(S^2 + K'^2/4 \right) \]
	has approximately / asymptotically a $\chi^2_2$ distribution. So we can do hypothesis testing! (for large samples). Though, the power of this test is quite low, especially for distributions that have similar 3rd and 4th moments like the normal distribution, so that they have skewness and excess kurtosis of $0$.\\
	
	Clearly, just using certain statistics, namely $\bar X$, $s_X$, skewness, kurtosis, etc. tell you information about the distribution, but there will always be an information loss. But if we use the order statistics, no information will be lost, you know \textbf{every} quantile of the measured distribution as good as you can get it with your sample. In fact, we can use the expectation of the order statistics given a certain sample size $n$ to make a \textbf{QQ-plot}. The follow test simulates this plot.
	
	\subsubsection{Shapiro-Francia test}
	Let $X_1, \ldots, X_n \sim N(0, 1)$ be an i.i.d. sample of the standard normal distribution. Let
	\[ X_{(1)} \leq \ldots \leq X_{(n)} \]
	be the corresponding order statistics. Let
	\[ X' := (X_{(1)}, \ldots, X_{(n)}), \qquad m := (EX_{(1)}, \ldots, EX_{(n)}) \]
	The \textbf{Shapiro-Francia test} uses the test statistic $W' := \operatorname{\hat {Corr}}$, i.e. $W'(X', m) := \operatorname{\hat {Corr}}(X', m)$. In practice this means we compute
	\[ W' = \frac{\sum_{i=1}^n (X_{(i)} - \bar X)(m_i - \bar m)}{\sqrt{\sum_{i=1}^n (X_{(i)} - \bar X)^2 \sum_{i=1}^n(m_i - \bar m)^2}} \]
	The test clearly works like this: we reject if $W' < c$ for appropriate critical value $c$.
	
	\subsubsection{Shapiro-Wilk test}
	Whereas the previous section essentially used $EX'$, this test also uses the `second moment', in the sense of the covariance matrix $\operatorname{Cov}(X') =: V$, which is assumed to be invertible / positive definite.
%	Recall that for linear $A : V \to W$ and all $v \in V$, $w \in W$, we have that for the \textbf{Hermitian adjoint} $A^*$ of $A$ that,
%	\[  \]
	Let $a := m^T V^{-1}/\norm{m^T V^{-1}}_2$.
	\begin{remark}
		Recall that for linear $A : \R^n \to \R^m$ where both spaces are equipped with the standard inner product and a vector $w \in \R^m$, that $w^T A := (A^T w)^T := (A^* w)^*$, where the transpose is essentially the index-pulling method that we get from choosing a basis, the standard basis. This vector satisfies that for $v \in \R^n$, $(w^T A) v = w^T (A v)$, i.e. it guarantees something about an inner product, namely that $\inprod{w^T A, v} = \inprod{w, Av} = \inprod{A^* w, v}$. Aha!
	\end{remark}
	Anyway, with such $a$, we get that $\norm{a}_2 = 1$, encoding `weighted averages' of covariances scaled by the mean. We can use this to compute the standard deviation, in the sense that
	\[ \inprod{a, X} = \inprod{m, V^{-1} X}/\norm{m^T V^{-1}}_2 = \sigma \]
	So we define the statistic $W := \frac{\inprod{a, X}}{\hat \sigma}$, where
	\[ \hat \sigma := \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar X)^2 \]
	Now we can show that $W \leq 1$, and $EW = 1$ under the null hypothesis, so we reject if $W < c$ for some appropriate quantile $c$, that is usually found through simulations.
	\begin{remark}~
		\begin{itemize}
			\item The test can encounter problems if many values are equal.
			\item The test rarely rejects $H_0$ when the sample size is small.
			\item  It rejects $H_0$ for tiny differences from normality when $n$ is large.
			\item A major problem is that normality is the null hypothesis.
			\item  In practice we are rather interested in the opposite problem.
			\item  For the construction it is crucial that data are normal under $H_0$ .
			\item  The limiting distribution (after standardization) is not $\chi^2$, but a weighted infinite
			sum of shifted squared standard normals. Critical regions are found by
			Monte-Carlo simulations.
		\end{itemize}
		In practice, the Shapiro-Wilk test is used together with other measures for normality,
		such as kurtosis, histograms and QQ-plots.
	\end{remark}
	
	\subsection{Permutation test (for distribution)}
	Let $X_1, \ldots, X_n$ and $Y_1, \ldots, Y_m$ be i.i.d. with means $\mu_X$, $\mu_Y$. We define the testing problem $H_0 : \mu_X = \mu_Y$ v.s. the negation. A simple idea for a test statistic is just $T := \bar X - \bar Y$ which has mean $0$, so we reject if $|T| > c$ for a critical value $c > 0$. But because we have not assumed anything about the distribution, this critical value (for the desired level) cannot be found. Well, if we can somehow approximate the c.d.f $P(T \leq t)$, we can find the quantiles.\\
	
	We can approximate this, however, by using the data that we already have. Additionally, under $H_0$, $X$ and $Y$ have identical distributions. More generally, we can have a testing problem with distributions $P, Q$ such that: $X_1, \ldots, X_n \sim P$, $Y_1, \ldots, Y_m \sim Q$, and we test
	\[ H_0 : P = Q \quad \text{v.s.}  \quad H_1 : P \neq Q \]
	 
	 \begin{definition}
	 	Let $n \in N^*$. A \textbf{permutation} of \textbf{arity} $n$ is a bijection $\phi : \{1, \ldots, n\} \to \{1, \ldots, n\}$.
	 \end{definition}
	 
	 The idea is that under $H_0$, permuting the data should not affect $T$ much if they are sampled from the same data, as $\bar X = \bar Y$ under $H_0$. So we can use \textbf{every} permutation (of which there exists finitely many, namely $(m+n)!$) of the vector $Z := (X_1, \ldots, X_n, Y_1, \ldots, Y_m)$, and find the $p$-value by testing how often a permutation goes past the mean in a certain direction (which under $H_0$ should be about half the time). We define $\hat p$, an estimator for the $p$-value, by first labelling $\phi_1, \ldots, \phi_N$ all permutations of $m + n$ elements, and
	 \[ \hat p(Z) := \frac{1}{N} \sum_{i=1}^n \ind{T(\phi_i(Z)) \geq T(Z)} \]
	 The test clearly becomes $\phi(Z) := \ind{\hat p(Z) < \alpha}$ at level $\alpha$. The Type I error, or empirical level, is
	 \[ P_0(\phi = 1) = P_0(\hat p(Z) < \alpha) \leq \alpha \]
	 How can we see the final step? Well, under $H_0$, we note that for any $i$,
	 \[ P(T(\phi_i(Z)) \geq T(Z)) = P(T(\phi_i(Z)) \leq T(Z)) = 1/2 \]
	 i.e.
	 \[ \sum_{i=1}^n \ind{T(\phi_i(Z)) \geq T(Z)} \sim B(N, 1/2) \]
	 
	 \subsection{Sign test (for median)}
	 \begin{definition}
	 	Let $X$ be a random variable. We say that $m \in \R$ is the \textbf{median} of $X$ if
	 	\[ P(X < m) \leq 1/2 \land P(X > m) \leq 1/2 \]
	 	If it is unique, we write $\operatorname{Med}(X) := m$. If $X$ has quantile function $Q_X$, $Q(1/2)$ is clearly a median.
	 \end{definition}
	 \begin{lemma}
	 	Let $X \sim N(\mu, \sigma^2)$, then $\operatorname{Med}(X) = \mu$.
	 \end{lemma}
	 \begin{proof}
	 	Consider that $P(X < \mu) = P(X > \mu) = 1/2$ so it is clearly a median. Now we show that is unique. Let $\varepsilon \neq 0$ be a deviation to the median, such that $\tilde m := m + \varepsilon$ is an alternative median. Then
	 	\begin{align*}
	 		P(X < \tilde m) &= P(m < X < \tilde m) + P(X < m) \\
	 		&=  \int_m^{\tilde m} f_X(x) \, dx + 1/2 \\
	 		&\geq (\tilde m - m) \sup_{x \in [m, \tilde m]} f_X(x) + 1/2 \\
	 		&> 1/2
	 	\end{align*}
	 	so $\tilde m$ is not a median, so by contradiction we are done.
	 \end{proof}
	 \begin{definition}[Sign test]
	 	Let $X_1, \ldots, X_n$ be iid that have a unique median, let $m_0 \in \R$. We test $H_0 : \operatorname{Med}(X_1) = m_0$ v.s. $H_1 : \operatorname{Med}(X_1) \neq m_0$. Define the test statistic
	 	\[ T := \sum_{i=1}^n \ind{X_i \leq m_0} \]
	 	Then we reject if $|T| > c$ for certain critical value $c$. This is called the \textbf{sign test} since
	 	\[ T = \sum_{i=1}^n \ind{\operatorname{sign}(X_i - m_0) = -1} \]
	 \end{definition}
	 \begin{lemma}
	 	In this sign test, $T \sim B(n, 1/2)$ under $H_0$.
	 \end{lemma}
	 \begin{proof}
	 	Clearly, $T$ is the sum of Bernoulli random variables $B_1, \ldots, B_n$, where under $H_0$ for each $i$
	 	\[ P(B_i = 1) = P(X_i \leq m_0) = 1/2 \]
	 	So then $T \sim B(n, 1/2)$.
	 \end{proof}
	 \begin{remark}
	 	For large $n$ one can use that $B(n, 1/2) \approx N(n/2, n/4)$.
	 \end{remark}
	 \subsection{Runs test (for distribution / stochastic ordering)}
	 \begin{definition}
	 	Let $X, Y$ be random variables. We say that $X \leq Y$, or \textbf{X is stochastically less than Y}, if for all $t \in \R$ we have
	 	\[ F_X(t) = P(X \leq t) \geq P(Y \leq t) = F_Y(t) \]
	 	i.e. the cdf of $X$ is bounded below by $Y$. This makes $\leq$ a partial ordering.
	 \end{definition}
	 \begin{theorem}
	 	Let $X$ and $Y$ be random variables with density $f_X$ and $f_Y$, respectively. tfae
	 	\begin{itemize}
	 		\item $X \leq Y$
	 		\item For all increasing functions $\phi : \R \to \R$, $E[\phi(X)] \leq E[\phi(Y)]$ (so in particular $EX \leq EY$)
	 		\item For all strictly increasing functions $\phi : \R \to \R$, $\phi(X) \leq \phi(Y)$
	 		\item For all $A \subseteq \R$ with increasing indicator, $P(X \in A) \leq P(Y \in A)$
	 	\end{itemize}
	 \end{theorem}
	 \begin{proof}
	 	Suppose $X \leq Y$. Let $\phi : \R \to \R$ be increasing. That is, for all $x, y \in \R$ such that $x \leq y$, we have $\phi(x) \leq \phi(y)$. Notice that for all nonnegative random variables $Z$
	 	\begin{align*}
	 		EZ &= \int_\R z f_Z(z) \, dz \\
	 		&= \int_\R \int_0^\infty \ind{z > t} f_Z(z) \, dt \, dz \\
	 		&= \int_0^\infty \int_\R \ind{z > t} f_Z(z) \, dz \, dt \\
	 		&= \int_0^\infty P(Z > t) \, dt
	 	\end{align*}
	 	Therefore, since $P(\phi(X) > t) \leq P(\phi(Y) > t)$, we have the desired result. Now suppose we have indeed $E[\phi(X)] \leq E[\phi(Y)]$ for all increasing functions. Clearly, if $\phi$ is strictly increasing it is in particular increasing. Let $t \in \R$, we want to show that $P(\phi(X) > t) \leq P(\phi(Y) > t)$. Let $g_n : \R \to \R$ be an increasing step function that converges to $g : \R \to \R$ defined by $x \mapsto \ind{u > t}$. Then by the assumption we have that $E[(g_n \circ \phi)(X)] \leq E[(g_n \circ \phi)(Y)]$. Then
	 	\[ P(\phi(X) > t) = \lim_{n \to \infty} E[(g_n \circ \phi)(X)] \leq \lim_{n \to \infty} E[(g_n \circ \phi)(Y)] = P(\phi(Y) > t) \]
	 	So we are done for this part. The next implications are trivial.
	 	%		We know that strictly increasing functions are invertible onto their image, and if $t$ is not in the image, the result is trivial so we can exclude that. So we need to instead just show that $P(X > \phi^{-1}(t)) \leq P(Y > \phi^{-1}(t))$. We know that
	 	%		\[ E[\phi(X)] = \int_0^\infty P(\phi(X) > p) \, dp = \int_0^t P(\phi(X) > p) \, dp + \int_t^\infty P(\phi(X) > p) \, dp \]
	 	%		\[ E[X] = \int_{-\infty}^\infty x f_X(x) \, dx = \int_0^\infty \int_t^\infty f_X(x) \, dx \, dt = \int_0^\infty P(X > t) \, dt  \]
	 	%		\[ E[\phi(X)] = \int_\R \phi(x) f_X(x) \, dx \leq \int_\R \phi(x) f_Y(x) \, dx = E[\phi(Y)] \]
	 \end{proof}
	 \begin{definition}[Rank of sample]
	 	Let $\pi : \{1, \ldots, n\} \to \{1, \ldots, n\}$ be a permutation (bijection) that orders a sample $z_1, \ldots, z_n$ into $z_{\phi(1)}, \ldots, z_{\phi(n)}$, so that the latter is sorted w.r.t $\leq$. We say that the \textbf{rank} of $z_i$ for all $i$ is $\phi(i)$. 
	 \end{definition}
	 \begin{definition}[Runs test]
	 	Let $X_1, \ldots, X_n \in \R$, $Y_1, \ldots, Y_m \in \R$ both from an iid sample (suppose $X \sim F$ and $Y \sim G$). We test $H_0 : F = G$ versus $H_1 : F \neq G$ by rejecting $H_0$ if $R < c$ for appropriate critical value $c$, where $R$ is the \textbf{amount of runs}, given by the following. We take the sample $X_1, \ldots, X_n, Y_1, \ldots, Y_m$ and find the ranks. Then we color the sequence $\pi(1), \ldots, \pi(n + m)$ by consecutive sequences, and count how many `runs' we have, this is how we define $R$.
	 \end{definition}
	 \begin{theorem}
	 	Under $H_0$, we have that
	 	\[ P(R = r) = huge monster help \]
	 \end{theorem}
	 \begin{remark}
	 	$R$ is asymptotically normal! So for large sample sizes, we can use that.
	 \end{remark}
	 \begin{remark}
	 	$R$ is a measure of stochastic ordering, so an alternative reformulation of \[H_1 : \text{distributions are stochastically ordered}\]
	 \end{remark}
	 
	 \subsection{Wilcoxon rank sum test (for distribution)}
	 \begin{definition}[Rank sum test]
	 	Let $X_1, \ldots, X_n \sim F$ and $Y_1, \ldots, Y_m \sim G$ be independent with continuous density, we test similarly again $H_0 : F = G$ v.s. $H_1 : F \neq G$. Let
	 	\[ W := \sum_{i=1}^n R(X_i) \]
	 \end{definition}
	 \begin{lemma}
	 	Clearly using well-known sum identities and bounding $R(\cdot)$
	 	\[ n(n+1)/2 \leq W \leq (n+m)(n+m+1)/2 - n(n+1)/2 \]
	 \end{lemma}
	 Recall that there are $(n+m)!$ permutations of arity $n + m$. Every outcome of
	 \[ (R(X_1), \ldots, R(X_n), R(Y_1), \ldots, R(Y_m)) \]
	 is equally likely, so $EW = n(n+m+1)/2$. Let $R_i := R(X_i)$. We can show this by observing that $R_i$ is uniform among $\{1, \ldots, n+m\}$, so $E[R_i] = \frac{n+m+1}2$. Therefore $EW=n(n+m+1)/2$. We can also show that
	 \[ \operatorname{Var}(R_i) = \frac{(n+m)^2 - 1}{12}, \qquad \operatorname{Cov}(R_i, R_j) = -\frac{n+m + 1}{12} \]
	 \begin{align*}
	 	\operatorname{Var}(W) &= \operatorname{Var}\left(\sum_{i=1}^n R_i\right) \\
	 	&= \sum_{i=1}^n \sum_{j=1}^n \operatorname{Cov}(R_i, R_j) \\
	 	&= \sum_{i=1}^n \frac{(n+m)^2 - 1}{12} - (n-1) \frac{n+m + 1}{12} \\
	 	&= \sum_{i=1}^n \frac{n^2 + 2nm + m^2 - 1}{12} - \frac{n^2 + nm + n - n - m - 1}{12} \\
	 	&= \sum_{i=1}^n m\frac{n + m + 1}{12} = nm\frac{n + m + 1}{12}
	 \end{align*}
	 Using this data we can use asymptotic normality for large samples and the CLT.
	 
	 \begin{theorem}
	 	Under $H_0$, $W$ has an approximate distribution of
	 	\[ N\left(\frac{n(n+m+1)}{2}, \frac{nm(n+m+1)}{12}\right) \]
	 \end{theorem}
	 
	 \subsection{Wilcoxon signed rank sum test (for median)}
	 Suppose that $X$ is a continuous random variable \textbf{with symmetric pdf around the mean}. We want to test $H_0 : m = m_0$ for some candidate median $m_0 \in \R$ versus $H_1 : m \neq m_0$ or similar. We define the \textbf{signed rank sum statistic} $W$ by
	 \[ W := \sum_{i=1}^n \ind{X_i \geq m_0} R_i \]
	 where the ranks $R_1, \ldots, R_n$ of the sample $X_1, \ldots, X_n$ are ordered along $|X_i - m_0|$, i.e. the distance to the proposed median defines the ranks. So in the end, we sum those ranks of the sample that exceed the median. Under $H_0$, we can expect $X_i \geq m_0$ to be true approximately half of the time, and the ranks are uniformly distributed among the sample, so like before under $H_0$ we can expect $W$ to hover around a certain number, we reject if $|W| > c$ (or something else depending on $H_1$) for an appropriate critical value $c$, or we compute the p-value instead with simulations.\\
	 
	 Like before, we can compute means and variances, they come out to be
	 \[ EW = \frac{n(n+1)}{4}, \qquad \operatorname{Var}(W) = \frac{n(n+1)(2n+1)}{24} \]
	 
	 \subsubsection{Ties}
	 What if a tie occurs when computing ranks? We can ignore the problem in two ways:
	 \begin{itemize}
	 	\item Remove the entries from the sample (especially useful when testing $H_0 : m = 0$ v.s. $H_1: m > 0$ and the ties are zero's)
	 	\item Find all the ranks that would be distributed amongst the ties, e.g. $R_1, \ldots, R_i, \ldots, R_j, \ldots, R_n$ are the ranks of the sample, and $R_i, \ldots, R_j$ give ties, then find $\tilde R := (R_i + \ldots + R_j)/(j-i)$ (the arithmetic mean) and assign that rank to every tied entry.
	 \end{itemize}
	 
	 \subsection{Summary}
	 \begin{remark}[Permutation test]
	 	A permutation test finds whether two samples are of the same distribution.\\
	 	Advantages:
	 	\begin{itemize}
	 		\item The test is exact: we can compute the quantiles under $H_0$ analytically, the test statistic is Bernoulli distributed
	 		\item Can be used for lots of test statistics instead of just $T := \bar X - \bar Y$
	 		\item Because of the above we can use many hypotheses
	 	\end{itemize}
	 	Disadvantages:
	 	\begin{itemize}
	 		\item Computationally expensive to compute the quantiles of $B(n, 1/2)$ for large $n$, as well as computing permutations (for large $n$, the amount of permutations grows superexponentially)
	 		\item Assumes that the distribution is invariant under permutation
	 	\end{itemize}
	 \end{remark}
	 \begin{remark}[Sign test]
	 	A sign test can check if the median significantly differs from a proposed median.\\
	 	Advantages: 
	 	\begin{itemize}
	 		\item Easy to compute
	 		\item Works for any data for which a sign can be computed (not necessarily numerical)
	 	\end{itemize}
	 	Disadvantage:
	 	\begin{itemize}
	 		\item Low power
	 		\item Something about paired samples?
	 	\end{itemize}
	 \end{remark}
	 \begin{remark}[Runs test]
	 	A runs test finds whether two samples are of the same distribution.\\
	 	Adv. \textbf{TODO}
	 \end{remark}
	 
	 \newpage
	 
	 \section{The bootstrap}
	 Let $\hat \theta$ be an estimator of $\theta \in \Theta \subseteq \R$. We would like to quantify the uncertainty of $\hat \theta$ given:
	 \begin{itemize}
	 	\item Realizations $(x_1, \ldots, x_n)$ of an iid sample $X_1, \ldots, X_n$ of random variables
	 	\item NO access to the distribution of $\hat \theta(X_1, \ldots, X_n)$ as an r.v.
	 	\item No access to any other data
	 	\item No parametric assumptions
	 \end{itemize}
	 \begin{definition}
	 	Given realizations $x_1, \ldots, x_n$ of an iid sample $X_1, \ldots, X_n$ and an estimator $\hat \theta$, we say that a \textbf{bootstrap sample} of $x_1, \ldots, x_n$ is an ordered tuple $x_1^*, \ldots, x_n^*$ of which each entry is contained in $\{x_1, \ldots, x_n\}$. I.e. it is a draw \textbf{with replacement} of $n$ items from $x_1, \ldots, x_n$.\\
	 	A \textbf{bootstrap estimate} $\hat \theta^* \in \R$ is then the realization of $\hat \theta$ under $x_1^*, \ldots, x_n^*$, i.e.
	 	\[ \hat \theta^* = \hat \theta(x_1^*, \ldots, x_n^*) \]
	 \end{definition}
	 \begin{definition}[Empirical distribution function]
	 	Let $X = (x_1, \ldots, x_n) \in \R^n$ be some sort of sample. We say that its associated \textbf{empirical distribution function} (EDF) is a function $\hat F_{X}(x) : \R \to [0, 1]$ (a type of cdf) such that
	 	\[ \hat F_{X}(x) := \frac{1}{n} \sum_{i=1}^n \ind{x_i \leq x} \]
	 \end{definition}
	 \begin{remark}
	 	We can also define a \textbf{family of random variables}, namely $\hat F_X(x)$ for each $x \in \R$ given an i.i.d. random sample $X_1, \ldots, X_n$ such that
	 	\[ \hat F_{X}(x) := \frac{1}{n} \sum_{i=1}^n \ind{X_i \leq x} \]
	 	almost the same thing, conflated very commonly with the EDF, the difference is that this is a family of random variables and the other simply a function $\R \to [0, 1]$.
	 \end{remark}
	 \begin{remark}
	 	In a measure theoretic sense, $\hat F_X$ is also a random variable, i.e. an estimator for $F$, the true cdf of $X_1, \ldots, X_n$. Notice that
	 	\[ E\hat F_{X}(x) = \frac{1}{n} \sum_{i=1}^n P(X_i \leq x) = \frac{1}{n} n F(X_i \leq x) = F(X_i \leq x) \]
	 	so in a sense, $E \hat F_{X} = F$, so $\hat F$ is unbiased for $F$. woah. And probably with the appropriate assumptions on $F$ we get that it is also consistent.
	 \end{remark}
	 \begin{remark}
	 	Given $\hat \theta_1^*, \ldots \hat \theta_B^*$, $B$ bootstrap estimates for $\theta$, we can now use the EDF to reason about the distribution of $\theta$ empirically. The \textbf{bootstrap assumption} is that the EDF of the bootstrap estimates is a good approximation of the EDF of $X_1, \ldots, X_n$.
	 \end{remark}
	 \begin{definition}
	 	Let $\hat \theta_1^*, \ldots \hat \theta_B^*$, be bootstrap estimates for $\theta$. We define the \textbf{standard bootstrap error}, $\hat{ \text{SE}_\text{boot}} > 0$, through
	 	\[ \hat{ \text{SE}_\text{boot}}^2  := \frac{1}{B-1}\sum_{i=1}^B \left( \hat \theta_i^* - \overline{\oldhat{\theta^*}} \right)^2 \]
	 	i.e. it is like $\hat \sigma^2$ but on bootstrap estimates.
	 \end{definition}
	 \begin{remark}
	 	The bootstrap standard error is a measure for the \textbf{variance} of $\hat \theta$, even though we had absolutely no idea about its underlying distribution before.
	 \end{remark}
	 
	 \subsection{Parametric bootstrap}
	 If we know the underlying model, we can also get a more nuanced bootstrap. Instead of sampling from measurements $X_1, \ldots, X_n \sim P_\theta$ where $\theta \in \Theta \subseteq \R^d$ is a parameter, we can sample $P_{\hat \theta}$ (if possible).
	 \begin{example}
	 	Let $X_1, \ldots, X_n \sim U[0, \theta]$. Let $\theta > 0$. The MLE is $\hat \theta := \max\{X_1, \ldots, X_n\}$. Suppose $\hat \theta$ has a realization, then we sample $X_1^*, \ldots, X_n^* \sim U[0, \hat \theta]$.
	 \end{example}
	 \subsection{Confidence intervals}
	 A lot of quantities that we want to compute can be written as functionals $T$ such that for all probability measures $P$, $T(P) \in \R$ is the measured quantity. For example, we can write
	 \[ \hat \mu(P) := \int_\Omega x dP, \qquad \hat \sigma^2(P) := \int_\Omega (x - \mu)^2 dP \]
	 Let $\hat P$ be the probability measure derived from the EDF of measurements $x_1, \ldots, x_n$. Then we can estimate $T(P)$ using $T(\hat P)$ clearly. How can we now construct confidence intervals for $T(P)$ using the estimator $T(\hat P)$? Well, just like before, we can bootstrap $\hat P$ by sampling it in some sort of way, and getting a bootstrap sample $x_1^*, \ldots, x_n^*$, of which the EDF induced another probability measure $P^*$.\\
	 The \textbf{bootstrap principle} then assumes that we can estimate the distribution of $T(\hat P) - T(P)$ (the error of $T(\hat P)$) using $T(P^*) - T(\hat P)$.
	 \begin{remark}
	 	In order to construct a one-sided $1 - \alpha$ confidence interval for $T(P)$, we could find $c$ s.t.
	 	\[ P(T(\hat P) - T(P) \geq c) = 1 - \alpha \]
	 	because then we have the $(1 - \alpha)$-CI $(-\infty, T(\hat P) + c)$ for $T(P)$. Using the bootstrap principle, we find this $c$ such that
	 	\[ P(T(P^*) - T(\hat P) \geq c) = 1 - \alpha \]
	 \end{remark}
	 And similarly for other intervals. Just pretend that $P^*$ is the actual probability distribution and good things will happen.
	 \begin{remark}
	 	In this framework with functionals $T$, we can similarly find a `bootstrap variance' or `bootstrap standard deviation' like so: if we have \textbf{bootstrap replicates} $T(x^{*1}), \ldots, T(x^{*B})$, we define
	 	\[ \hat{\operatorname{Var}}^*(T(x)) := \frac{1}{B-1} \sum_{i=1}^B \left(T(x^{*i}) - \frac{1}{B} \sum_{j=1}^B T(x^{*j}) \right)^2, \quad \hat{\operatorname{SD}}^*(T(x)) := \sqrt{\hat{\operatorname{Var}}^*(T(x))} \]
	 \end{remark}
	 \begin{example}
	 	We apply the bootstrap to the median. We measure $X_1, \ldots, X_n$, generate $B$ bootstrap estimates $M_n^{*i}$ for $M_n$, where $i = 1, \ldots, B$, and like in the above we can compute the bootstrap variance:
	 	\[ \hat{\operatorname{Var}}^*(M_n) = \frac{1}{B-1} \sum_{i=1}^B \left(M_n^{*i} - \frac{1}{B} \sum_{j=1}^B M_n^{*j} \right)^2 \]
	 	And what about a confidence interval? We suppose that due to the CLT and large enough sample size we can use that
	 	\[ M_n^* \sim N(M_n, \hat{\operatorname{Var}}^*(M_n)) \]
	 	By doing so, we can simply use $z$-quantiles.
	 \end{example}
	 \begin{definition}
	 	Let $[a, b] \subseteq \R$ be a $1 - \alpha$ confidence interval for a quantity $T(P)$. We say that the \textbf{converage probability} of the confidence interval is the probability
	 	\[ P(T(P) \in [a, b]) \]
	 \end{definition}
	 \begin{remark}
	 	More generally, let $T(P_n)$ be a quantity of interest. We estimate using observed data $T(\hat P_n)$. If the central limit theorem holds, and $T(\hat P_n)$ is unbiased, (which can be different criteria), we obviously have that
	 	\[ \frac{T(\hat P_n) - T(P_n)}{\sqrt{\operatorname{Var}(T(\hat P_n))}} \to N(0, 1) \]
	 	of which we can easily compute quantiles. The quantity
	 	\[ \operatorname{Var}(T(\hat P_n)) \approx \operatorname{Var}(T(P^*_n)) \]
	 	using the bootstrap (assuming consistency). Under the assumption that $T(P^*_n)$ is normally distributed with mean $T(\hat P)$, we can compute approximate confidence intervals
	 	\begin{align*}
	 		1 - \alpha &\approx P\left(\frac{T(\hat P_n) - T(P_n)}{\sqrt{\operatorname{Var}(T(\hat P_n))}} \in [-z_{1 - \alpha/2}, z_{1 - \alpha/2}]\right) \\
	 		&\approx P\left(\frac{T(\hat P_n) - T(P_n)}{\underbrace{\sqrt{\operatorname{Var}(T(P_n^*))}}_{=: \sigma^*_n}} \in [-z_{1 - \alpha/2}, z_{1 - \alpha/2}]\right) \\
	 		&= P\left(T(\hat P) - z_{1 - \alpha/2}\sigma^*_n \leq T(P) \leq T(\hat P) + z_{1 - \alpha/2}\sigma^*_n\right)
	 	\end{align*}
	 \end{remark}
	 \begin{definition}
	 	The \textbf{normal bootstrap interval} for an estimator $T(\hat P_n)$ that is \textbf{asymptotically normal} given bootstrap samples
	 	\[ \sigma^*_n = \sqrt{\operatorname{Var}(T(P^*_n))} \]
	 	is the $1 - \alpha$ approximate confidence interval (for $T(p_n)$)
	 	\[ [T(\hat P) - z_{1 - \alpha/2}\sigma^*_n, T(\hat P) + z_{1 - \alpha/2}\sigma^*_n] \]
	 \end{definition}
	 \begin{remark}
	 	Recall that for a certain $X \sim F$ where $F$ has quantile function $Q : (0, 1) \to \R$ satisfies (with $\alpha < \beta$)
	 	\[ P(Q(\alpha) \leq X \leq Q(\beta)) = P(X \leq Q(\beta)) - P(X \leq Q(\alpha)) = \beta - \alpha \]
	 	So that in particular
	 	\[ P(Q(\alpha / 2) \leq X \leq Q(1 - \alpha/2)) = 1 - \alpha \]
	 	which can be used to make confidence intervals. Note that if $EX = 0$ and the pdf is symmetric, then it holds that $Q(\alpha/2) = -Q(1 - \alpha/2)$, but not in general.
	 \end{remark}
	 \begin{definition}
	 	The \textbf{basic pivot bootstrap interval} for an estimator $T(\hat P_n)$ such that
	 	\[ T(\hat P_n) - T(P_n) \approx T(P_n^*) - T(\hat P_n) \tag{asymptotically} \]
	 	is constructed by taking the quantile function $Q$ of $T(P_n^*) - T(\hat P_n)$, and then taking the $1 - \alpha$ approximate confidence interval (for $T(P_n)$)
	 	\[ [T(\hat P_n) - Q(1 - \alpha/2), T(\hat P_n) - Q(\alpha/2)] \]
	 \end{definition}
	 \begin{definition}
	 	The \textbf{percentile bootstrap interval} for an estimator $T(\hat P_n)$ such that
	 	\[ T(P_n^*) \approx T(\hat P_n) \tag{asymptotically} \]
	 	is constructed by taking the quantile function $Q$ of $T(P_n^*)$, and then taking the $1 - \alpha$ approximate confidence interval (for $T(P_n)$)
	 	\[ [Q(\alpha/2), Q(1 - \alpha/2)] \]
	 \end{definition}
	 \begin{remark}
	 	These intervals are cool and all, but do they achieve good coverage? This is hard to compute, and generally speaking we can only talk about asymptotic behaviour, but we can do simulations. It turns out that the coverage is usually not very close to the desired $1 - \alpha$ expected probability. This is because the above intervals may not come from a distribution that is asymptotically close to the real distribution. There is simply a bias that occurs for lower sample sizes, in the form of variance. `Dividing out' the variance would solve this problem, which leads to \text{studentized bootstrap intervals}.
	 \end{remark}
	 \begin{definition}
	 	The \textbf{studentized bootstrap interval} for an estimator $T(\hat P_n)$ of $T(P_n)$ given a standard deviation operator
	 	\[ \sigma(P) := \sqrt{\operatorname{Var}(T(P))} \]
	 	and a quantile function $Q : (0, 1) \to \R$ of the quantity
	 	\[ \frac{T(P_n^*) - T(\hat P)}{\sigma(P_n^*)} \]
	 	is the $1 - \alpha$ approximate confidence interval (for $T(p_n)$)
	 	\[ [T(\hat P) - Q(1 - \alpha/2) \sigma(\hat P_n), T(\hat P) - Q(\alpha/2) \sigma(\hat P_n)] \]
	 \end{definition}
	 \begin{remark}
	 	The above follows from
	 	\begin{align*}
	 		1 - \alpha &\approx P\left(Q(\alpha / 2) \leq \frac{T(P_n^*) - T(\hat P)}{\sigma(P_n^*)} \leq Q(1 - \alpha/2) \right) \\
	 		&\approx P\left(Q(\alpha / 2) \leq \frac{T(\hat P_n) - T(P)}{\sigma(\hat P_n)} \leq Q(1 - \alpha/2) \right) \\
	 		&= P\left(T(\hat P_n) - Q(\alpha / 2)\sigma(\hat P_n) \geq T(P) \geq T(\hat P_n) - Q(1 - \alpha/2)\sigma(\hat P_n) \right)
	 	\end{align*}
	 \end{remark}
	 \subsection{Bootstrap failure}
	 \begin{example}
	 	Let $X_1, \ldots, X_n \sim U[0, 1]$ be i.i.d, let $T_n := \min \{X_1, \ldots, X_n\}$. Then (in distribution)
	 	\[ nT_n \to \operatorname{Exp}(1) \]
	 	\begin{proof}
	 		We need to check that $F_n \to F$ pointwisely, where $F_n$ is the cdf of $nT_n$, and $F$ is the cdf of a $\operatorname{Exp}(1)$ variable, the latter is given by
	 		\[ F(x) = 1 - e^{-x} \]
	 		We compute
	 		\begin{align*}
	 			F_n(x) &= P(nT_n \leq x) \\
	 			&= P\left(T_n \leq \frac{x}{n} \right) \\
	 			&= 1 - P\left(T_n > \frac{x}{n} \right) \\
	 			&= 1 - \prod_{i=1}^n P\left(X_i > \frac{x}{n} \right) \\
	 			&= 1 - \prod_{i=1}^n \left(1 - P\left(X_i \leq \frac{x}{n} \right) \right) \\
	 			&= 1 - \prod_{i=1}^n \left(1 - \frac{x}{n} \right)\\
	 			&= 1 - \left(1 - \frac{x}{n} \right)^n \to 1 - e^{-x}
%	 			&= 1 - \prod_{i=1}^n \ind{\frac{x}{n} \leq 1} \frac{n - x}{n} \\
	 		\end{align*}
	 		There is more subtlety, though, because the c.d.f is actually given by
	 		\[ P\left(X_i \leq \frac{x}{n} \right) = \begin{cases*}
	 			0 & if $\frac{x}{n} < 0$ \\
	 			\frac{x}{n} & if $0 \leq \frac{x}{n} \leq 1$ \\
	 			1 & if $\frac{x}{n} > 1$
	 		\end{cases*} \]
	 		So we restrict to $x > 0$ and then for large enough $n$ we have $\frac{x}{n} < 1$, and in the case that $x < 0$ the cdfs trivially coincide.
	 	\end{proof}
	 	The result is that for appropriate $n$, $nT_n$ is approximately $\operatorname{Exp}(1)$, and therefore
	 	\[ 1 - e^{-nx} \approx P(nT_n \leq nx) = P(T_n \leq x) \]
	 	and therefore $T_n$ is approximately $\operatorname{Exp}(n)$.\\
	 	Now suppose one bootstrap sample is $X_1^*, \ldots, X_n^*$, and we compute $T_n^* = \min \{X_1^*, \ldots, X_n^*\}$. We get that for all $i = 1, \ldots, n$
	 	\[ P(X_i^* = M_n) = \frac{1}{n} \]
	 	because we are drawing with replacement, assuming the original data $X_1, \ldots, X_n$ contains no duplicates (which happens wp 1). Therefore
	 	\[ P(X_1^* \neq M_n, \ldots, X_n^* \neq M_n) = P(X_1^* \neq M_n)^n = \left(1 - \frac{1}n\right)^n \to e^{-1} \]
	 	whereby
	 	\[ P(M_n^* = M_n) \to 1 - e^{-1} \neq 0 \]
	 \end{example}
	 \subsection{Hypothesis testing}
	 Let's say we want to test
	 \[ H_0 : \theta \in \Theta_0 \qquad \text{v.s.} \qquad H_1 : \theta \notin \theta_0 \]
	 with a test statistic $T : \R^n \to \R$, and our rejection region is $T > c_\alpha$ for appropriate critical value $c_\alpha$, such that
	 \[ \forall \theta \in \Theta_0 : P_{\theta}(T(X) > c_\alpha) \leq \alpha  \iff \sup_{\theta \in \Theta_0} P_{\theta}(T(X) > c_\alpha) \leq \alpha \]
	 \begin{remark}[p-value]
	 	By definition, the $p$-value is the smallest significance level for which we would still reject. So for measured data $x \in \R^n$
	 	\[ p(x) = \inf \{\alpha \in (0, 1) \mid T(x) > c_\alpha \} \]
	 	so that
	 	\[ T(x) > c_\alpha \iff p(x) \leq \alpha \]
	 	But also by monotonicity and the decreasing property of the test
	 	\[ T(x) > c_\alpha \iff \sup_{\theta \in \Theta_0} P_\theta(T(X) > T(x)) < \sup_{\theta \in \Theta_0} P_\theta(T(X) > c_\alpha) \leq \alpha \]
	 	so that
	 	\[ p(x) = \inf \{\alpha \in (0, 1) \mid \sup_{\theta \in \Theta_0} P_\theta(T(X) > T(x)) < \alpha \} = \sup_{\theta \in \Theta_0} P_\theta(T(X) > T(x)) \]
	 	Intuituively, this means that the \textbf{smallest level for which we would still reject} is the \textbf{upper bound of the probability of measuring} $T(x)$ or higher.\\
	 \end{remark}
	 \textbf{Continuing}, we find an estimator $\hat \theta$ and use the approximate distribution $P_{\hat \theta} =: \hat P_\theta$. We generate i.i.d. bootstrap samples $(X^{*1}, \ldots, X^{*B}) \in (\R^n)^B$. We define the estimator of the $p$-value:
	 \[ \hat p := \hat P(T(X) \geq T(x)) = \frac{1}{B} \sum_{i=1}^B \ind{T(X^{*i}) \geq T(x)} \]
	 An alternative:
	 \[ \hat p := \frac{1 + \sum_{i=1}^B \ind{T(X^{*i}) \geq T(x)}}{B + 1}  \]
	 which behaves asymptotically identical but slightly overestimates the $p$-value for small values of $B$, to not reject too often.
	 
	 \subsection{Asymptotic behaviour}
	 \begin{theorem}
	 	Let $(X_n)_{n \in \N^*}$ be a sequence of i.i.d. random variables in $\R^d$. Let $\mu := EX_1$ and $\Sigma := \operatorname{Cov}(X_1)$. Let for every $n \in \N*$, there be a iid bootstrap sample of $X_1, \ldots, X_n$, namely for empirical distribution $\hat P_n$,
	 	\[ X^*_{n, 1}, \ldots, X^*_{n, n} \sim \hat P_n \]
	 	(draws with replacement). Define additionally the sequences
	 	\[ \bar X_n := \frac{1}{n} \sum_{i=1}^n X_i, \qquad \bar X^*_n := \frac{1}{n} \sum_{i=1}^n X^*_{n, i} \]
	 	Then the conditional distribution of $\sqrt{n}(\bar X^*_n - \bar X_n)$ on $X_1, \ldots, X_n$ converges weakly to $N(0, \Sigma)$.
	 \end{theorem}
	 \begin{remark}
	 	This can be shown using the multivariate CLT and additionally using that the bootstrap samples are random variables in a sense.
	 \end{remark}
	 \begin{remark}
	 	Surprisingly enough, not much else can be said, namely if the quantities that we estimate give consistent estimators. As we have seen before, this may as well not be the case.
	 \end{remark}
	 
	 \subsection{Regression}
	 Simple linear regression model is as follows. We have $X_1, \ldots, X_n$ and $(Y_1, \ldots, Y_n)$ i.i.d. samples. The hypothesis is that there exist $\varepsilon_1, \ldots, \varepsilon_n$ independent with $E\varepsilon_i = 0$ and $\operatorname{Cov}(X_i, \varepsilon_i)$ for each $i$, such that there exist $\beta_0, \beta_1 \in \R$
	 \[ Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i \]
	 Suppose we have estimators $\hat \beta_1$ and $\hat \beta_0$ typically given by
	 \[ \beta_1 := \hat{\operatorname{Corr}}(X, Y) = \frac{\sum_{i=1}^n (X_i - \bar X)(Y_i - \bar Y)}{\sum_{i=1}^n (X_i - \bar X)^2}, \qquad \beta_0 := \bar Y - \hat \beta_1 \bar X \]
	 The linear model becomes
	 \[ \hat Y_i := \hat \beta_0 + \hat \beta_1 X_i \]
	 The residuals are
	 \[ e_i := Y_i - \hat Y_i \]
	 We want to now use bootstrapping to find variances of $\beta_1, \beta_0$. Simple. We create bootstrap samples $X_1^*, \ldots, X_n^*$ and $Y_1^*, \ldots, Y_n^*$, but we draw with replacement from $(X_1, Y_1), \ldots, (X_n, Y_n)$, so that for each $i, j = 1, \ldots, n$
	 \[ P(X^*_i = X_j, Y^*_i = Y_j) = \frac{1}{n} \]
	 If we repeat this procedure $B$ times we get bootstrap samples $(X_i^{*j})$ for every $i=1, \ldots, n$, $j = 1, \ldots, B$. For each sample, find $\hat \beta_0^{*j}$ and $\hat \beta_1^{*j}$, and now we have empirical distributions. We can compute $\operatorname{Var}(\beta_i^*)$ as well as $\operatorname{MSE}$, and using the quantiles of $\hat \beta_i^*$ we can even construct a confidence interval:
	 \[ [2\hat \beta_i - q_{i, 1 - \alpha/2}^*, 2\hat \beta_i - q_{i, \alpha/2}^*] \]
	 
	 \newpage
	 
	 \section{Bayesian statistics}
	 \begin{definition}
	 	\letprobalways Let $A, B \in \Sigma$ be events. Then the \textbf{conditional probability} of $A$ given $B$ (when $P(B) > 0$) is given by
	 	\[ P(A|B) := \frac{P(A \cap B)}{P(B)} \]
	 \end{definition}
	 \begin{theorem}[Bayes' theorem]
	 	\letprobalways Let $A, B \in \Sigma$ be events, such that $P(A) > 0$ and $P(B > 0)$, then
	 	\[ P(B|A) = \frac{P(A|B)P(B)}{P(A)} \]
	 \end{theorem}
	 \begin{definition}
	 	Let $X, Y$ be random variables with densities $f_X, f_Y$ respectively, and joint density $f_{X, Y} : \R^2 \to \R$. We define the \textbf{conditional density} of $X$ given $Y$ as a function $f_{X|Y} : D \to \R$ such that
	 	\[ f_{X|Y}(x|y) := \frac{f_{X, Y}(x, y)}{f_Y(y)} \]
	 	where $D \subseteq \R^2$ is the domain where $(x, y) \mapsto f_Y(y)$ is nonzero.
	 \end{definition}
	 \begin{theorem}[Bayes' theorem for densities]
	 	Let $X, Y$ be random variables with densities $f_X, f_Y$ respectively, and joint density $f_{X, Y} : \R^2 \to \R$. Then for every $x, y \in \R$ such that $f_X(x) \neq 0$ and $f_Y(y) \neq 0$ we have
	 	\[ f_{Y|X}(y|x) = \frac{f_{X|Y}(x|y)f_Y(y)}{f_X(x)} \]
	 \end{theorem}
	 \begin{remark}
	 	Notice that very importantly,
	 	\[ f_{X}(x) = \int_\R f_{X|Y}(x|y) f_Y(y) \, dy \]
	 	such that we can write Bayes' theorem as
	 	\[ f_{Y|X}(y|x) = \frac{f_{X|Y}(x|y)f_Y(y)}{\int_\R f_{X|Y}(x|y) f_Y(y) \, dy} \]
	 	which immediately makes it clear why it is a valid probability density.
	 \end{remark}
	 Let $X_1, \ldots, X_n \sim P_\theta$ be an i.i.d. sample, for which $\theta$ is unknown but $\Theta \subseteq \R$, the parameter space, is known, and must be a compact set. Suppose that $\theta$ was drawn from a distribution with density $\pi$, called the \textbf{prior}. Then the \textbf{posterior} given data is the conditional density $\pi(\cdot|X_1, \ldots, X_n)$. It can be computed through:
	 \[ \pi(\theta|X_1, \ldots, X_n) = \frac{f_{X|\theta}(X_1, \ldots, X_n|\theta) \pi(\theta)}{\int_\Theta f_{X|\theta}(X_1, \ldots, X_n|\theta) \pi(\theta) \, d\theta} \]
	 The posterior gives you in a sense the likelihood of $\theta$ given our measurements.
	 \begin{example}
	 	If the prior is chosen uniformly, let's say $\pi(\theta) = \ind{0 \leq \theta \leq 1}$, we get
	 	\[ \pi(\theta|X_1, \ldots, X_n) = \frac{f_{X|\theta}(X_1, \ldots, X_n|\theta) \pi(\theta)}{\int_\Theta f_{X|\theta}(X_1, \ldots, X_n|\theta) \pi(\theta) \, d\theta} = \frac{f_{X|\theta}(X_1, \ldots, X_n|\theta)}{\int_0^1 f_{X|\theta}(X_1, \ldots, X_n|\theta) \, d\theta} = f_{X|\theta}(X_1, \ldots, X_n|\theta) \]
	 	i.e. we get our familiar friend the likelihood back for the posterior.
	 \end{example} 
	 \begin{definition}
	 	Let $\alpha, \beta > 0$. We say that some random variable $X$ with density $f$ is \textbf{beta distributed} if
	 	\[ f(x) = \frac{x^{\alpha - 1} (1-x)^{\beta - 1}}{\int_0^1 y^{\alpha - 1} (1-y)^{\beta - 1} \, dy} \ind{0 \leq x \leq 1} \]
	 	We write $X \sim \operatorname{Beta}(\alpha, \beta)$.
	 \end{definition}
	 \begin{example}
	 	Let $Y \sim U[0, 1]$, let $X_1, \ldots, X_n \sim U[0, 1]$ be i.i.d. and let
	 	\[ Z := \sum_{i=1}^n \ind{X_i \leq Y} \]
	 	Then what is the distribution of $Z|Y$? If we know, then we can compute an updated distribution $Y|Z$ for $Y$, that gives you more information about $Y$ given $Z$. Since $X_i$ are chosen uniformly, we know that $X_i|Y=y \sim \operatorname{Ber}(y)$, so the sum is $B(n, y)$. Then
	 	\begin{align*}
	 		f_{Y|Z}(y|z) &= \frac{f_{Z|Y}(z|y) f_Y(y)}{\int_\R f_{Z|Y}(z|y) f_Y(y) \, dy} \\
	 		&= \frac{\binom{n}{z} y^{z} (1 - y)^{n - z} \ind{0 \leq y \leq 1}}{\int_\R \binom{n}{z} y^{z} (1 - y)^{n - z} \ind{0 \leq y \leq 1} dy} \\
	 		&= \frac{y^{z} (1 - y)^{n - z} \ind{0 \leq y \leq 1}}{\int_\R y^{z} (1 - y)^{n - z} \ind{0 \leq y \leq 1} dy}
	 	\end{align*}
	 	which implies that $Y|Z \sim \operatorname{Beta}(Z + 1, n + 1 - Z)$.
	 \end{example}
	 \begin{example}
	 	Let's say we have a linear regression model $Y = X\beta + \varepsilon$, where
	 	\[ \varepsilon \sim N(0, \sigma^2 \id) \]
	 	We observe $Y$, know $X$ and $\sigma^2 > 0$. Can we use Bayesian statistics to find an approximation of $\beta$? Suppose that $\beta \sim N(0, \tau^2 \id)$ with $\tau^2 > 0$, this is the prior. Recall that $Y \sim N(X\beta, \sigma^2 \id)$. What is the distribution of $\beta|Y$? Well
	 	\begin{align*}
	 		f_{Y|\beta}(y, \beta) f_\beta(\beta) &= \frac{1}{\sqrt{2 \pi \sigma^n}} \exp\left(-\frac{\norm{y - X\beta}^2}{2\sigma^2}\right) f_\beta(\beta) \\
	 		&= \frac{1}{\sqrt{4 \pi^2 \sigma^n \tau^n}} \exp\left(-\frac{\norm{y - X\beta}^2}{2\sigma^2} - \frac{\norm{\beta}^2}{2\tau^2}\right) \\
	 		&\propto \exp\left(-\frac{\tau^2\norm{y}^2 - 2\tau^2\inprod{X^T y, \beta} + (\tau^2 + \sigma^2)\norm{\beta}^2}{2(\sigma \tau)^2}\right) \\
	 		&\propto \exp\left(-\frac{\tau^2\norm{y}^2 - \inprod{2\tau^2 X^T y + (\tau^2 + \sigma^2) \beta, \beta}}{2(\sigma \tau)^2}\right) \\
	 		&\propto \exp\left(-\frac{\norm{y}^2 - \inprod{2 X^T y + (1 + (\sigma/\tau)^2) \beta, \beta}}{2\sigma^2}\right)
	 	\end{align*}
	 \end{example}
	 \subsection{Estimation}
	 Let's say we have a posterior $\pi(\theta|X)$, how do we estimate?
	 \begin{definition}[Posterior mean]
	 	The \textbf{posterior mean estimator} assigns to each measurement $X$ the estimation
	 	\[ \hat \theta(X) := \int_\Theta \theta \pi(\theta|X) \, d\theta \]
	 	i.e. the mean of the posterior. woah.
	 \end{definition}
	 \begin{definition}[Posterior mode]
	 	The MAP or posterior mode estimator assigns to each measurement $X$ the estimation
	 	\[ \hat \theta(X) := \operatorname*{argsup}_{\theta \in \Theta} \pi(\theta|X) \]
	 \end{definition}
	 \begin{definition}
	 	A \textbf{loss function} is a metric on $\Theta$ without triangle inequality. That is, it is a map $L : \Theta^2 \to [0, \infty)$ such that for all $\theta_1, \theta_2 \in \Theta$ we have
	 	\begin{itemize}
	 		\item $L(\theta_1, \theta_2) = 0 \iff \theta_1 = \theta_2$
	 		\item $L(\theta_1, \theta_2) = L(\theta_2, \theta_1)$
	 	\end{itemize}
	 \end{definition}
	 \begin{definition}
	 	The following loss functions are common:
	 	\begin{itemize}
	 		\item Absolute loss: $L(\theta, \hat \theta) := |\theta - \hat \theta|$
	 		\item Quadratic loss: $L(\theta, \hat \theta) := \left(\theta  - \hat \theta\right)^2$
	 		\item All-or-nothing loss: $L(\theta, \hat \theta) := \ind{\theta \neq \hat \theta}$
	 	\end{itemize}
	 \end{definition}
	 \begin{definition}
	 	The \textbf{Bayes risk} of a loss function $L$ given a parameter $\theta$ a random variable in $\Theta$ and an estimate $\hat \theta \in \Theta$, is defined to be the expectation of the loss $L$, that is, the Bayes risk $R(\theta, \hat \theta)$ is defined by
	 	\[ R(\theta, \hat \theta) := E[L(\theta, \hat \theta)] \]
	 \end{definition}
%	 \begin{example}
%	 	For a quadratic loss, we can find $R$ as follows:
%	 	\begin{align*}
%	 		R(\theta, \hat \theta) = E[(\theta - \hat \theta)^2] = \operatorname{MSE}(\hat \theta)
%	 	\end{align*}
%	 \end{example}
	 \begin{lemma}
	 	With the quadratic loss function, the estimator that minimizes the Bayes risk is the posterior mean.
	 \end{lemma}
	 \begin{proof}
	 	Notice that
	 	\[ R(\theta, \hat \theta) = E[(\hat \theta - \theta)^2] = \int_\Theta (\hat \theta - \theta)^2 \pi(\theta|X) \, d\theta \]
	 	Minimizing amounts to finding a critical point first, so let's do that
	 	\[ \frac{d}{d\hat \theta} \int_\Theta (\hat \theta - \theta)^2 \pi(\theta|X) \, d\theta = 2\int_\Theta (\hat \theta - \theta) \pi(\theta|X) \, d\theta = 0 \]
	 	which gives
	 	\[ \int_\Theta \hat \theta \pi(\theta|X) \, d\theta = \hat \theta = \int_\Theta \theta \pi(\theta|X) \, d\theta \]
	 	is the only critical point. The second derivative is positive, so we have a convex function $R$ for fixed $\theta$, which attains a global minimum at the critical point. Notice that said critical point is the posterior mean.
%	 	Recognize that the posterior mean $\hat \theta := E[\pi(\theta|X)]$
	 \end{proof}
	 \begin{lemma}
	 	For the absolute loss, the best estimator is the posterior median, i.e. the median of $\theta|X$.
	 \end{lemma}
	 \begin{proof}
	 	Suppose that $\Theta = [a, b]$ for $a < b$ (can be generalized for finite unions of compact intervals and even improper integrals by taking limits). We require that $\hat \theta \in (a, b)$, otherwise we get slightly into trouble (but the $R$ is a.s. not minimized then)
	 	\begin{align*}
	 		R(\theta, \hat \theta) &= E[|\hat \theta - \theta|] \\
	 		&= \int_\Theta |\hat \theta - \theta| \pi(\theta|X) \, d\theta \\
	 		&= \int_a^{\hat \theta} (\hat \theta - \theta) \pi(\theta|X) \, d\theta - \int_{\hat \theta}^b (\hat \theta - \theta) \pi(\theta|X) \, d\theta
	 	\end{align*}
	 	We differentiate using the Leibnitz rule
	 	\begin{align*}
	 		\frac{\partial}{\partial \hat \theta} R(\theta, \hat \theta) &= \frac{d}{d \hat \theta} \left(\int_a^{\hat \theta} (\hat \theta - \theta) \pi(\theta|X) \, d\theta - \int_{\hat \theta}^b (\hat \theta - \theta) \pi(\theta|X) \, d\theta\right) \\
	 		&= \frac{d}{d \hat \theta}\int_a^{\hat \theta} (\hat \theta - \theta) \pi(\theta|X) \, d\theta + \frac{d}{d \hat \theta}\int^{\hat \theta}_b (\hat \theta - \theta) \pi(\theta|X) \, d\theta \\
	 		&= \int_a^{\hat \theta} \pi(\theta|X) \, d\theta + \int^{\hat \theta}_b \pi(\theta|X) \, d\theta = 0
	 	\end{align*}
	 	This clearly gives you the median of $\theta|X$, now what is remaining is to show that it is indeed a minimum, but once again the original function is convex, using the Leibnitz rule again we get a second derivative of $2 \pi(\theta|X) > 0$.
%	 	Since $\Theta$ is a compact set, we can split it into an almost disjoint union of compact sets, $\Theta = \Theta_+ \cup \Theta_-$, such that for all $\theta \in \Theta_+$ we have $\hat \theta - \theta \geq 0$, and for $\Theta_-$ we have $\hat \theta - \theta \leq 0$, so that
%	 	\[ \int_\Theta |\hat \theta - \theta| \pi(\theta|X) \, d\theta = \int_{\Theta_+} (\hat \theta - \theta) \pi(\theta|X) \, d\theta - \int_{\Theta_-} (\hat \theta - \theta) \pi(\theta|X) \, d\theta \]
	 \end{proof}
	 \begin{lemma}
	 	For discrete posteriors and the all-or-nothing loss, the posterior mode is the Bayes estimator.
	 \end{lemma}
	 \begin{definition}
	 	Let $\pi(\cdot|X)$ be the posterior. Let $\Theta$ be the parameter space. A $(1-\alpha)$ \textbf{credible set} $C$ satisfies
	 	\[ \int_C \pi(\theta|X) \, d\theta \geq 1 - \alpha \]
	 	Notice that this definition depends on data sample $X$.
	 \end{definition}
	 \begin{remark}
	 	Compare credible sets with confidence intervals: if $C$ is a credible set then
	 	\[ P(\underbrace{\theta}_{\text{an r.v.}} \in C|\underbrace{X}_\text{data dependency}) \geq 1 - \alpha \]
	 	If $C$ is a confidence interval then
	 	\[ \inf_{\theta \in \Theta} P_\theta(\underbrace{\theta}_{\text{fixed value}} \in C) \geq 1 - \alpha \]
	 \end{remark}
	 \begin{definition}
	 	Let $C$ be a credible set, we say that the \textbf{frequentist coverage} of $C$ given \textbf{true value} $\theta \in \Theta$ is $P_\theta(\theta \in C)$. We say that the (overall) frequentist coverage of $C$ is the infimum over all $\theta \in \Theta$.
	 \end{definition}
	 \begin{example}
	 	Let $\theta \sim N(0, \tau^2)$ be the prior, let $X|\theta \sim N(\theta, 1)$, then we find the posterior:
	 	\[ \theta|X \sim N\left(\frac{\tau^2}{1 + \tau^2}X, \frac{\tau^2}{1 + \tau^2}\right) \]
	 	Note that in this case the posterior mean is then $\frac{\tau^2}{1 + \tau^2}X$. Note that as $\tau \to \infty$ the posterior mean approaches $X$, which is an unbiased estimator of $\theta$. Anyway, we find a credible set. We want to find $C$ such that
	 	\[ P(\theta \in C|X) = 1 - \alpha \]
	 	Suppose $C$ is of the form $\left[\frac{\tau^2}{1 + \tau^2}X - c, \frac{\tau^2}{1 + \tau^2}X + c\right]$, where $c$ is to be determined, then for standard normal $Z$
	 	\begin{align*}
	 		P(\theta \in C|X) &= P\left(-c \sqrt{\frac{1 + \tau^2}{\tau^2}} \leq Z \leq c \sqrt{\frac{1 + \tau^2}{\tau^2}}\right) \\
	 		&= P\left(Z \leq c \sqrt{\frac{1 + \tau^2}{\tau^2}}\right) - P(Z \leq -c \sqrt{\frac{1 + \tau^2}{\tau^2}}) \\
	 		&= 2\Phi(c \sqrt{\frac{1 + \tau^2}{\tau^2}}) - 1 = 1 - \alpha \\
	 		\implies \Phi\left(c \sqrt{\frac{1 + \tau^2}{\tau^2}}\right) &= 1 - \alpha/2 \\
	 		\implies c &= z_{1 - \alpha/2} \sqrt{\frac{\tau^2}{1 + \tau^2}}
	 	\end{align*}
	 	Notice that this is very very similar to the standard confidence interval, just slightly narrower and centered differently.
	 \end{example}
	 \newpage
	 \section{Time series}
	 \subsection{Definitions}
	 \begin{definition}
	 	A \textbf{time series model} is a sequence $(X_t)_{t \in \N^*}$ of random variables, that are not necess
	 \end{definition}
	 \begin{definition}
	 	A time series $(X_t)_{t \in \N^*}$ is said to be \textbf{white noise} if for all $i \neq j$ we have $\operatorname{Cov}(X_i, X_j) = 0$, and have the same first two moments. We write $(X_t) \sim \operatorname{WN}(\mu, \sigma^2)$.
	 \end{definition}
	 \begin{definition}
	 	A \textbf{random walk} is a time series that is the sum of i.i.d. random variables.
	 \end{definition}
	 \begin{definition}
	 	The \textbf{autocovariance function} of a timeseries $(X_t)_{t \in \N^*}$ is a function $\gamma_X : \N \times \N \to \R$ such that
	 	\[ (s, t) \mapsto \operatorname{Cov}(X_s, X_t) = E[(X_s - \mu_X(s))(X_t - \mu_X(t))] \]
	 	where $\mu_X : \N \to \R$ is the \textbf{mean function} of $X$, given by $t \mapsto E[X_t]$.
	 \end{definition}
	 \subsection{Stationarity}
	 \begin{definition}
	 	A time series is called \textbf{stationary} if
	 	\begin{itemize}
	 		\item $\mu_X$ is constant
	 		\item The forward autocovariance is stationary, that is, for all $\tau$, $\gamma_X(t + \tau, t)$ is independent of $t$ / constant
	 	\end{itemize}
	 \end{definition}
	 \begin{definition}
	 	A time series is called \textbf{strictly stationary} if the c.d.f is invariant under timeshift, that is, for all $\tau \in \Z$, $n \in \N^*$ and $t_1, \ldots, t_n$, we have
	 	\[ F_{t_1 + \tau, \ldots, t_n + \tau} = F_{t_1, \ldots, t_n} \]
	 	Where
	 	\[ F_{t_1, \ldots, t_n}(x_1, \ldots, x_n) := P(X_{t_1} \leq x_1, \ldots, X_{t_n} \leq x_n) \]
	 \end{definition}
	 \begin{theorem}
	 	The following properties hold
	 	\begin{itemize}
	 		\item If $X_t$ is strictly stationary, then they all come from the same distribution
	 		\item If $X_t$ is strictly stationary, then $(X_t, X_{t+h}) = (X_1, X_{1 + h})$ for all $t, h$
	 		\item Strict stationarity implies stationarity if the second moment of all time series variables $X_t$ is finite
	 		\item The converse is not true
	 		\item An i.i.d. timeseries is strictly stationary
	 	\end{itemize}
%	 	Strictly stationary implies stationary.
	 \end{theorem}
	 \begin{proof}
	 	Consider the case for a certain $i, j$ and $n=1$, then we get
	 	\[ P(X_j \leq x) = F_{i + (j - i)}(x) = F_i(x) = P(X_i \leq x) \implies X_i = X_j \]
	 	This covers the second case as well? Anyway, given this we easily get that $\mu_X$ is constant, and because of that also the autocovariance. The converse not being true can be seen by an example but I don't bother. The last property we prove. Let $(X_t)_{t \in \N^*}$ be an i.i.d. sequence of random variables. Let $t_1, \ldots, t_n$, $x_1, \ldots, x_n$ and $\tau \in \Z$ be arbitrary. Then
	 	\[ P(X_{t_1 + \tau} \leq x_1, \ldots, X_{t_n + \tau} \leq x_n) = \prod_{i=1}^n P(X_{t_i + \tau} \leq x_i) = \prod_{i=1}^n P(X_{t_i} \leq x_i) = P(X_{t_1} \leq x_1, \ldots, X_{t_n} \leq x_n) \]
	 \end{proof}
	 \begin{definition}
	 	A timeseries $(X_t)$ is said to be $q$-dependent for some $q \in \N$ if for all $s, t \in \N^*$ we have the implication
	 	\[ |s - t| > q \implies X_s, X_t \text{ are independent} \]
	 \end{definition}
	 \begin{example}
	 	If $X_t = g(Z_t, Z_{t-1}, \ldots, Z_{t-q})$ for some i.i.d. noise series $(Z_t)$, this is $q$-dependent. It is however also strictly stationary by the i.i.d. assumption.
	 \end{example}
	 \subsection{Gaussian processes}
	 \begin{definition}
	 	A timeseries process $(X_t)$ is said to be \textbf{Gaussian} if for every $t_1, \ldots, t_n \in \N^*$, we have that $X := (X_{t_1}, \ldots, X_{t_n})$ is multivariate normal.
	 \end{definition}
	 \begin{lemma}
	 	Let $(X_t)$ be Gaussian, then
	 	\[ X_t \text{ stationary} \implies X_t \text{ strictly stationary} \]
	 	In particular, we already had the reverse implication, so we have equivalence.
	 \end{lemma}
	 \begin{proof}
	 	Let $X_t$ be Gaussian and stationary. Let $t_1, \ldots, t_n \in \N^*$, $\tau \in \Z$, and $x_1, \ldots, x_n \in \R$ be arbitary. We want to show that
	 	\[ F_{t_1, \ldots, t_n} = F_{t_1 + \tau, \ldots, t_n + \tau} \]
	 	Consider that $X := (X_{t_1}, \ldots, X_{t_n}) \sim N(\mu, \Sigma)$ for certain $\mu \in \R^n$ and $\Sigma \in \R^{n \times n}$ a covariance matrix. If we can show that $\tilde X := (X_{t_1 + \tau}, \ldots, X_{t_n + \tau}) \sim N(\mu, \Sigma)$, which is already by the Gaussian assumption multivariate normal (so it suffices to show that mean and covariance coincide), we are done. Consider that
	 	\[ \mu = (EX_{t_1}, \ldots, EX_{t_n}) = (\mu_X(t_1), \ldots, \mu_X(t_n)) = (\tilde \mu, \ldots, \tilde \mu) \]
	 	for a certain $\tilde \mu \in \R$ by weak stationarity. The same reasoning can be used to conclude that $EX = E\tilde X$. Now is just suffices to show that the covariance matrix coindices. Consider that
	 	\[ \Sigma^i_j = \operatorname{Cov}(X_{t_i}, X_{t_j}) = \gamma_X(t_i, t_j) = \gamma_X(t_j + (t_i - t_j), t_j) \]
	 	Let $\tilde \Sigma := \operatorname{Cov}(\tilde X)$, then
	 	\[ \tilde \Sigma^i_j = \operatorname{Cov}(X_{t_i + \tau}, X_{t_j + \tau}) = \gamma_X(t_i + \tau, t_j + \tau) = \gamma_X(t_j + \tau + (t_i - t_j), t_j + \tau) = \gamma_X(t_j + (t_i - t_j), t_j) = \Sigma^i_j \]
	 	so $\tilde \Sigma = \Sigma$ and we are done.
	 \end{proof}
	 \subsection{Auto-things and linear processes}
	 \begin{definition}
	 	Let $(X_t)$ be stationary, then we define the function $\gamma_X : \N \to \R$ by
	 	\[ \tau \mapsto \gamma_X(t + \tau, t) = \operatorname{Cov}(X_{t + \tau}, X_t) \]
	 	which is defined independent of $t$ by stationarity, we might as well take $t=|\tau| + 1$ for well-definedness or smth. We additionally define the \textbf{autocorrelation function}, $\rho_X : \N \to [-1, 1]$, by
	 	\[ \tau \mapsto \frac{\gamma_X(\tau)}{\gamma_X(0)} \]
	 \end{definition}
	 \begin{theorem}
	 	The following hold
	 	\begin{itemize}
	 		\item $\gamma_X(0) \geq 0$
	 		\item For all $\tau \in \Z$, $|\gamma_X(\tau)| \leq \gamma_X(0)$
	 		\item For all $\tau \in \Z$, $\gamma_X(\tau) = \gamma_X(-\tau)$
	 		\item $\gamma_X$ is positive definite, that is for all $a \in \R^n$ and arbitary $n$ we have
	 		\[ \sum_{i=1}^n \sum_{j=1}^n a^i \gamma_X(i - j) a^j \geq 0 \]
	 	\end{itemize}
	 \end{theorem}
	 \begin{proof}~
	 	\begin{itemize}
	 		\item $\gamma_X(0) = \gamma_X(t, t) = \operatorname{Cov}(X_t, X_t) = \operatorname{Var}(X_t) \geq 0$
	 		\item Let $\tau \in \Z$, then
	 		\begin{align*}
	 			|\gamma_X(\tau)| &= |\gamma_X(t + \tau, t)| \\
	 			&= |\operatorname{Cov}(X_{t + \tau}, X_t)| \\
	 			&\leq \sqrt{\operatorname{Var}(X_{t+\tau})\operatorname{Var}(X_{t})} \\
	 			&= \sqrt{\operatorname{Var}(X_{t})\operatorname{Var}(X_{t})} \\
	 			&= \operatorname{Var}(X_t) = \gamma_X(0)
	 		\end{align*}
	 		\item Let $\tau \in \Z$, then
	 		\[ \gamma_X(-\tau) = \gamma_X(t - \tau, t) = \gamma_X(t - \tau + \tau, t + \tau) = \gamma_X(t, t + \tau) = \gamma_X(t + \tau, t) = \gamma_X(\tau) \]
	 	\end{itemize}
	 \end{proof}
	 \begin{definition}
	 	A time series $(X_t)_{t \in \N}$ is called \textbf{a linear process} if there exists white noise $(W_t)_{t \in \Z} \sim \operatorname{WN}(0, \sigma^2)$ and absolutely summable coefficients $(\psi_j)_{j \in \Z}$ and $\mu \in \R$ such that
	 	\[ X_t = \mu + \sum_{j \in \Z} \psi_j W_{t - j} \]
	 \end{definition}
	 \begin{lemma}
	 	A linear process is weakly stationary with
	 	\[ \mu_X(t) = \mu \]
	 	\[ \gamma_X(\tau) = \sigma^2 \sum_{j \in \Z} \psi_j\psi_{\tau + j} \]
	 \end{lemma}
	 \begin{proof}
	 	We simply show the formulae, which directly implies weak stationarity. Consider that
	 	\begin{align*}
	 		\mu_X(t) &= E[X_t] \\
	 		&= \mu + \sum_{j \in \Z} \psi_j E[W_{t - j}] = \mu
	 	\end{align*}
	 	Let $t \in \N$ and $\tau \in Z$, then
	 	\begin{align*}
	 		\gamma_X(\tau) &= \gamma_X(t + \tau, t) \\
	 		&= \operatorname{Cov}(X_{t + \tau}, X_t) \\
	 		&= \operatorname{Cov}\left( \mu + \sum_{j \in \Z} \psi_j W_{t + \tau - j}, \mu + \sum_{j \in \Z} \psi_j W_{t - j} \right) \\
	 		&= \operatorname{Cov}\left( \sum_{i \in \Z} \psi_i W_{t + \tau - i}, \sum_{j \in \Z} \psi_j W_{t - j} \right) \\
	 		&= \sum_{i \in \Z} \sum_{j \in \Z} \psi_i \psi_j \operatorname{Cov}\left( W_{t + \tau - i}, W_{t - j} \right) \\
	 		&= \sum_{i \in \Z} \sum_{j \in \Z} \psi_{\tau + i} \psi_j \operatorname{Cov}\left( W_{t - i}, W_{t - j} \right) \\
	 		&= \sum_{i \in \Z} \sum_{j \in \Z} \psi_{\tau + i} \psi_j \operatorname{Var}\left(W_{t - i} \right) \delta^i_j \\
	 		&= \sigma^2 \sum_{i \in \Z} \psi_{\tau + i} \psi_i \qedhere
	 	\end{align*}
	 \end{proof}
	 \subsection{Estimation on stationary processes}
	 The sample mean of a sample $X_1, \ldots, X_n$ of a stationary timeseries
	 \[ \bar X := \frac{1}{n} \sum_{i=1}^n X_i \]
	 is unbiased for $\mu := \mu_X(t)$, this is always true if a sample has the same first moment. We can compute the variance:
	 \begin{align*}
	 	\operatorname{Var}(\bar X) &= \operatorname{Var}\left(\frac{1}{n} \sum_{i=1}^n X_i\right) \\
	 	&= \frac{1}{n^2} \operatorname{Var}\left(\sum_{i=1}^n X_i\right) \\
	 	&= \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n \operatorname{Cov}(X_i, X_j) \\
	 	&= \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n \gamma_X(i, j) \\
	 	&= \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n \gamma_X(j + (i - j), j) \\
	 	&= \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n \gamma_X(i - j) \\
	 	&= \frac{1}{n^2} \left(n \gamma_X(0) + 2 \sum_{i=1}^{n-1} (n - i) \gamma_X(i)\right) \\
	 \end{align*}
	 Now if the timeseries is uncorrelated, we get that $\operatorname{Var}(\bar X) = \frac{\gamma_X(0)}{n}$ which is neat, very similar to the normal estimator. We can show that
	 \[ \lim_{n \to \infty} n \operatorname{Var}(\bar X_n) = \sum_{j \in \Z} \gamma_X(j) \]
	 so that we have some sort of asymptotic behaviour, and indeed the normal variance as the sample size increases goes to 0. Now if we dealt with a \textbf{linear process} we can look at
	 \[ \sum_{\tau \in \Z} \gamma_X(\tau) = \sum_{\tau \in \Z} \sum_{j \in \Z} \sigma^2 \psi_j\psi_{\tau + j} = \sigma^2 \left(\sum_{j \in \Z} \psi_j\right)^2 \]
	 which gives rise to the theorem
	 \begin{theorem}
	 	Let $(X_t)$ be a linear process whose coefficients sum to a nonzero real number, then distributionally
	 	\[ \sqrt{n}(\bar X_n - \mu) \to N(0, V) \]
	 	where
	 	\[ V := \sigma^2 \left(\sum_{j \in \Z} \psi_j\right)^2 \]
	 \end{theorem}
 \end{document}