Anytime I see an $n!$ in some formula, my instinct is to look for the symmetric group on $n$ letters coming in somewhere. I have never done this seriously with the $n!$ in Taylor's theorem.

Question: Is there some way to see the $n!$ in Taylor's theorem coming naturally from a symmetry group?

Possible lead:

Here is a definition of $f^{(n)}(a)$ which does not depend on finding earlier derivatives: Let $g: \mathbb{R}^n \rightarrow \mathbb{R}$ be defined by $g(x_1,x_2,x_3, ..., x_n)$ is the lead coefficient of the unique $n^{th}$ degree polynomial passing through $(a, f(a)), (x_1, f(x_1)),(x_2, f(x_2)),...,(x_n, f(x_n))$. Then $f^{(n)}(a)$ is $1/n!$ times the limit of $g(x_1,x_2,x_3, ..., x_n)$ as $(x_1,x_2,x_3, ..., x_n)$ approaches $(a,a,a,...,a)$. Could the $1/n!$ be related to the symmetry of $g$ under exchange of coordinates?

There must be many ways to think of this. Here's one: 

The symmetric group is involved with homogeneous polynomials of degree $n$ because they correspond to symmetric multilinear functions of $n$ variables, and division by $n$ factorial appears when recovering the former from the latter.
For example, a homogeneous polynomial map $f:V\to W$ of degree $3$ determines a map $F:V\times V\times V\to W$ by
$$
F(x,y,z)=f(x+y+z)-f(x+y)-f(x+z)-f(y+z)
+f(x)+f(y)+f(z)-f(0),
$$
and $f(x)$ can be recovered as 
$$\frac{F(x,x,x)}{6}.
$$

The same formula for $F(x,y,z)$ can be applied to other functions $f(x)$. Applied to a polynomial of degree $<3$ it gives zero. Applied to a polynomial of degree $\le3$ it gives the multilinear function that corresponds to the purely cubic part of $f$. Applied to any function at all it gives a symmetric function of $(x,y,z)$ that vanishes when $x$ or $y$ or $z$ is zero.

The foundation of all differential calculus is the process that takes $f$, subtracts $f(0)$, and then makes the best linear approximation of the result. If you apply this process in all three variables to $f(x+y+z)$, you are first making $F(x,y,z)$ and then (assuming that $f$ was sufficiently smooth) making best linear approximation in all variables. If you then set $x=y=z$ and divide by $3$ factorial, you get the third term in the Taylor series of $f$.

Of course, if you perform the linearizations in the three variables one after another then you see this trilinear map as a derivative of a  derivative of a derivative. 

The reason I look at it this way is that a homotopical categorical analogue of this plays a big role in my "functor calculus". The analogue there of division by $n$ factorial is a homotopy orbit spectrum for an action of the symmetric group. 

In a little more detail: 

Let $V$ and $W$ be model categories in which filtered homotopy colimits commute with finite homotopy limits. Assume that $W$ is stable, meaning that the final object is equivalent to initial object and that homotopy pushout squares are the same as homotopy pullback squares. For example, $W$ might be the category of spectra and $V$ might be spaces or spectra. Consider functors $f:V\to W$ that are homotopy-invariant (preserve weak equivalences). Say that $f$ has degree $\le 1$ if it preserves homotopy pushout squares. Say that it has degree $\le d$ if it takes those $(d+1)$-dimensional square diagrams all of whose $2$-dimensional faces are homotopy pushouts to homotopy pushout cubes, those in which the last object is the homotopy colimit of the others. Call $f$ reduced if it takes the trivial object $\star$ to itself (up to equivalence). Call it linear if it has degree $\le 1$ and is reduced. There is a linearization process that takes a reduced functor and makes the universal linear functor under it: basically $hocolim \Omega^nf\Sigma^n$. This can be easily generalized to make something called $P_1$ (first polynomial approximation) which takes any functor $f$, reduced or not, and makes the universal degree $\le 1$ functor $P_1f$ under it. (If we first reduce $f$ by taking the homotopy fiber of $f\to f(\star)$ and then linearize, it's the same as first doing $P_1$ and then reducing.) It can be further generalized to make something called $P_d$ ($d$th polynomial approximation) which takes any functor $f$ and makes the universal degree $\le d$ functor $P_df$ under it. The homotopy fiber of the canonical map $P_df\to P_{d-1}f$ is always a degree $\le d$ functor such that $P_{d-1}$ of it is trivial. Call such functors homogeneous of degree $d$. When $d=1$ this is the same as linear. 

There is the functor of $d$ variables $(X_1,\dots,X_d)\mapsto f(X_1+\dots +X_d)$, where $+$ means (derived) coproduct. This is symmetric in the sense that there are isomorphisms when you permute the inputs, satisfying the obvious identities. Reduce it in all variables simultaneously, by first making a cubical diagram consisting of the objects $f(X_S)$, coproduct of all the $X_i$ for $i\in S$, and then looking at the homotopy fiber of the canonical  map from $f(X_1+\dots +X_d)$ to the holim of all the rest. I call this reduced symmetric functor of $d$ variables the $d$th crosseffect of $f$. If we simultaneously linearize it in all variables, we get a symmetric multilinear functor.

If $f$ has degree $\le (d-1)$ then its $d$th crosseffect is contractible. If $f$ is homogeneous of degree $d$ then its $d$th crosseffect is already multilinear, and this way of making a symmetric multilinear functor from a homogeneous degree $d$ functor can be inverted: $f(X)$ is the homotopy orbit object for the action of the $d$th symmetric group on $F(X,\dots ,X)$, where the action on $F(X,\dots ,X)$ is the one that you get from the symmetry on $F(X_1,\dots ,X_d)$.

For general $f$ the multilinearization of the crosseffect is the same symmetric multilinear functor that corresponds in this way to the $d$th homogeneous part of $f$, i.e. to the homotopy fiber of $P_df\to P_{d-1}f$.

Again, by linearizing in one variable at a time you can work out a sense in which the multilinear functor corresponding to the $d$th homogeneous part is in fact an iterated derivative. But in this context it is important to have the option of performing linearization in all variables simultaneously, so to speak, in order that the result should be symmetric in the sense that is required for reconstructing the corresponding homogeneous functor. (Here symmetry is a structure, not a property.)
