I recently came across a strange and beautiful connection between the Gaussian $e^{-x^2}$ and the method of least squares. It turns out that the square in $e^{-x^2}$ and the square in ``least squares'' is the same square.
Let $(x_i,y_i)$ (with $1\leq i \leq n$) be the data set, and assume that for each $x$, the $y$'s are normally distributed with mean $\mu (x-\bar x)$$\mu(x)=\alpha x+\beta$ and variance $\sigma^2$. Then, the likelihood of generating our data (assuming that the data points are independent) is $$\prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(y_i-\mu(x_i))^2}{2\sigma^2}\right) =\left(\frac{1}{\sqrt{2\pi \sigma^2}}\right)^n \exp\left( \frac{-1}{2\sigma^2} \sum_{i=1}^n (y_i - \alpha x+\beta)^2 \right)$$ We would obviously want to choose the parameters $\alpha,\beta$ so that the likelihood is maximized, and this is accomplished by minimizing $$\sum_{i=1}^n (y_i - \alpha x+\beta)^2.$$ In other words, the least squares approximation is the one that makes the data set most likely to happen.