Can a convergent perturbation series converge to the wrong answer for this type of integral? Consider a convergent integral depending on a parameter $g$ of the form $$I(g) = \int_{-\infty}^{\infty} dx e^{-f(x)+g h(x)}$$ and suppose that the series
$$S(g) = \sum_{n=0}^\infty \frac{g^n}{n!} \int_{-\infty}^{\infty}dx h(x)^n e^{-f(x)}$$ converges.
Does it necessarily follow that $I(g) = S(g)$ in the range of $g$ for which we have both $I(g)$ finite and the series $S(g)$ convergent?
In particular, if the answer is no, please give an example of an $f(x)$ and $g(x)$ where it does not hold. I have in mind $f,g,h$ real and will assume so below, but I think I'll be happy with answers that don't assume that.

Here are my thoughts. My mindset is to try to find a counterexample and see what obstacles there are to a counterexample.
First, if we try to find a counterexample (i.e. if we try to find an $I$ and $S$ that are not equal), we must be wary of the Fubini theorem. See for example this answer on swapping integrals and sums; if the integral of the sum of the absolute values of the terms is finite, then one can swap the integral and the sum.
If I understand correctly, applying this concept to this problem, $I(g)$ and $S(g)$ will be guaranteed equal should the quantity
$$J(g) = \int_{-\infty}^{\infty} dx e^{-f(x)+|g h(x)|}$$ be finite.
Thus, if we are trying to find an $I(g)$ and $S(g)$ that are not equal in order to have a nifty counterexample, we want the above $J(g)$ to diverge.
Some of the simplest ideas, like $\int_{-\infty}^{\infty} dx e^{-x^2 + (-gx^4)}$, sadly lead to series $S(g)$ with a radius of convergence of zero, so they can't be used as nifty counterexamples, since they trivially agree with $I(g)$ at $g=0$, the only point in the relevant range of $g$.
 A: There appears to be a split between the case where $f$ and $h$ are real and when they can be complex. In particular, $S(g)$ and $I(g)$ will agree where they are well-defined for the real case, but there are counterexamples in the complex case.

I will use this answer on the moments of the absolute value of $h(x)$ to answer my own question for the real case: when $S(g)$ has a nonzero radius of convergence $R$, then for $|g| < R$ and $f(x)$ and $h(x)$ both real, then $S(g)$ and $I(g)$ match.

First we need a quick lemma. Note that by Cauchy-Schwarz $$b(n+1)^2 \equiv \left(\int_{-\infty}^{+\infty} |h(x)|^{n+1}e^{-f(x)}dx\right)^2 \leq \int_{-\infty}^{+\infty} |h(x)|^{n}e^{-f(x)}dx\cdot\int_{-\infty}^{+\infty} |h(x)|^{n+2}e^{-f(x)}dx$$
so $$b(n+1)^2 \leq b(n) b(n+2) \leq \max(b(n), b(n+2))^2,$$
so $$b(n+1) \leq \max(b(n), b(n+2))  \tag{1}\label{eq1}$$

Now consider $S(g)$. When $|g|<R$, $S(g)$ absolutely converges by the ratio test, so in this range,
$$\sum_{n=0}^\infty \frac{|g|^n}{n!} \left|\int_{-\infty}^{\infty}dx h(x)^n e^{-f(x)}\right| < \infty $$
It follows then that summing only the even terms
$$\sum_{n=0}^\infty \frac{|g|^{2n}}{2n!} \left|\int_{-\infty}^{\infty}dx h(x)^{2n} e^{-f(x)}\right| =  \sum_{n=0}^\infty \frac{|g|^{2n}}{2n!} \int_{-\infty}^{\infty}dx |h(x)|^{2n} e^{-f(x)} < \infty  \tag{2}\label{eq2}$$
is finite, where in the middle I used the reality of $h$ and $f$ to pull the absolute value sign into the integral.

Note that if we defined odd terms with absolute values inside the integral, the lemma $\eqref{eq1}$ yields the first inequality below, the second following from non-negativity of the integrals:
$$\sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \int_{-\infty}^{\infty}dx |h(x)|^{2n+1} e^{-f(x)} \leq \sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \max\left(\int_{-\infty}^{\infty}dx |h(x)|^{2n} e^{-f(x)}, \int_{-\infty}^{\infty}dx |h(x)|^{2n+2} e^{-f(x)}\right) \leq \sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \left(\int_{-\infty}^{\infty}dx |h(x)|^{2n} e^{-f(x)} + \int_{-\infty}^{\infty}dx |h(x)|^{2n+2} e^{-f(x)}\right)$$
The ratio test applied to the last sum above gives a finite answer by virtue of the limiting ratio being the same as the limiting ratio for the even only series in $\eqref{eq2}$. Thus,
$$\sum_{n=0}^\infty \frac{|g|^{2n+1}}{(2n+1)!} \int_{-\infty}^{\infty}dx |h(x)|^{2n+1} e^{-f(x)} < \infty \tag{3}\label{eq3}$$

Putting this all together by adding $\eqref{eq2}$ and $\eqref{eq3}$, for $|g|<R$, we have
$$\sum_{n=0}^\infty \frac{|g|^{n}}{n!} \int_{-\infty}^{\infty}dx |h(x)|^{n} e^{-f(x)} < \infty $$
However, this is a sufficient condition for Fubini's theorem to allow the swapping of the sum and integral sign for the original problem, so we have at last
$$S(g) \equiv \sum_{n=0}^\infty \frac{g^{n}}{n!} \int_{-\infty}^{\infty}dx h(x)^{n} e^{-f(x)} =   \int_{-\infty}^{\infty}dx \sum_{n=0}^\infty \frac{g^{n}}{n!}  h(x)^{n} e^{-f(x)} = \int_{-\infty}^{\infty}dx e^{-f(x)+g h(x)} \equiv I(g)$$
Thus, when $S(g)$ converges, it is equal to $I(g)$. More carefully stated, when $S(g)$ has a radius of convergence $R$, then $S(g)=I(g)$ when $|g| < R$.

The above assumed $f(x)$ and $g(x)$ real in e.g. $\eqref{eq2}$. The complex case allows for co-opting a famous example from probability as a clear example of $S(g) \neq I(g)$.
Consider the complex functions
$$f(x) = -x^{1/4}+\ln(\sin(x^{1/4})) \text{ for } x\geq 0$$
and
$$g(x) = ix$$
I will take $f(x)$ to be $-\infty$ for $x\leq 0$ for simplicity in the following.
Then the case I am considering is
$$I(g) = \int_{0}^{\infty} e^{-x^{1/4}}\sin(x^{1/4})e^{igx}$$
which is the Fourier transform of the continuous $L^{1}$ function $e^{-x^{1/4}}\sin(x^{1/4}) \theta(x)$, where $\theta(x)$ is Heaviside's step function. This is a unique, invertible Fourier transform. In particular, this Fourier transform cannot merely be the identically $0$ function in $g$, as that would be the Fourier transform of $0$.
However, infamously,
$$\int_{0}^{\infty} e^{-x^{1/4}}\sin(x^{1/4})x^n =0 $$
for all non-negative integers $n$, and so
$$S(g) = 0$$
$S(g)$ is identically zero and hence converges everywhere, but it nevertheless does not agree with $I(g)$ which cannot be zero for all values of $g$ by the invertibility of Fourier transforms. Thus, this gives an example where $S(g)$ converges but to the wrong answer.
