Unify conventions for divergent and undefined moments

[mmikhasenko]

There seems to be no single convention for moments that diverge or do not exist.

For example, current behavior mixes Inf and NaN across heavy-tailed distributions:

mean(Levy()) == Inf
var(TDist(2)) == Inf

while skewness/kurtosis often return NaN.

Related prior discussion:

• Boundary parameter handling #283 mentions choosing between NaN and throwing for skewness/kurtosis in boundary-parameter cases. This issue is about the same convention question for divergent/nonexistent moments more generally.

Code references:

• Cauchy: mean, var, skewness, and kurtosis all return NaN.
• Levy: mean/var return Inf, while skewness/kurtosis return NaN.
• TDist: var and kurtosis can return Inf, while mean, skewness, and kurtosis can return NaN in other parameter regimes.
• Pareto and InverseGamma: lower moments return Inf when divergent, but higher standardized moments return NaN.
• MvTDist: mean, var, and cov use NaN when moments do not exist.

Could we document and apply a consistent rule package-wide?

For example: use Inf when the corresponding raw/central moment diverges to infinity,
and reserve NaN or an error for genuinely undefined standardized quantities.