# 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 StatisticsDistributionPropertiesDumpx$10142 near {StatisticsDistributionPropertiesDumpx$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.

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?

• FWIW: Mean[] and Expectation[]/NExpectation[] don't work, but the explicit integral works nicely with an appropriate choice of method: Table[(δ/Sqrt[2 π]) NIntegrate[Exp[-(γ + δ Log[x/(1 - x)])^2/2]/(1 - x), {x, 0, 1}, Method -> "DoubleExponential"], {γ, -10, 10}, {δ, 1, 10}]. Oct 24, 2015 at 11:11
• @J.M. Wow this works very efficiently, also with all higher moments (multiplying integrand by $x^k$)! Still, it would be good to incorporate this somehow, I believe an average user would hardly come up with this. Could you please make this an answer? Oct 24, 2015 at 11:17
• As you have noted, the moments of the Johnson $S_B$ system are extremely complicated. As such, the method of moments is not generally used for fitting $S_B$ systems. Instead, a method of percentile points is usually used, which equates percentile points of the observed and fitted curves. The latter is not an exact methodology. I have come across some new alternative fitting methods in recent years. One can also always use Pearson fitting functions as a substitute, which is inevitably a much simpler strategy. Oct 24, 2015 at 14:47
• @wolfies, extremely complicated indeed; I encountered special functions I haven't seen in a while in the papers the OP linked to. Oct 24, 2015 at 15:00
• For your perusal, here is Egon Pearson's original paper. Oct 24, 2015 at 15:08

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 // MatrixForm • Oh my so all I had to do is increase precision! Oct 3, 2018 at 20:43

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.

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}


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]]]]


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]]

• This must be doubtlessly accepted, although it must be called a workaround :D Oct 24, 2015 at 11:52
• Agreed, a mere workaround. Hopefully this gets noticed by the appropriate folks. Oct 24, 2015 at 11:58
• Actually (although this is not part of the question per se) it would be great if you could also indicate how to pass this to the ParameterEstimator->"MethodOfMoments" option of the EstimatedDistribution... Oct 24, 2015 at 12:02
• Ah, that looks like a harder question; nothing obvious in the docs. Consider asking it as a separate question. Oct 24, 2015 at 12:07
• OK will try to formulate. And many thanks once more, you helped me out a lot! Oct 24, 2015 at 12:08

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} *)

• Good to know! Is there a way to pass this to the ParameterEstimator -> "MethodOfMoments" option of EstimatedDistribution or FindDistributionParameters? (I've asked this in another question) Oct 24, 2015 at 16:16
• Dear ilian, would it not have been more expedient to have NExpectation[] inherit NIntegrate[]'s options? Nevertheless, thanks for pointing this syntax out. Oct 24, 2015 at 16:22