The sinh-arcsinh (SHASH) distribution has a pdf as follows:

$f(x) = {\delta cosh(\omega)\over \sqrt{1+({x-\theta \over \sigma})^2}}\phi[sinh(\omega)]$ where $\omega=\gamma+\delta sinh^{-1}({x-\theta \over \sigma})$

and $\phi()$ is the standard normal pdf.

Note: $-\infty<\gamma,x,\theta<\infty; 0<\delta, \sigma$

Does anyone know what the expected value and variance of it is? I can see it is closely related to the Johnson SU distribution which has the mean, variance, median, and quantiles analytically defined. It also has the Normal distribution as a special case when $\gamma = 0$ and $ \delta = 1$. I know some distributions don't have these available in closed form analytic expressions....but I don't know if that's the case here. If available, it would be very helpful if someone knows what they are!

Update: It appears the moments are published, but with the caveat that they depend on the modified Bessel function of the second kind ($P_\nu)$ below. E.g.

$E(X_{\gamma, \delta})=-sinh(\gamma/ \delta)P_{1/\delta}$

$Var(X_{\gamma, \delta})={1 \over 2}(cosh({2\gamma \over \delta})P_{2/\delta} - 1) - \mu_{\gamma, \delta}^2$

However -- this doesn't take into account the center, scale in the full version ($\theta, \sigma$). Can someone help modify the mean and var with center and scale included?


1 Answer 1


Indeed, the moments of this distribution are already calculated in the original paper


For the case with location and scale parameters, you just need to use the usual properties of the location-scale family. Let $Z = \mu + \sigma X$, where $X \sim SHASH(\gamma,\delta)$, then

$$E[Z] = \mu + \sigma E[X].$$


$$Var[Z] = \sigma^2 Var[X].$$

[Here] is an implementation of this distribution in R. Using their command rsas, you can simulate from a distribution with specific parameter values and approximate its mean and variance using mean() and var().

sim <- rsas(1e5,3,2,-1,1.5)

Alternatively, you can use their function dsas to approximate the mean using numerical integration:

mean_sas <- function(mu,sigma,epsilon,delta){
  tempf <- Vectorize(function(x) x*dsas(x,mu,sigma,epsilon,delta))
  val <- integrate(tempf,-Inf,Inf)$value


Similarly for the variance:

var_sas<- function(mu,sigma,epsilon,delta){
  tempf1 <- Vectorize(function(x) x*dsas(x,mu,sigma,epsilon,delta))
  tempf2 <- Vectorize(function(x) x^2*dsas(x,mu,sigma,epsilon,delta))
  val <- integrate(tempf2,-Inf,Inf)$value - (integrate(tempf1,-Inf,Inf)$value)^2


Finally, the same ideas can be applied to [their alternative version] of the SHASH distribution.

  • $\begingroup$ +1 for including code! $\endgroup$
    – JPJ
    Commented Feb 9, 2020 at 21:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.