What's the best way to approximate
Ordinarily, I would use a normal approximation with integer correction...
...if the skewness and excess kurtosis were (closer to) 0, but that's not the case here.
I have to perform multiple approximations for different discrete distributions with different values of
This is an interesting question, which doesn't really have a good solution. There a few different ways of tackling this problem.
Assume an underlying distribution and match moments - as suggested in the answers by @ivant and @onestop. One downside is that the multivariate generalisation may be unclear.
Saddlepoint approximations. In this paper:
Gillespie, C.S. and Renshaw, E. An improved saddlepoint approximation. Mathematical Biosciences, 2007.
We look at recovering a pdf/pmf when given only the first few moments. We found that this approach works when the skewness isn't too large.
Laguerre expansions:
Mustapha, H. and Dimitrakopoulosa, R. Generalized Laguerre expansions of multivariate probability densities with moments. Computers & Mathematics with Applications, 2010.
The results in this paper seem more promising, but I haven't coded them up.
Fitting a distribution to data using the first four moments is exactly what Karl Pearson devised the Pearson family of continuous probability distributions for (maximum likelihood is much more popular these days of course). Should be straightforward to fit the relevant member of that family then use the same type of continuity correction as you give above for the normal distribution.
I assume you must have a truly enormous sample size? Otherwise sample estimates of skewness and especially kurtosis are often hopelessly imprecise, as well as being highly sensitive to outliers.In any case, I highly recommend you have a look at L-moments as an alternative that have several advantages over ordinary moments that can be advantageous for fitting distributions to data.
You could try to use skew normal distribution and see if excess kurtosis for your particular data sets is sufficiently close to the excess kurtosis of the distribution for given skewness. If it is, you can use the skew normal distribution cdf to estimate the probability. If not, you would have to come up with a transformation to the normal/skew pdf similar to the one used for the skew normal distribution, which would give you control over both skewness and excess kurtosis.