Let X_i be individually and independently distributed by a probabilistic distribution with a mean of \mu and a variance of \sigma^2.
Note that the distribution is not necessarily a normal distribution.
The central limit theorem states that as the number of random variables n \rightarrow \infty, the mean of those random variables \bar{X} = \frac{1}{N} \sum_{n=1}^{\infty} X_i approaches a normal distribution.
\begin{align*}
\bar{X} \sim \mathcal{N}(\bar{X}|\mu, \frac{\sigma^2}{n})
\end{align*}
The following equations should be noted here,
\begin{align*}
f_{\bar{X}}(\bar{x}) &= f_{\frac{X_1}{n} + \cdots + \frac{X_n}{n}}(\frac{x_1}{n} + \cdots + \frac{x_n}{n}) \\
&= f_{\frac{X_1}{n}}(\frac{x_1}{n}) \cdots f_{\frac{X_n}{n}}(\frac{x_n}{n}) \\
&= f_{X_1}(x_1) \cdots f_{X_n}(x_n) \\
\end{align*}
Let's think about the moment generating function of \bar{X}.
\begin{align*}
M_{\bar{X}}(\theta) &= E[e^{\bar{X} \theta}] \\
&= \int e^{\bar{x} \theta} f_{\bar{X}}(\bar{x}) d\bar{x} \\
&= \int \cdots \int e^{\frac{\theta}{n} (x_1 + \cdots + x_n)} f_{X_1}(x_1) \cdots f_{X_n}(x_n) dx_1 \cdots dx_n \\
&= \int e^{\frac{\theta}{n} x_1} f_{X_1}(x_1) dx_1 \cdots \int e^{\frac{\theta}{n} x_n} f_{X_n}(x_n) dx_n \\
&= M_{X}(\frac{\theta}{n}) \cdots M_{X}(\frac{\theta}{n}) \\
&= \left\{ M_{X}(\frac{\theta}{n}) \right\}^n \\
\end{align*}
Here we will use Maclaurin series of e^t and \ln(1 + t), let's remember
\begin{align*}
e^t = 1 + \frac{1}{1!} t + \frac{1}{2!} t^2 + \cdots \\
\ln(1 + t) = t - \frac{1}{2} t^2 + \frac{1}{3} t^3 \cdots \\
\end{align*}
\begin{align*}
\ln M_{\bar{X}}(\theta) &= n \ln M_{X}(\frac{\theta}{n}) \\
&= n \ln E[e^{\frac{\theta}{n} X}] \\
&= n \ln E\left[
1 + \frac{1}{1!} \left( \frac{\theta}{n} \right)X + \frac{1}{2!} \left( \frac{\theta}{n} \right)^2 X^2 + \cdots
\right]\\
&= n \ln \left\{
1 + \frac{1}{1!} \left( \frac{\theta}{n} \right)E[X] + \frac{1}{2!} \left( \frac{\theta}{n} \right)^2 E[X^2] + \cdots
\right\}\\
&= n \left\{
\left(
\frac{1}{1!} \left( \frac{\theta}{n} \right)E[X] + \frac{1}{2!} \left( \frac{\theta}{n} \right)^2 E[X^2] + \cdots
\right)
- \frac{1}{2} \left(
\frac{1}{1!} \left( \frac{\theta}{n} \right)E[X] + \frac{1}{2!} \left( \frac{\theta}{n} \right)^2 E[X^2] + \cdots
\right)^2
+ \cdots
\right\}\\
&\overset{n \rightarrow \infty}{=} \theta E[X] \\
&= \mu \theta
\end{align*}
Let's think about the moment generating function M_Y(\theta) of a random variable Y which distributes over \mathcal{N}(X|\mu, \bar{\sigma}^2) where \bar{\sigma}^2 = \frac{\sigma^2}{n}.
(see: this link for the derivation of the equations below.)
\begin{align*}
\ln M_Y(\theta) &= E[e^{\theta Y}] \\
&= \mu \theta + \frac{\bar{\sigma}^2 \theta^2}{2} \\
&= \mu \theta + \frac{\sigma^2 \theta^2}{2 n} \\
&\overset{n \rightarrow \infty}{=} \mu \theta
\end{align*}
Therefore,
\begin{align*}
M_{\bar{X}}(\theta) = M_Y(\theta)
\end{align*}
this indicates \bar{X} distiributes as \mathcal{N}(\mu, \frac{\sigma^2}{n}) .
Discussion