com.quantarray.anaheim.numerics.NormalDistribution2.scala Maven / Gradle / Ivy

Go to download
package com.quantarray.anaheim.numerics

/**
 * Bivariate normal distribution.
 */
object NormalDistribution2
{
  private[this] def N: Double => Double = StandardNormalDistribution.cdf

  private[this] val TwoPi = 2 * math.Pi
  private[this] val RootTwoPi = math.sqrt(TwoPi)

  private[this] val W = Array(
    // Gauss Legendre points and weights, n =  6
    Array(0.1713244923791705, 0.3607615730481384, 0.4679139345726904),

    // Gauss Legendre points and weights, n = 12
    Array(0.4717533638651177e-1, 0.1069393259953183, 0.1600783285433464,
      0.2031674267230659, 0.2334925365383547, 0.2491470458134029),

    // Gauss Legendre points and weights, n = 20
    Array(0.1761400713915212e-1, 0.4060142980038694e-1, 0.6267204833410906e-1,
      0.8327674157670475e-1, 0.1019301198172404, 0.1181945319615184,
      0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
      0.1527533871307259)
  )

  private[this] val X = Array(
    // Gauss Legendre points and weights, n =  6
    Array(0.9324695142031522, 0.6612093864662647, 0.2386191860831970),

    // Gauss Legendre points and weights, n = 12
    Array(0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
      0.5873179542866171, 0.3678314989981802, 0.1252334085114692),

    // Gauss Legendre points and weights, n = 20
    Array(0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
      0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
      0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
      0.7652652113349733e-1)
  )

  def cdf(x: Double, y: Double, corr: Double): Double =
  {
    var lg: Int = 0
    var ng: Int = 0

    if (math.abs(corr) < 0.3)
    {
      ng = 0
      lg = 3
    }
    else if (math.abs(corr) < 0.75)
    {
      ng = 1
      lg = 6

    }
    else
    {
      ng = 2
      lg = 10
    }

    var h = -x
    var k = -y
    var hk = h * k
    var bvn = 0.0

    if (math.abs(corr) < 0.925)
    {
      val hs = (h * h + k * k) / 2.0
      val asr = math.asin(corr)
      for (i <- 0 until lg)
      {
        val sn1 = math.sin(asr * (1.0 - X(ng)(i)) / 2.0)
        bvn += W(ng)(i) * math.exp((sn1 * hk - hs) / (1.0 - sn1 * sn1))
        val sn2 = math.sin(asr * (1.0 + X(ng)(i)) / 2.0)
        bvn += W(ng)(i) * math.exp((sn2 * hk - hs) / (1.0 - sn2 * sn2))
      }

      bvn = bvn * asr / (2 * TwoPi) + N(-h) * N(-k)
    }
    else
    {
      if (corr < 0.0)
      {
        k = -k
        hk = -hk
      }

      var xs = 0.0

      if (math.abs(corr) < 1.0)
      {
        val as = (1.0 - corr) * (1.0 + corr)
        var a = math.sqrt(as)
        val bs = (h - k) * (h - k)
        val c = (4.0 - hk) / 8.0
        val d = (12.0 - hk) / 16.0
        var asr = -(bs / as + hk) / 2.0

        if (asr > -100.0)
          bvn = a * math.exp(asr) * (1.0 - c * (bs - as) * (1.0 - d * bs / 5.0) / 3.0 + c * d * as * as / 5.0)

        if (-hk < 100.0)
        {
          val b = math.sqrt(bs)
          val sp = RootTwoPi * N(-b / a)
          bvn = bvn - math.exp(-hk / 2.0) * sp * b * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0)
        }

        a = a / 2.0

        for (i <- 0 until lg)
        {
          for (is <- -1 until 1 by 2)
          {
            xs = a * (is * X(ng)(i) + 1.0)
            xs = xs * xs
            val rs = math.sqrt(1.0 - xs)
            asr = -(bs / xs + hk) / 2.0
            if (asr > -100.0)
            {
              val sp = 1.0 + c * xs * (1.0 + d * xs)
              val ep = math.exp(-hk * (1.0 - rs) / (2.0 * (1.0 + rs))) / rs
              bvn += a * W(ng)(i) * math.exp(asr) * (ep - sp)
            }
          }
        }

        bvn = -bvn / TwoPi
      }

      if (corr > 0.0)
      {
        if (k > h)
          h = k
        bvn += N(-h)
      }

      if (corr < 0.0)
      {
        xs = N(-h) - N(-k)
        if (xs < 0.0)
          xs = 0.0
        bvn = -bvn + xs
      }
    }

    if (bvn <= 0.0)
      return 0.0

    if (bvn >= 1.0)
      return 1.0

    bvn
  }
}