마지막으로 현재까지 내가 인터넷에서 발견한 가장 정밀한 정규 누적분포 구하는 SQL을 보도록 하자. 내용은 아래에서 퍼왔다.

http://home.online.no/~pjacklam/notes/invnorm/

영어로 되어 있으니까 나도 해석은 못하고 프로그래밍된 것만 파악한 거.

정밀도에 대해서는 아래와 같이 써져 있는데 숫자만 봤다.

The absolute value of the relative error is less than 1.15 × 10^-9 in the entire region

전체적으로는 소수 이하 10째자리까지는 아니지만, 어쨌거나 정밀도가 매우 높은 편이다.

Coefficients in rational approximations.
   a(1) <- -3.969683028665376e+01
   a(2) <-  2.209460984245205e+02
   a(3) <- -2.759285104469687e+02
   a(4) <-  1.383577518672690e+02
   a(5) <- -3.066479806614716e+01
   a(6) <-  2.506628277459239e+00

   b(1) <- -5.447609879822406e+01
   b(2) <-  1.615858368580409e+02
   b(3) <- -1.556989798598866e+02
   b(4) <-  6.680131188771972e+01
   b(5) <- -1.328068155288572e+01

   c(1) <- -7.784894002430293e-03
   c(2) <- -3.223964580411365e-01
   c(3) <- -2.400758277161838e+00
   c(4) <- -2.549732539343734e+00
   c(5) <-  4.374664141464968e+00
   c(6) <-  2.938163982698783e+00

   d(1) <-  7.784695709041462e-03
   d(2) <-  3.224671290700398e-01
   d(3) <-  2.445134137142996e+00
   d(4) <-  3.754408661907416e+00

   Define break-points.

   p_low  <- 0.02425
   p_high <- 1 - p_low

   Rational approximation for lower region.

   if 0 < p < p_low
      q <- sqrt(-2*log(p))
      x <- (((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
            ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1)
   endif

   Rational approximation for central region.

   if p_low <= p <= p_high
      q <- p - 0.5
      r <- q*q
      x <- (((((a(1)*r+a(2))*r+a(3))*r+a(4))*r+a(5))*r+a(6))*q /
           (((((b(1)*r+b(2))*r+b(3))*r+b(4))*r+b(5))*r+1)
   endif

   Rational approximation for upper region.

   if p_high < p < 1
      q <- sqrt(-2*log(1-p))
      x <- -(((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
             ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1)
   endif

   The relative error of the approximation has
   absolute value less than 1.15 × 10−9.  One iteration of
   Halley’s rational method (third order) gives full machine precision.

   if 0 < p < 1
      e <- Phi(x) - p                       [use this line
      e <- 0.5 * erfc(-x/sqrt(2)) - p        OR this; not both]
      u <- e * sqrt(2*pi) * exp(x^2/2)
      x <- x - u/(1 + x*u/2)
   endif

무시무시한 공식에 영어까지 섞여 있다. 어쨌거나 저걸 이해하는 것보다는 실제 퍼다 날라서 써먹을 수 있는 게 핵심이니. 위의 공식을 컴퓨터의 각 언어들로 구현한 것들이 링크로 잔뜩 있는데 내가 필요한 것은 오라클용 SQL...

http://home.online.no/~pjacklam/notes/invnorm/impl/smit/stdnormal_cdf.prc

링크 클릭해보면 알겠지만, Stored Procedure 로 되어 있어서 함수 형태로 쓰기는 좋은데... 금융권에서는 SP 형태로는 못쓰는 곳이 많다. 로직이 SP 에 숨어져 있으면 유지보수가 힘들고 하니 정책을 그렇게 정한 것이니, 저걸 바로 쓸 수는 없는 것.

결국 저 SP를 보고 SQL 로 컨버젼 했고, 컨버젼한 결과는 아래에 있다.

SELECT  X AS NRML_DSTRBTN_INP_VAL
      , ROUND(
              CASE WHEN ABS_X <= THRS_SQRT2 THEN Y2_1
                   ELSE
                       CASE WHEN X < 0 THEN Y2_2 ELSE 1 - Y2_2 END
              END
              , 10 -- 소수 10째자리에서 반올림
             )
        AS NRML_CUMUL_DSTRBTN
      , ROUND(STNDRD_NRML_DSTRBTN, 10) AS STNDRD_NRML_DSTRBTN
FROM
    (
        SELECT    0.5 + X
                        * (((( A1 * Z + A2 ) * Z + A3 ) * Z + A4 ) * Z + A5 )
                        / (((( B1 * Z + B2 ) * Z + B3 ) * Z + B4 ) * Z + B5 )
                  AS Y2_1
                , CASE WHEN ABS_X <= SQRT32 THEN
                              Z
                            * (
                                  (((((((( C1 * Y + C2 ) * Y + C3 ) * Y + C4 ) * Y + C5 ) * Y + C6 ) * Y + C7 ) * Y + C8 ) * Y + C9 )
                                / (((((((( D1 * Y + D2 ) * Y + D3 ) * Y + D4 ) * Y + D5 ) * Y + D6 ) * Y + D7 ) * Y + D8 ) * Y + D9 )
                              )
                       ELSE
                              Z
                            * (
                                  (((( P1 * Y + P2 ) * Y + P3 ) * Y + P4 ) * Y + P5 )
                                / (((( Q1 * Y + Q2 ) * Y + Q3 ) * Y + Q4 ) * Y + Q5 )
                              )
                  END
                  AS Y2_2
                , X AS X
                , ABS_X AS ABS_X
                , STNDRD_NRML_DSTRBTN
                , THRS_SQRT2 AS THRS_SQRT2
          FROM
              (
                  SELECT    1.161110663653770E-002 AS A1
                          , 3.951404679838207E-001 AS A2
                          , 2.846603853776254E+001 AS A3
                          , 1.887426188426510E+002 AS A4
                          , 3.209377589138469E+003 AS A5
                          , 1.767766952966369E-001 AS B1
                          , 8.344316438579620E+000 AS B2
                          , 1.725514762600375E+002 AS B3
                          , 1.813893686502485E+003 AS B4
                          , 8.044716608901563E+003 AS B5
                          , 2.15311535474403846E-8 AS C1
                          , 5.64188496988670089E-1 AS C2
                          , 8.88314979438837594E00 AS C3
                          , 6.61191906371416295E01 AS C4
                          , 2.98635138197400131E02 AS C5
                          , 8.81952221241769090E02 AS C6
                          , 1.71204761263407058E03 AS C7
                          , 2.05107837782607147E03 AS C8
                          , 1.23033935479799725E03 AS C9
                          , 1.00000000000000000E00 AS D1
                          , 1.57449261107098347E01 AS D2
                          , 1.17693950891312499E02 AS D3
                          , 5.37181101862009858E02 AS D4
                          , 1.62138957456669019E03 AS D5
                          , 3.29079923573345963E03 AS D6
                          , 4.36261909014324716E03 AS D7
                          , 3.43936767414372164E03 AS D8
                          , 1.23033935480374942E03 AS D9
                          , 1.63153871373020978E-2 AS P1
                          , 3.05326634961232344E-1 AS P2
                          , 3.60344899949804439E-1 AS P3
                          , 1.25781726111229246E-1 AS P4
                          , 1.60837851487422766E-2 AS P5
                          , 6.58749161529837803E-4 AS P6
                          , 1.00000000000000000E00 AS Q1
                          , 2.56852019228982242E00 AS Q2
                          , 1.87295284992346047E00 AS Q3
                          , 5.27905102951428412E-1 AS Q4
                          , 6.05183413124413191E-2 AS Q5
                          , 2.33520497626869185E-3 AS Q6
                          , ABS_X AS ABS_X
                          , CASE WHEN ABS_X <= THRS_SQRT2 THEN X2
                                 ELSE
                                     CASE WHEN ABS_X > SQRT32
                                              THEN EXP( -X2 / 2 ) / 2 * SQRT2 / ABS_X
                                          ELSE EXP( -X2 / 2 ) / 2
                                     END
                            END Z
                          , CASE WHEN ABS_X <= SQRT32
                                     THEN ABS_X / SQRT2
                                 ELSE 2 / X2
                            END AS Y
                          , X AS X
                          , THRS_SQRT2 AS THRS_SQRT2
                          , SQRTPI AS SQRTPI
                          , SQRT2 AS SQRT2
                          , SQRT32 AS SQRT32
                          , 1 / (SQRTPI * SQRT2) * EXP(-X2 / 2) AS STNDRD_NRML_DSTRBTN
                  FROM      DUAL
                          , (
                                SELECT    SQRT(2)           AS SQRT2
                                        , SQRT(32)          AS SQRT32
                                        , SQRT(ACOS(-1))    AS SQRTPI
                                        , ABS(:InpVal)      AS ABS_X
                                        , :InpVal           AS X
                                        , POWER(:InpVal, 2) AS X2
                                        , SQRT(2) * 15/32   AS THRS_SQRT2
                                FROM    DUAL
                            ) B
              )
    )
;

결과를 보자.

엑셀에서의 결과로 검증해보자.

맞는 거 같다.

하나 더.

자 이제 저 소스코드만 퍼가고, :InpVal 를 입력해주면 표준정규분포 확률밀도, 그리고 누적정규분포까지 계산해서 보여준다. 그럼 끝.