What is the joint distribution of sample mean and sample variance of normal distribution? $X_i \sim N( \mu,\sigma^2)$, define $\overline X =\dfrac{1}{n} \sum\limits_{i = 1}^n X_i $, $S^2 = \dfrac{1}{n - 1}\sum\limits_{n = 1}^n \left( {X_i - \overline X} \right)^2$. What is the distribution of 
$$
\sqrt n \left( \begin{array}{c}
\overline X  - \mu \\
S^2 - \sigma ^2
\end{array} \right)
$$

This is a proof that $\overline X$ and $S^2$ are independent of each other, as mentioned in Michael Hardy's comment.
Construct $y_1=\bar x=\frac{1}{n}(x_1+\cdots+x_n), y_2=x_2-\bar x, \ldots, y_n=x_n-\bar x$. Adding them all up, we can solve for $x$'s in terms of $y$'s
$$\begin{align} 
x_1&=y_1-y_2-\cdots-y_n\\
x_2&=y_1+y_2\\
\vdots\\
x_n&=y_1+y_n
\end{align}$$
So the Jacobian of transformation is
$$J_n=\left| {\begin{array}{*{20}{c}}
1&-1&-1&\cdots&-1\\
1&1&0&\cdots&0\\
1&0&1&\cdots&0\\
\vdots&\vdots&\vdots&\vdots&\vdots\\
1&0&0&\cdots&1
\end{array}} \right|. $$
Expanding the determinant using the elements of the bottom line, we have $J_n=1+J_{n-1}$, where the $1$ on the right side is obtained again by expanding recursively using the bottom line. Since $J_2=2$, we have immediately $J_n=n$.
Thus, we can write the transformation formula for joint pdf of $Y_1,Y_2,\ldots,Y_n$ from that of $X_1,X_2,\ldots,X_n$:
$$\begin{align} f_{Y_1Y_2\cdots Y_n}(y_1,y_2,\ldots,y_n)&=f_{X_1X_2\cdots X_n}(x_1,x_2,\ldots,x_n)|J_n|\\
&=n\left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n\exp\left\{-\frac{1}{2\sigma^2}\sum\limits_{i = 1}^n(x_i-\mu)^2\right\}.
\end{align}$$
To proceed, observe that $$\begin{align} 
\sum\limits_{i=1}^n(x_i-\mu)^2&=\sum\limits_{i=1}^n(x_i-\bar x)^2+n(\bar x-\mu)^2\\
&=(-y_2-\cdots-y_n)^2+\sum\limits_{i=2}^n y_i^2+n(y_1-\mu)^2,
\end{align}$$
where the second equation is based on the solution of $x$'s in terms of $y$'s.
Plugging the sum of squares back to the joint pdf, we get
$$\begin{align}
f_{Y_1Y_2\cdots Y_n}(y_1,y_2,\ldots,y_n)&=n\left(\frac{1}{\sqrt{2\pi}\sigma}\right)^n\exp\left\{-\frac{1}{2\sigma^2}\biggl[\bigl(\sum\limits_{i=2}^n y_i\bigr)^2+\sum\limits_{i=2}^n y_i^2+n(y_1-\mu)^2\biggr]\right\}\\
&=f(y_1)g(y_2,\cdots,y_n).
\end{align}$$
Since the final joint pdf is a product of a function of only $y_1$ and a function of only $y_2,\cdots,y_n$, we have $Y_1$ is independent of $(Y_2,\cdots,Y_n)$. Since $Y_1$ is just $\overline X$ and $S^2$ is a function of $Y_2,\cdots, Y_n$ (see the expansion of sum of squares above), we established the fact that sample mean is independent of sample variance.
