Exactly how you count symmetries depends on your normalizations, and in particular on whether you use divided (i.e. "exponential") power series or ordinary power series. I believe strongly that divided powers are the way to go. To establish notation, I will first review the preliminaries, that I'm sure you already know.

So let me consider Dyson series for integrals of the form:
$$\int_{x\in X} \exp\left( \sum_{n\geq 2} \frac{c_n x^n}{n!} \right) {\rm d}x$$
In actual physics examples, $X$ is an infinite-dimensional space and the measure ${\rm d}x$ does not exist, but no matter. Each coefficient $c_n$ is a symmetric tensor $X^{\otimes n} \to \mathbb R$ or $\mathbb C$, and we suppose that $c_2$ is invertible as a map $X \to X^*$ — in infinite-dimensional settings, this is a very problematic supposition, and leads to questions of renormalization, which I will not address here. (Incidentally, in all actual physical quantum field theories, $c_2$ is not a priori invertible, because of gauge symmetry; so I am assuming that you have fixed that however suits your fancy.) Finally, Dyson series / Feynman diagrams are by definition *perturbative*, meaning that you need some perturbation parameter. There are various choices for how to do this; ultimately what's important is that ratios $\lvert c_n\rvert / \lvert c_2 \rvert^{n/2}$ are infinitesimal for $n \geq 3$. In your example, $c_2$ is normalized to unity, and $c_4 = \lambda \ll 1$. Another option is to set all coefficients $\lvert c_n\rvert \sim \hbar^{-1}$, where $\hbar \ll 1$.

In any case, the Dyson series is correctly calculated by expanding
$$ \int_{x\in X} \exp\left( \sum_{n\geq 2} \frac{c_n x^n}{n!} \right) {\rm d}x = \int_{x\in X} \exp \left( \frac{c_2x^2}2\right) \sum_{m=0}^\infty \left( \sum_{n\geq 3} \frac{c_n x^n}{n!} \right)^m {\rm d}x, $$
using "Wick's theorem" for $b_\ell$ a totally symmetric tensor $X^{\otimes \ell} \to \mathbb R$:
$$ \int_{x\in X} \exp \left( \frac{c_2x^2}2\right) \frac{b_\ell x^\ell}{\ell!} = \begin{cases} 0, & \ell \text{ odd} \\ \text{normalization} \times \frac{b_\ell (c_2/2)^{\ell/2}}{(\ell/2)!}, & \ell \text{ even} \end{cases},$$
(where $\text{normalization}$ involves $\det c_2$, factors of $\sqrt{2\pi}$, and doesn't make really sense in infinite dimensions), and recognizing what all of these exponential generating functions are counting.

When all the dust has settled, what these are counting is (all: disconnected, empty, etc., are allowed) graphs with trivalent-and-higher vertices, mod symmetries. If you throw a $\log$ out in front of the whole integral, then the resulting exponential generating function counts connected graphs, still mod symmetries.

In any case, to actually answer your question for the Dyson series above: the value of each graph is computed by placing a $c_n$ at each vertex of valence $n$ and a $(c_2)^{-1}$ on each edge, and contracting these tensors according to the graph. The symmetry factor is precisely the number of automorphisms of the graph, defined as follows:

**Definition:** A *Feynman graph* is a collection $H$ of *half edges* along with two partitions: one (*edges*) into sets of size precisely $2$, and the other (*vertices*) into sets of size at least $3$. (Since each graph will be weighted by the coefficients $c_n$, if you want graphs for a theory with most $c_n = 0$, you can restrict to only graphs with the prescribed valences. If $X$ splits as a direct sum $X = X_1 \oplus X_2$, you can reasonably define graphs with half-edges colored by $X_1,X_2$. If you want to compute an integral with an "observable", you should consider graphs with prescribed "external half-edges".) This is the correct definition, because it gives the correct notion of iso/automorphism. In particular, an *isomorphism of Feynman graphs* is a bijection of half edges that induces bijections on the partitions. An *automorphism* is an isomorphism from a Feynman graph to itself. Then the symmetry factor is precisely the number of automorphism in this sense.

(A different notion of "graph" more commonly used in mathematics is that of an "adjacency matrix", which is a $\mathbb Z_{\geq 0}$-valued matrix indexed by the vertices of the graph. The natural notion of isomorphism of such things is a bijection of vertices that induces an equality of adjacency matrices, but this is not the right one for counting symmetries of Feynman diagrams, as, for example, it gives "$1$" for each of the figure-eight and tadpole diagrams.)

For small graphs, the automorphism group splits as a direct product. For example, for the figure eight, the automorphism group is $(\mathbb Z/2) \ltimes (\mathbb Z/2)^2$, where the $(\mathbb Z/2)^2$ acts as flips of the two lobes of the figure eight, and the left-hand $\mathbb Z/2$ switches the two lobes. In a "$\phi^3$" theory, the theta graph has automorphism group $(\mathbb Z/2) \ltimes S_3$, where $S_3$, the symmetric group on three objects, acts to permute the edges, and the $\mathbb Z/2$ switches the two vertices. So the figure eight has $2\times 2^2 = 8$ symmetries, and the theta has $2\times 3! = 12$ symmetries. For small graphs, you can really just read off these semidirect-product decompositions: self-loops contribute factors of $\mathbb Z/2$, each collection of $k$ parallel edges contributes a factor of $S_k$, if you permute vertices, that's another factor, etc. For large graphs, the computations become much harder, and I think I can find a graph whose symmetry group is any finite group you want (but don't quote me on that thought). A slow way to do the computation is to consider all possible permutations of half edges, and see if they induce automorphism. This is slow because there are $(\text{lots})!$ many such permutations. The whole point of working with Dyson series is to avoid this. Fortunately, no one ever computes large graphs, because even just contracting all the tensors is so hard for those.

Finally, I cannot help but mention one more fact you probably know. These Dyson series almost never have positive radius of convergence in the perturbation parameters. In particular, even after dividing by symmetry factors, there are still a lot of graphs, and so the coefficients grow as $n!$. A good exercise is to compute the Dyson series for the conditionally-convergent Riemann integral $\int_{\mathbb R} \exp \frac i \hbar \bigl( \frac{x^2}2 + \lambda \frac{x^3}6 \bigr){\rm d}x$. (The $i$ is there to make the integral conditionally converge; I tend to work with $\hbar$ as my perturbation parameter.) Then the Dyson series grows as $\sum \frac{(6n!)}{(2n)!\,(3n)!} \alpha^n$, where $\alpha$ is linear in $\lambda\hbar$ (something like $\alpha = \lambda\hbar / 6$). The point is that by Stirling's formula $\frac{(6n!)}{(2n)!\,(3n)!} \approx \beta^n n!$ for some positive $\beta$. In any case, this is not at all surprising. Integrals of the form $\int \exp \sum c_nx^n/n!$ are essentially never analytic at $0,\infty$ in the coefficients — for example, Gauss's formula $\int_{\mathbb R}\exp(-a^{-1}x^2/2){\rm d}x = \sqrt{2\pi a}$ has a ramified point at $a = 0,\infty$, and so the Riemann integral extends from its domain of convergence (${\rm Re}(a) > 0$) to a double-valued complex function in $a$.