Why doesn't a simple mean give the position of a centroid in a polygon?

The centroid of a polygon is indeed its center of mass -- but the mass of a polygon is uniformly distributed over its surface, not only at the vertices. You're right that if the mass were split evenly among the vertices only, the centroid would be the arithmetic mean of the coordinates of the vertices.

It just so happens that both definitions are equivalent (mass evenly distributed over the surface vs mass at the vertices only) for simple shapes like triangles and rectangles.

Laurent is right, a polygon is not just the vertices, but the whole region. The arithmetic mean of the vertices would give the centroid if the (equal) masses were concentrated at the vertices. While that answers your question, perhaps for the future, the below might be useful.

In case of a mass distributed over a region (not just a polygon), Green's Theorem might be helpful in calculating the area and the centroid:

In your case,

If $\displaystyle D$ is the region of the polygon, then the x-coordinate of the center of mass is given by the area of the polygon times

$\displaystyle \iint_{D} x dxdy$ which by Green's theorem is same as

$\displaystyle \oint_{C} \frac{x^2}{2} dy$, where the line integral is taken over the perimeter of the polygon.

For centroid we choose $\displaystyle M(x,y) = \frac{x^2}{2}$ and $\displaystyle L(x,y) = 0$

The area of the polygon is given by

$\displaystyle \iint_{D} 1 dxdy$ and can be written as a line integral.

For area we choose $\displaystyle M(x,y) = \frac{x}{2}$ and $\displaystyle L(x,y) = \frac{-y}{2}$

Yes, yes, but why? The current answers (as well as Wikipedia) do not contain enough detail to understand these formulas immediately. So let's start with the very basic definition of the area centroid: $$ \vec{C} = \frac{\iint \vec{r}(x,y)\,dx\,dy}{\iint dx\,dy} = \frac{(\iint x\,dx\,dy,\iint y\,dx\,dy)}{\iint dx\,dy} = \frac{(m_x,m_y)}{A}=(C_x,C_y) $$ Allright, let's get rid of the double integrals in the first place, by employing Green's theorem : $$ \iint \left( \frac{\partial M}{\partial x} - \frac{\partial L}{\partial y} \right) dx\,dy = \oint \left( L\,dx + M\,dy \right) $$ At the edges of the (convex) polygon we have: $$\begin{cases} x = x_i + (x_{i+}-x_i)\,t \\ y = y_i + (y_{i+}-y_i)\,t \end{cases} \quad \mbox{with} \quad \begin{cases} i = 0,1,2,\cdots,n-1 \\ i+=i+1\mod n \end{cases} \quad \mbox{and} \quad 0 \le t \le 1 $$ Then by substitution of $M(x,y) = x$ and $L(x,y) = 0$ we have: $$ A = \iint dx\,dy = \oint x\,dy = \sum_{i=0}^{n-1} \int_0^1 \left[x_i + (x_{i+}-x_i)\,t\right](y_{i+}-y_i)\,dt=\\ \sum_{i=0}^{n-1}(y_{i+}-y_i)\left[x_i\left.t\right|_0^1 + (x_{i+}-x_i)\frac{1}{2}\left.t^2\right|_0^1\right] = \frac{1}{2}\sum_{i=0}^{n-1}(x_{i+}+x_i)(y_{i+}-y_i)=\\ \frac{1}{2}\sum_{i=0}^{n-1} (x_iy_{i+}-x_{i+}y_i) $$ The last move by telescoping.
The main integral for the $x$-coordinate of the centroid is,
with $M(x,y) = x^2/2$ and $L(x,y) = 0$: $$ m_x = \iint x\,dx\,dy = \oint \frac{1}{2}x^2 \,dy = \frac{1}{2}\sum_{i=0}^{n-1}(y_{i+}-y_i)\int_0^1\left[x_i + (x_{i+}-x_i)\,t\right]^2\,dt=\\ \frac{1}{2}\sum_{i=0}^{n-1}(y_{i+}-y_i)\left[x_i^2\left.t\right|_0^1+2x_i(x_{i+}-x_i)\frac{1}{2}\left.t^2\right|_0^1 +(x_{i+}-x_i)^2\frac{1}{3}\left.t^3\right|_0^1\right]=\\ \frac{1}{2}\sum_{i=0}^{n-1}(y_{i+}-y_i)\left[x_i^2+x_i(x_{i+}-x_i)+\frac{1}{3}(x_{i+}-x_i)^2\right]=\\ \frac{1}{6}\sum_{i=0}^{n-1}(y_{i+}-y_i)\left[x_{i+}^2+x_ix_{i+}+x_i^2\right]=\\ \frac{1}{6}\sum_{i=0}^{n-1}\left[x_ix_{i+}y_{i+}+x_i^2y_{i+}-x_{i+}^2y_i-x_ix_{i+}y_i\right]\quad\Longrightarrow\\ m_x = \frac{1}{6}\sum_{i=0}^{n-1}(x_i+x_{i+})(x_iy_{i+}-x_{i+}y_i) $$ The last two moves after telescoping again.
The main integral for the $y$-coordinate of the area centroid is,
with $M(x,y) = 0$ and $L(x,y) = -y^2/2$: $$ m_y = \iint y\,dx\,dy = \oint -\frac{1}{2}y^2 \,dx = -\frac{1}{2}\sum_{i=0}^{n-1}(x_{i+}-x_i)\int_0^1\left[y_i + (y_{i+}-y_i)\,t\right]^2\,dt $$ Which is similar to the main integral for the $x$-coordinate of the centroid: $$ m_x = \iint x\,dx\,dy = \oint \frac{1}{2}x^2 \,dy = \frac{1}{2}\sum_{i=0}^{n-1}(y_{i+}-y_i)\int_0^1\left[x_i + (x_{i+}-x_i)\,t\right]^2\,dt $$ It is seen that everything is the same if we just exchange $x$ and $y$, except for the minus sign, hence: $$ m_y = -\frac{1}{6}\sum_{i=0}^{n-1}(y_i+y_{i+})(y_ix_{i+}-y_{i+}x_i)=\frac{1}{6}\sum_{i=0}^{n-1}(y_i+y_{i+})(x_iy_{i+}-x_{i+}y_i) $$ Combining the partial results found gives the end result, as displayed in the question.