I'm just going to add a few more details to the paper that Stefan linked to.
First of all, the matrix $H_1$ displayed in the paper is incorrect; every element should be multiplied by a factor of $\frac{1}{8}$.
$$ H_1 = \frac{1}{8} \begin{bmatrix}
4 & 4 & 0 & 0 \\
1 & 6 & 1 & 0 \\
0 & 4 & 4 & 0 \\
0 & 1 & 6 & 1
\end{bmatrix}$$
Second, the derivation of the new face point formula for quads in the first half of the paper is confusing, so I will fill in some of the details.
The matrix $G_1 = HGH^T$ contains the 16 control points for one sub-quad of the original quad with control points $G$. The entries of $G_1$ are referred to with $q_{i,j}$. Of most interest are the elements $q_{1,1}$, $q_{1,2}$ and $q_{2,2}$, which represent control points for a corner, edge, and interior of the $4\times 4$ control point grid for the sub-quad. All the other elements have similar formulas due to the symmetry of the subdivision process.
By carrying out the matrix multiplication, we get:
$q_{2,2} = \frac{P_{1,1} +P_{3,1} +P_{3,1} +P_{3,3}}{64} + \frac{6(P_{1,2} +P_{2,1} +P_{3,2} +P_{2,3})}{64} + \frac{36P_{2,2}}{64}$
We can see from this formula that, when relocating a control point, its new location is determined by the $3\times 3$ grid of points surrounding it. Here we see that in this grid, corner points have $\frac{1}{64}$ influence, edge points have $\frac{6}{64} = \frac{3}{32}$ influence, and the center ("old") point has $\frac{36}{64} = \frac{9}{16}$ influence. These weights are consistent with the ones given at the bottom of this blog post.
With a bit of rearranging, we get:
$q_{2,2} = \frac{1}{4} \big[ \frac{1}{4}(\frac{P_{1,1} + P_{1,2} + P_{2,1} + P_{2,2}}{4} + \frac{P_{1,2} + P_{2,2} + P_{1,3} + P_{2,3}}{4} + \frac{P_{2,1} + P_{2,2} + P_{3,2} + P_{3,1}}{4} + \frac{P_{2,2} + P_{2,3} + P_{3,2} + P_{3,3}}{4}) + 2\frac{1}{4}(\frac{P_{1,2} + P_{2,2}}{2} + \frac{P_{2,1} + P_{2,2}}{2} + \frac{P_{2,3} + P_{2,2}}{2} + \frac{P_{3,2} + P_{2,2}}{2}) + P_{2,2}\big]$
The first term is merely the average of the four surrounding face points; in other words, it is $F$. The second term is the average of the four surrounding edge midpoints multiplied by $2$; in other words, it is $2R$. Finally, $P_{2,2}$ is simply the original control point $P$.
And so we see this is the same as the original formula for quads: $\frac{1}{4}(F + 2R + P)$.
The paper does not provide a derivation of the general case for an $n$-sided polygon, but it seems like a reasonable extension based on intuition.
As an aside, here is an explanation for the other points, and their meanings as seen on wikipedia:
$q_{1,1} = \frac{P_{1,1} +P_{1,2} +P_{2,1} +P_{2,2}}{4}$
This is simply the formula for adding a new control point at the center of an existing quad.
$q_{1,2} = \frac{P_{1,1} +P_{1,2} +P_{2,1} +P_{2,3}}{16} + \frac{6(P_{1,2} +P_{2,1})}{16}$
Rearranging this, we see it is just the average of the two neighbouring face points and the edge endpoints:
$q_{1,2} = \frac{1}{4}(\frac{P_{1,1} +P_{2,1} +P_{1,2} +P_{2,2}}{4} + \frac{P_{1,2} +P_{2,2} +P_{1,3} +P_{2,3}}{4} + P_{1,2} + P_{2,1})$
