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

Question

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 1

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

Answer 2

I'll preface this answer first with a complaint:

NExpectation[] and NProbability[] are not sufficiently resilient obviously adjustable.

Ideally, these two functions are an "interface" to NIntegrate[], allowing the user to formulate his expression purely in distributional terms. Unfortunately, when one hits cases like this, the things one might usually fiddle with in NIntegrate[] seem to be absent in NExpectation[] and NProbability[]!

Options[NExpectation]
   {AccuracyGoal -> ∞, Compiled -> Automatic, Method -> Automatic,
    PrecisionGoal -> Automatic, WorkingPrecision -> MachinePrecision}

Options[NIntegrate]
   {AccuracyGoal -> ∞, Compiled -> Automatic, EvaluationMonitor -> None,
    Exclusions -> None, MaxPoints -> Automatic, MaxRecursion -> Automatic,
    Method -> Automatic, MinRecursion -> 0, PrecisionGoal -> Automatic,
    WorkingPrecision -> MachinePrecision}

See the difference? If something breaks while evaluating NExpectation[] or NProbability[], a number of the things that otherwise can be adjusted in NIntegrate[] aren't there. This forces the user to fall back on NIntegrate[], and ruefully wonder why he even thought of trying the fancy syntactic sugar in the first place.

(But, see the addenda below.)

---

Having said all this, the $S_B$ distribution is apparently one of those distributions that require the use of NIntegrate[] proper. To help us along, we display the corresponding probability density function:

PDF[JohnsonDistribution["SB", γ, δ, 0, 1], x]
   Piecewise[{{δ/(E^((γ + δ Log[x/(1 - x)])^2/2) Sqrt[2 π] (1 - x) x), 0 < x < 1}}, 0]

Thus, to assemble the $k$-th moment, we multiply this expression with $x^k$. The logistic function and the denominator present in the PDF will give the default quadrature method a spot of trouble, so we switch to a method that is relatively robust towards endpoint singularities: the double exponential quadrature of Takahashi and Mori:

SetAttributes[sbMoment, Listable];
sbMoment[k_Integer?NonNegative, γ_?NumericQ, δ_?NumericQ, opts___] := 
         Module[{prec = Precision[{γ, δ}]},
                If[prec === ∞, prec = MachinePrecision];
                (δ/Sqrt[2 π]) NIntegrate[x^(k - 1)
                   Exp[-(γ + δ Log[x/(1 - x)])^2/2]/(1 - x), {x, 0, 1},
                   Method -> "DoubleExponential", opts, WorkingPrecision -> prec]]

Examples:

sbMoment[Range[0, 5], -10, 1]
   {0.999999999999918, 0.999925163391599, 0.999850341996451,
    0.9997755358060877, 0.9997007448121356, 0.9996259690062346}

sbMoment[Range[0, 2], -5, 3, WorkingPrecision -> 20]
   {1.0000000000000000000, 0.83615201160671172434, 0.70122109224906876673}

---

Addendum

As ilian notes in his answer, it's actually possible to pass NIntegrate[]'s options to NExpectation[]; nevertheless, the syntax is not as transparent as I'd like. With that,

SetAttributes[sbMoment, Listable];
sbMoment[k_Integer?NonNegative, γ_?NumericQ, δ_?NumericQ, opts___] := 
         Module[{prec = Precision[{γ, δ}]},
                If[prec === ∞, prec = MachinePrecision];
                NExpectation[\[FormalX]^k, \[FormalX] \[Distributed]
                             JohnsonDistribution["SB", γ, δ, 0, 1], 
                             Method -> {"NIntegrate",
                             FilterRules[{Method -> "DoubleExponential",
                                          opts}, Options[NIntegrate]]},
                             Sequence @@ FilterRules[{opts,
                                                      WorkingPrecision -> prec}, 
                                                     Options[NExpectation]]]]

---

Addendum II

As it turns out, there is a formulation of the Johnson $S_B$ moments that leads to an even more efficient evaluation routine.

Draper, in his paper, gives an equivalent, but more computationally efficient, integral formula for the moments:

$$\small \mathtt{Moment[JohnsonDistribution["SB"},\gamma,\delta,0,1\mathtt{]},k\mathtt{]}=\frac1{\sqrt{2\pi}}\int\limits_{-\infty}^\infty \frac{\exp\left(-\frac{t^2}{2}\right)}{\left(\exp\left(-\frac{t-\gamma}{\delta}\right)+1\right)^k}\;\mathrm dt$$

The good thing about this is that this is the precise sort of integral that the trapezoidal rule is very efficient at evaluating; see this paper and its references for further discussion. If speed is truly critical, one might be able to write a Compile[]-d function for evaluating the trapezoidal sum involved, but NIntegrate[] already gave results that were just as good as the results from the previous two versions, and in much less time:

SetAttributes[sbMoment, Listable]
sbMoment[k_Integer?NonNegative, γ_?NumericQ, δ_?NumericQ, opts___] := 
 Module[{prec = Precision[{γ, δ}]}, 
        If[prec === ∞, prec = MachinePrecision];
        NIntegrate[PDF[NormalDistribution[], t]/(Exp[(γ - t)/δ] + 1)^k, {t, -∞, ∞},
                   Method -> {"Trapezoidal", "SymbolicProcessing" -> 0},
                   opts, WorkingPrecision -> prec]]

Answer 3

As a minor addition to J.M.'s excellent answer,

If something breaks while evaluating NExpectation[] or NProbability[], a number of the things that otherwise can be adjusted in NIntegrate[] aren't there.

Options can be passed to NIntegrate, for example try something like

Table[NExpectation[X, X \[Distributed] JohnsonDistribution["SB", γ, δ, 0, 1], 
       Method -> {"NIntegrate", Method -> "DoubleExponential"}, WorkingPrecision -> 20], 
      {γ, -10, 10}, {δ, 1, 10}]

also for higher moments

NExpectation[X^Range[0, 5], X \[Distributed] JohnsonDistribution["SB", -5, 3, 0, 1], 
     Method -> {"NIntegrate", Method -> "DoubleExponential"}, WorkingPrecision -> 20]

(* {1.0000000000000000000, 0.83615201160671172434, 0.70122109224906876673,
    0.58972935904106407569, 0.49730884869043345102, 0.42046180967574075322} *)