Covariance matrix - wrong results either in Matlab or in formula? The definition of a covariance for discrete variables (taken from https://en.wikipedia.org/wiki/Covariance) is:
$$\operatorname{cov} (X,Y)=\frac{1}{n}\sum_\limits{i=1}^n (x_i-E(X))(y_i-E(Y))$$
Therefore, for matrix $A = \begin{bmatrix}1 & 2 & 3\\2 & 4 & 6\end{bmatrix}$ I would expect $\operatorname{cov}(A) = \begin{bmatrix}0.5 & 0.5 & 0.75\\ 0.5 & 2.0 & 1.5 \\ 0.75 & 1.5 & 4.5\end{bmatrix}$ because:
$E\left(\begin{bmatrix}1\\2\end{bmatrix}\right) = \frac{1+2}{2} = 1.5$
$E\left(\begin{bmatrix}2\\4\end{bmatrix}\right) = \frac{2+4}{2} = 3$
Using the above mentioned formula for covariance of the first and second column I obtain $0.5$:
$\operatorname{cov}\left(\begin{bmatrix}1\\2\end{bmatrix},\begin{bmatrix}2\\4\end{bmatrix}\right) = \frac{1}{n=2}[(1-1.5)(2-3)+(2-1.5)(4-3)] = \frac{1}{2}(0.5 + 0.5) = 0.5$
Matlab however gives
> A = [1 2 3; 2 4 6];
> cov(A)

ans =

    0.5000    1.0000    1.5000
    1.0000    2.0000    3.0000
    1.5000    3.0000    4.5000

The variances in diagonal entries are correct, but for all the covariances they don't match (for the above example I get $1.0$ instead of $0.5$). What am I missing...?
 A: Matlab is using Bessel's correction.

In statistics, Bessel's correction, named after Friedrich Bessel, is the use of $n - 1$ instead of $n$ in the formula for the sample variance and sample standard deviation, where $n$ is the number of observations in a sample.

(And the same applies to covariances.)
Many statistics texts, including nearly all texts for students who don't want to understand the mathematical theory but who want to know how to use its results, state that the variance is defined as
$$
\frac 1 {n-1} \sum_{i=1}^n (x_i - \bar x)^2 \quad\text{where } \bar x = \frac 1 n \sum_{i=1}^n x_i. \tag 1
$$
However, in probability theory, the variance of a discrete random variable having $n$ possible values $x_1,\ldots, x_n$ is
\begin{align}
\sum_{i=1}^n p_i (x_i - \bar x)^2 & \quad \text{where } p_i \text{ is the probability assigned to }x_i \\
& \quad \text{and } \bar x = \sum_{i=1}^n  p_i x_i.
\end{align}
When $p_i = 1/n$ for every $i$, then this leads to the variance
$$
\frac 1 n \sum_{i=1}^n (x_i - \bar x)^2. \tag 2
$$
The reason for the seeming discrepancy is that $(1)$ and $(2)$ are two different things having two different purposes.  The purpose of $(1)$ is to estimate the variance of a population when $x_1,\ldots,x_n$ is a random sample from that population.  If the terms were $(x_i-\mu)^2$, where $\mu$ is the mean of the whole population, then one would divide by $n$, but if the terms are $(x_i-\bar x)^2$, where $\bar x$ is the mean only of the sample $x_1,\ldots,x_n$ rather than of the whole population rather than only of the small sample $x_1,\ldots,x_n$, then on average one would underestimate the population variance.
Using the "unbiased" estimate (which gives the right value on average) with $n-1$ rather than $n$ in the denominator, is of debatable merit.  For example, the mean squared error of estimation is actually smaller with the biased estimate that uses $n$ rather than $n-1$ in the denominator, and there are yet other criticisms.  But that's what's going on here.
