It is always possible to do MLE. The general approach is to assume a distribution type (aka family) for the underlying random process. Say the distribution PDF is ##f(x;\vec a)## where ##\vec a## is a vector of unknown parameters. The functional form of ##F## has to be assumed and the parameters ##a_j## are to be estimated.
We have a set of observations which in this case are the pairs ##(M_k,N_k)##. To do MLE, we just need to find a formula for the likelihood of each observation, which will be a probability or a probability density for ##N_k## given ##M_k##. Then we construct the overall likelihood which is the product of the likelihoods of the individual daily observations:
$$L(\vec N;\vec M, \vec a) = \prod_{k=1}^n L(N_k; M_k,\vec a)$$
We take partial derivatives of this wrt each ##a_j##, set them to zero and so get ##m## equations in ##m## unknowns, where ##m## is the dimensionality of ##\vec a## (number of parameters). We solve the system to find ##\vec a##.
In case of this casino problem the underlying distribution is that of the payout on a single play. The only thing we know about ##f## is that its integral on ##[0,\infty)## must be 1. We make a pragmatic choice based on what shape we might expect and how many observations we have. The fewer observations we have, the smaller should be the number of parameters, to avoid overfitting.
Say we decide the observation set is large enough to justify ##m## parameters. Then a suitable functional form would be
$$f(x;\vec a)=a_m \delta(x) + C\exp\left(-x^m + \sum_{j=1}^{m-1} a_j x^j\right)$$
The ##\delta(x)## is the Dirac delta function. We use it to allow a nonzero probability ##a_m## of no payout.
The parameter ##C## is chosen to make the integral 1, so it is calculated from the other parameters, rather than estimated. This distribution is always positive and will tend asymptotically to zero as ##x\to\infty##, which is what we would expect for a payoff.
Given that, the likelihood of the observation ##(M_k,N_k)## is
$$L(N_k;M_k,\vec a) = \int_0^{D_k-B_1} \int_0^{D_k - B_2}
...
\int_0^{D_k - B_{N_k}} \prod_{r=1}^{N_k} f(x_r;\vec a)
dx_{N_k}...dx_1
$$
where ##D_k\triangleq N_k-M_k+1## and ##B_j\triangleq \sum_{r=1}^{j-1} x_r##. Note that ##B_1=0## and ##B_2=x_1##.
In solving the MLE, we apply the constraint that ##a_m\in [0,1]## as that is the probability of a zero payout. All the other ##a_j## can be any real number.
Once we have solved the MLE, we have a full estimated specification of the PDF, so we can integrate that to find the expected payoff from a single play.
Given the problem has been specified in recent posts to have continuous rather than discrete support, talk of a 'paytable' is misplaced. Tables only make sense for discrete sample spaces. We are estimating a function, not a table, and we estimate it by assuming a suitable form, then estimating its parameters.
The minimum number of parameters that can be used in this method is one, in which case the parameter is the probability of no payout on a play and the distribution of nonzero payouts is exponential with parameter 1. As the number of observations increases (the number of days for which ##(M_k,N_k)## was observed), the number of parameters can be increased, but I prefer to err on the side of under-fitting, so would only add another parameter after a load more observations.