Function ncdft(t値 As Double, df As Double, 非心度 As Double)
   ncdft = cumtnc(t値, df, 非心度)
End Function


Function cumtnc(t As Double, df As Double, pnonc As Double)
'**********************************************************************
'
'     SUBROUTINE CUMTNC(T,DF,PNONC,CUM,CCUM)
'
'                 CUMulative Non-Central T-distribution
'
'DCDFLIB： http://odin.mdacc.tmc.edu/anonftp/page_2.html#DCDFLIB
'Section of Computer Science, Department of Biomathematics, University of Texas M.D. Anderson Hospital.
'の cumtnc.f を excel vba に移植したものである。excel の統計関数を使っている。
'移植に当たって，構造化し goto 文をなくした。 excel 97 を使って確認している。2000/6/11
'
'                              Function
'
'
'     Computes the integral from -infinity to T of the non-central
'     t-density.
'
'
'                              Arguments
'
'
'     T --> Upper limit of integration of the non-central t-density.
'           ｔ値                                   T is DOUBLE PRECISION
'
'     DF --> Degrees of freedom of the non-central t-distribution.
'           自由度                                 DF is DOUBLE PRECISIO
'
'     PNONC --> Non-centrality parameter of the non-central t distibutio
'           非心パラメータ                          PNONC is DOUBLE PRECI
'
'     CUM <-- Cumulative t-distribution.
'           返値：累積非心ｔ分布（確率）
'
'
'                              Method
'
'     Upper tail    of  the  cumulative  noncentral t   using
'     formulae from page 532  of Johnson, Kotz,  Balakrishnan, Coninuous
'     Univariate Distributions, Vol 2, 2nd Edition.  Wiley (1995)
'
'     This implementation starts the calculation at i = lambda,
'     which is near the largest Di.  It then sums forward and backward.
'***********************************************************************
'     .. Parameters ..

      Dim one As Double, zero As Double, half As Double, two As Double
      Dim onep5 As Double
      one = 1#: zero = 0#: half = 0.5: two = 2#: onep5 = 1.5
      Const conv As Double = 0.0000001
      Const tiny As Double = 0.0000000001
'     ..
'     .. Local Scalars ..
      Dim alghdf As Double, b As Double, bb As Double, bbcent As Double
      Dim bcent As Double, cent As Double, d As Double, dcent As Double
      Dim dpnonc As Double, dum1 As Double, dum2 As Double, e As Double
      Dim ecent As Double, halfdf As Double, lambda As Double, lnomx As Double
      Dim lnx As Double, omx As Double, pnonc2 As Double, s As Double
      Dim scent As Double, ss As Double, sscent As Double, t2 As Double
      Dim term As Double, tt As Double, twoi As Double, x As Double
      Dim xi As Double, xlnd As Double, xlne As Double

      Dim cum As Double, ccum As Double
      Dim qrevs As Boolean

'     Case pnonc essentially zero

      If (Abs(pnonc) <= tiny) Then
          cumtnc = 1 - Application.TDist(t, df, 1)
          Exit Function
      End If

      qrevs = t < zero
      If (qrevs) Then
          tt = -t
          dpnonc = -pnonc

      Else
          tt = t
          dpnonc = pnonc
      End If

      pnonc2 = dpnonc * dpnonc
      t2 = tt * tt

      If (Abs(tt) <= tiny) Then
        '  CALL cumnor(-pnonc,cum,ccum)
          cumtnc = Application.NormSDist(-pnonc)
          Exit Function

      End If

      lambda = half * pnonc2
      x = df / (df + t2)
      omx = one - x

      lnx = Log(x)
      lnomx = Log(omx)

      halfdf = half * df
      alghdf = Application.GammaLn(halfdf)

'     ******************** Case i = lambda

      cent = CLng(lambda)

      If (cent < one) Then cent = one

'     Compute d=T(2i) in log space and offset by exp(-lambda)

      xlnd = cent * Log(lambda) - Application.GammaLn(cent + one) - lambda

      dcent = Exp(xlnd)

'     Compute e=t(2i+1) in log space offset by exp(-lambda)

      xlne = (cent + half) * Log(lambda) - Application.GammaLn(cent + onep5) - lambda
      ecent = Exp(xlne)

      If (dpnonc < zero) Then ecent = -ecent

'     Compute bcent=B(2*cent)

     ' CALL bratio(halfdf,cent+half,x,omx,bcent,dum1,ierr)
      bcent = Application.BetaDist(x, halfdf, cent + half)
'     compute bbcent=B(2*cent+1)

  '    CALL bratio(halfdf,cent+one,x,omx,bbcent,dum2,ierr)
       bbcent = Application.BetaDist(x, halfdf, cent + one)
'     Case bcent and bbcent are essentially zero
'     Thus t is effectively infinite

      If ((bcent + bbcent) < tiny) Then
          If (qrevs) Then
              cumtnc = zero
          Else
              cumtnc = one
          End If

          Exit Function

      End If

'     Case bcent and bbcent are essentially one
'     Thus t is effectively zero

  '    If ((dum1 + dum2) < tiny) Then
        If ((1 - (bcent + bbcent)) < tiny) Then
         ' CALL cumnor(-pnonc,cum,ccum)
           cumtnc = Application.NormSDist(-pnonc)
          Exit Function

      End If

'     First term in ccum is D*B + E*BB

      ccum = dcent * bcent + ecent * bbcent

'     compute s(cent) = B(2*(cent+1)) - B(2*cent))

      scent = Application.GammaLn(halfdf + cent + half) - Application.GammaLn(cent + onep5) - alghdf + halfdf * lnx + (cent + half) * lnomx
      scent = Exp(scent)

'     compute ss(cent) = B(2*cent+3) - B(2*cent+1)

      sscent = Application.GammaLn(halfdf + cent + one) - Application.GammaLn(cent + two) - alghdf + halfdf * lnx + (cent + one) * lnomx
      sscent = Exp(sscent)

'     ******************** Sum Forward

      xi = cent + one
      twoi = two * xi

      d = dcent

      e = ecent

      b = bcent

      bb = bbcent

      s = scent

      ss = sscent
      Do
        b = b + s
        bb = bb + ss

        d = (lambda / xi) * d
        e = (lambda / (xi + half)) * e

        term = d * b + e * bb

        ccum = ccum + term

        s = s * omx * (df + twoi - one) / (twoi + one)

        ss = ss * omx * (df + twoi) / (twoi + two)

        xi = xi + one
        twoi = two * xi

      Loop While (Abs(term) > conv * ccum)

'     ******************** Sum Backward

      xi = cent
      twoi = two * xi

      d = dcent

      e = ecent

      b = bcent

      bb = bbcent

      s = scent * (one + twoi) / ((df + twoi - one) * omx)

      ss = sscent * (two + twoi) / ((df + twoi) * omx)

      Do
        b = b - s
        bb = bb - ss

        d = d * (xi / lambda)

        e = e * ((xi + half) / lambda)

        term = d * b + e * bb

        ccum = ccum + term

        xi = xi - one

      If (xi < half) Then Exit Do  'go to 30

      twoi = two * xi

      s = s * (one + twoi) / ((df + twoi - one) * omx)

      ss = ss * (two + twoi) / ((df + twoi) * omx)

      Loop While (Abs(term) > conv * ccum)


      If (qrevs) Then
          cum = half * ccum
          ccum = one - cum

      Else
          ccum = half * ccum
          cum = one - ccum
      End If

'     Due to roundoff error the answer may not lie between zero and one
'     Force it to do so

      cumtnc = Application.Max(Application.Min(cum, one), zero)


End Function