More efficient method to compute moments of the Johnson $S_B$ distribution
I'll preface this answer first with a complaint:
NExpectation[]
andNProbability[]
are notsufficiently resilientobviously 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]]
As a minor addition to J.M.'s excellent answer,
If something breaks while evaluating
NExpectation[]
orNProbability[]
, a number of the things that otherwise can be adjusted inNIntegrate[]
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} *)