```
override def cumulativeProbability(x: Double): Double = {
    // An approximation to the true CDF since numerical integration diverges.
    // Source and notation: https://www.sciencedirect.com/science/article/pii/S209012321000069X.
    val standardizedX = (x - location) / scale

    if (standardizedX < -3) 0.0
    else if (-3 <= standardizedX && standardizedX < -1) H1(standardizedX)
    else if (-1 <= standardizedX && standardizedX < 1) H2(standardizedX)
    else if (1 <= standardizedX && standardizedX < 3) H3(standardizedX)
    else 1.0
  }

private def H1(x: Double): Double = {
    val first = (3 * SHAPESQ + 9) * gI(-3 / shape, x)
    val second = shape * (3 + 2 / 3 * SHAPESQ) * math.exp(-9 / (2 * SHAPESQ))
    val third = -shape * (9 + 2 / 3 * SHAPESQ + 3 * shape * x + SHAPESQ * x * x / 3) * math.exp(-x * x / 2)
    (first + second + third) / (8 * SQRT2PI)
  }

  private def H2(x: Double): Double = {
    val first = (3 * SHAPESQ + 9) * gI(-3 / shape, -1 / shape)
    val second = 8 * gI(-1 / shape, x)
    val third = shape * (3 + 2 / 3 * SHAPESQ) * math.exp(-9 / (2 * SHAPESQ))
    val fourth = -shape * (1 + 2 * SHAPESQ) * math.exp(-1 / (2 * SHAPESQ))
    val fifth = 2 / 3 * shape * (-9 + 2 * SHAPESQ + SHAPESQ * x * x) * math.exp(-x * x / 2)
    (first + second + third + fourth + fifth) / (8 * SQRT2PI)
  }

  private def H3(x: Double): Double = {
    val first = (3 * SHAPESQ + 9) * gI(-3 / shape, -1 / shape)
    val second = 8 * gI(-1 / shape, 1 / shape)
    val third = (7 - 3 * SHAPESQ) * gI(1 / shape, x)
    val fourth = shape * (3 + 2 / 3 * SHAPESQ) * math.exp(-9 / (2 * SHAPESQ))
    val fifth = -shape / 3 * (27 - 9 * shape * x + SHAPESQ * x * x + 2 * SHAPESQ) * math.exp(-x * x / 2)
    (first + second + third + fourth + fifth) / (8 * SQRT2PI)
  }

  // Integrate exp(-x^2 / 2) from a to b.
  private def gI(a: Double, b: Double): Double = SQRT2PI / 2 * Erf.erf(a / SQRT2PI, b / SQRT2PI)
```