numerical integration - More efficient method to compute moments of the Johnson $S_B$ distribution

Here is a very specific feature request. I need

Mean[JohnsonDistribution["SB", γ, δ, 0, 1]]

When I issue e.g.

Table[N[Mean[JohnsonDistribution["SB", γ, δ, 0, 1]]],
      {γ, -10, 10}, {δ, 1, 10}] // TableForm

I get several messages of the form

NIntegrate::ncvb: "NIntegrate failed to converge to prescribed accuracy after
9 recursive bisections in Statistics`DistributionPropertiesDump`x$10142 near
{Statistics`DistributionPropertiesDump`x$10142} = {1.}. NIntegrate obtained
0.9960487815187503` and 0.03645902871201284` for the integral and error estimates."

which seemingly indicates that Mathematica tries to find the mean using numerical integration.

I could live with this somehow, but actually, I also need several higher moments of that distribution (because I want to estimate a distribution with them as the Maximum Likelihood method seems to be too sensitive to fluctuations in my empirical data). And with higher moments, this seems to be even worse - when I am trying to compute an estimated distribution with the method of moments, Mathematica goes on computing for a long time and I stopped it, since I need very many instances of it, and it would be impractical.

So much for the motivation. Now here's what I've been able to find out about this:

In the Appendix to the initial paper where Johnson introduced his distributions, he gives a series expansion for the mean of his $S_B$ distribution based on papers of Mordell in 1920 and 1933 (I could not find any online file for the older one). Johnson also gives expressions of higher moments in terms of derivatives of the mean with respect to one of the parameters, and provides sort of an algorithm to compute them. His formula for the mean reads

$$\sqrt{2\pi}\delta e^{-\frac{\gamma^2}2}\int_{-\infty}^\infty\frac{e^{\pi i\psi t^2-2\pi\nu t}}{1+e^{2\pi t}}dt$$ where $\psi=2\pi i\delta^2$ and $\nu=-\gamma\delta$.

I cannot judge how efficient his procedures might be, but what I know is that in recent years, there has been an explosion of interest in Mordell integrals from the side of number theorists, and they produced lots of new information about it.

In particular, for example, in Theorem 2.1 here one may find the following:

$$\frac12\int_{-\infty}^\infty\frac{e^{\pi i\tau w^2-2\pi zw}}{\cosh(\pi w)}dw=\eta(z,\tau)+e^{\frac{\pi i z^2}\tau}\sqrt{\frac i\tau}\eta(\frac z\tau,-\frac1\tau)$$ where $$\eta(z,\tau)=\sum_{n=0}^\infty\left(\frac{-4}n\right)e^{-\pi inz}e^{-\frac{\pi in^2\tau}4},$$ with $\left(\frac{-4}n\right)$ the Jacobi symbol.

The latter series has very good convergence properties; at least, I am sure it will give a much more efficient method to handle the moments of the Johnson $S_B$ distribution than the numerical integration which Mathematica seemingly uses now.

Will it be possible to incorporate this?

Answer

Using a real value with desired precision for the last parameter:

mom[k_Integer] := Outer[Moment[JohnsonDistribution["SB", ##, 0, 
  SetPrecision[1., 20]], k] &, Range[-10, 10], Range[1, 10]]
mom[1] // MatrixForm