I have a table of points at which a moment generating function is evaluated (for points $t_0,t_1,t_2,\ldots,t_n$ I know $M(t_0), M(t_1), M(t_2),\ldots,M(t_n)$).
I've approximated these tabular function with a rational function.
And I want to recreate PDF from this MGF (because it is unknown).
But clearly not every function can be a moment generating function.
So I wonder what conditions should I add in order to find a rational function that can be a moment generating function?
One possible approach here is as follows.
Note that the m.g.f. $M_{a,b}$ of the gamma distribution with parameters $a,b>0$ is given by the formula $$M_{a,b}(t)=(1-bt)^{-a}$$ for real $t<1/b$ (with $M_{a,b}(t)=\infty$ for real $t\ge1/b$). So, the m.g.f. $M_{a,b}$ is rational on the interval $(-\infty,1/b)$ if $a$ is a positive integer.
For real $c,d>0$, letting now $$N_{c,d}(t)=(1+dt)^{-c}$$ for real $t>-1/b$ (with $N_{c,d}(t)=\infty$ for real $t\le-1/d$), we see that $N_{c,d}$ is the m.g.f. of (the distribution of) the random variable (r.v.) $-X$ such that $X$ is a gamma r.v. with parameters $c,d$.
So, any finite mixture $$\tilde M:=\sum_{i=1}^k p_i M_{a_i,b_i} +\sum_{j=1}^l q_j N_{c_j,d_j} \tag{1}\label{1}$$ of such m.g.f.'s -- with any positive real $p_i$'s and $q_j$'s such that $\sum_{i=1}^k p_i+\sum_{j=1}^l q_j =1$, any positive real numbers $b_i$ and $d_j$, and any positive integers $a_i$ and $c_j$ -- will be a rational m.g.f. on the interval $(-1/\max_{j=1}^l d_j,1/\max_{i=1}^k b_i)$.
So, you will have $k-1+k+l-1+l=2k+2l-2$ positive real "degrees of freedom" and $k+l$ positive integral "degrees of freedom" to fit a rational m.g.f. of the form \eqref{1} to your $n+1$ "data points" $(t_0,M(t_0)),\dots,(t_n,M(t_n))$.
Now, the p.d.f. corresponding to the m.g.f. $\tilde M$ in \eqref{1} is $$\tilde f:=\sum_{i=1}^k p_i f_{a_i,b_i} +\sum_{j=1}^l q_j g_{c_j,d_j},$$ where $f_{a,b}$ and $g_{c,d}$ are the p.d.f.'s corresponding to the m.g.f.'s $M_{a,b}$ and $N_{c,d}$, respectively: $$f_{a,b}(x)=\frac{x^{a-1}e^{-x/b}1(x>0)}{\Gamma(a)b^a},\quad g_{c,d}(x)=\frac{(-x)^{c-1}e^{x/d}1(x<0)}{\Gamma(c)d^c}$$ for real $x$.
To illustrate how flexible model \eqref{1} is, suppose that your data points were $(t_j,M(t_j))$ for $j=0,\dots,8$, where $t_j:=t_0+jh$, $t_0:=-19/12$, $h:=1/3$, and $M$ is the m.g.f. of the mixture of two normal distributions:
one, taken with weight $2/5$, is with mean $-2$ and standard deviation $1/2$;
the other one, taken with weight $1-2/5=3/5$, is with mean $3/2$ and standard deviation $1$.
Let $\tilde M$ be of the form \eqref{1}, with $k=1$, $l=1$, $p_1=0.56$, $q_1=1-p_1$, $a_1=5$, $b_1=0.338$, $c_1=15$, and $d_1=0.128$.
Shown below are the $9$ data points and the graph of $\tilde M$ over the interval $[t_0,t_0+8h]=[-19/12,13/12]$.
We see that even with $k=1$ and $l=1$ in \eqref{1} we get an almost perfect fit.
