Proofs for statistics and normal distributions I am studying statistics, and it it, we are given many different results about what kind of estimations we can make and what kind of distributions these estimations have.
For example, assume $Y_1,\ldots,Y_n$ er iid normally distributed with some mean and variance.... then our estimation of the mean is the average, and this average is the realisation of a random variable which itself is normally distributed which mean $ = \mu$ and variance $= \sigma^2/n$. 
But our book has no proof of this, and it's not part of the course, yet I'd like to know it anyways, but I have some trouble finding proofs of this on the internet. Does any of you know of certain sites where the proofs are given?
Every proof regarding a first year course in statistics, and proof for the theory of normal distributions especially, will work. 
 A: If I understand you question correctly this Wikipedia page will give a better explanation than what I can put into the post. (ie the sum of normally distributed random variables)
If you want to know where they come from, or at least there importance, I'd recommend getting your head around central limit theorem. 
A: Do you know about moment generating functions? If so you can use this proof
https://onlinecourses.science.psu.edu/stat414/node/173
A: Because your supposition is that the $Y_i$s are normal, the Central Limit Theorem is not relevant:  the distribution of their sum is exactly normal, not merely asymptotically normal.
In fact, all you need to do is show that, for two independent normal random variables $X \sim \operatorname{Normal}(\mu_x, \sigma_x^2)$ and $Y \sim \operatorname{Normal}(\mu_y, \sigma_y^2)$, their sum is $$W = X+Y \sim \operatorname{Normal}(\mu_w = \mu_x + \mu_y, \sigma_w^2 = \sigma_x^2 + \sigma_y^2);$$ that is, the means and variances add.  This is easily accomplished through the use of generating functions or convolution.  The first method, using the MGF, follows from the relationship $$M_W(t) = \operatorname{E}[e^{Wt}] = \operatorname{E}[e^{(X+Y)t}] = \operatorname{E}[e^{Xt}e^{Yt}] \overset{\text{ind}}{=} \operatorname{E}[e^{Xt}]\operatorname{E}[e^{Yt}] = M_X(t)M_Y(t).$$  The MGF for a normal distribution is $$\begin{align*} M_X(t) &= \int_{x=-\infty}^\infty \frac{e^{tx}}{\sqrt{2\pi}\sigma_x} e^{-(x-\mu_x)^2/(2\sigma_x^2)} \, dx \\ &= \frac{1}{\sqrt{2\pi}\sigma_x}\int_{x=-\infty}^\infty \exp\left( -\frac{(x-(\mu_x +\sigma_x^2 t))^2}{2\sigma_x^2} + \frac{(\mu_x + \sigma_x^2 t)^2 - \mu_x^2}{2\sigma_x^2}\right) \, dx \\ &= \exp\left(\frac{(\mu_x + \sigma_x^2 t)^2 - \mu_x^2}{2\sigma_x^2} \right) = \exp\left(\mu_x t + \frac{\sigma_x^2 t^2}{2} \right), \end{align*} $$ where the last line comes from separating out the exponential factor that is independent of $x$, and noting that the remaining exponential factor is simply the integrand of a normal distribution with mean $\mu_x + \sigma_x^2 t$ and variance $\sigma_x^2$, hence equals $1$.  Then we immediately find $$M_W(t) = \exp\left(\mu_x t + \frac{\sigma_x^2 t^2}{2} \right) \exp\left(\mu_y t + \frac{\sigma_y^2 t^2}{2} \right) = \exp\left((\mu_x + \mu_y) t + \frac{(\sigma_x^2 + \sigma_y^2) t^2}{2} \right)$$ from which it follows that $W$ is normal with the claimed parameters.
The consequence of this above case is easily generalized via induction to the case where we have $$W = \sum_{i=1}^n X_i, \quad X_i \sim \operatorname{Normal}(\mu_i, \sigma_i^2),$$ with independent (but not necessarily identically distributed) $X_i$s, of which your case corresponds to the choice $\mu_i = \mu_j$, $\sigma_i = \sigma_j$--that is, the identically distributed case.

Addendum.  It seems you are asking about the claim that if $Y_1, \ldots, Y_n$ are IID realizations of some normal distribution with unknown mean $\mu$, then the estimator $$\hat \mu = \bar Y$$ is also normal with mean $\mu$.  In fact, this is adequately addressed in my answer above, since $$\bar Y = \frac{1}{n}\sum_{i=1}^n Y_i = \frac{W}{n}$$ using the notation I used above.  Thus we can easily see that for $\mu_1 = \mu_2 = \ldots = \mu_n$, $W/n \sim \operatorname{Normal}((n \mu)/n, (n \sigma^2)/n^2)$ from the simple scaling transformation of the random variable $W$.
A: $\newcommand{\var}{\operatorname{var}}\newcommand{\E}{\operatorname{E}}$If $Y_1,\ldots,Y_n$ are uncorrelated (and if they're independent, then they're uncorrelated) and each has expected value $\mu$ and variance $\sigma^2$, then
\begin{align}
\var\left( \frac{Y_1+\cdots+Y_n} n \right) & = \underbrace{\frac 1 {n^2}\var(Y_1+\cdots+Y_n) = \frac 1 {n^2} \left( \var(Y_1)+\cdots+\var(Y_n) \right)}_{\text{These are equal because of uncorrelatedness.}} \\[10pt]
= {} &  \frac 1 {n^2} (\sigma^2+\cdots+\sigma^2) = \frac 1 {n^2} n\sigma^2 = \frac{\sigma^2} n
\end{align}
and
$$
\E\left( \frac{Y_1+\cdots+Y_n} n \right) = \frac 1 n \E(Y_1+\cdots+Y_n) = \frac 1 n (\mu+\cdots+\mu) = \mu.
$$
None of the above relies an any assumption about the distribution of $Y_1,\ldots,Y_n$ other than that they all have variance $\sigma^2$ and expected value $\mu$ and they are uncorrelated.  In particular, it is not assumed that they are normally distributed.
If in addition they are normally distributed, then the sample mean is also normally distributed.  That can be shown by computing a convolution or by finding a characteristic function.  The latter method requires showing that if two distributions have the same characteristic function then they are the same distribution.
If $Y_1,\ldots,Y_n$ are not normally distributed by all of the assumptions above hold, and in addition they are independent and identically distributed, then the sample mean will still be nearly normally distributed if $n$ is large enough.  How large is large enough depends on the distribution of each $Y_k$.  For a very skewed distribution you'd need a bigger $n$.  Google the term "central  limit theorem".
