'Begin Description
' 非心F分布累積確率を求める。使用法 =ncdff(F値, df1, df2, 非心度パラメータ)
'Function ncdff は Function cumfnc を呼び出しているだけです。
'spss と名前を合わせました。= を押したときに日本語でF値, df1, df2, 非心度が
'でるようにしました。df1 は分子，df2 は分母の自由度です。どちらを呼び出して
'も同じですが，忘れたときにはncdff のほうが便利ですのでこちらをよびだしましょう。
'Function qsmall は精度を決めてます。
'End Description


Option Explicit

Function ncdff(F値 As Double, df1 As Double, df2 As Double, 非心度 As Double)
   ncdff = cumfnc(F値, df1, df2, 非心度)
End Function

Function cumfnc(f As Double, dfn As Double, dfd As Double, pnonc As Double)


'**********************************************************************
'
'               F -NON- -C-ENTRAL F 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.
'の cumfnc.f を excel vba に移植したものである。excel の統計関数を使っている。
'移植に当たって，構造化し goto 文をなくした。 excel 97 を使って確認している。2000/5/24
'
'                              Function
'
'
'     COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
'     DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC
'
'
'                              Arguments
'
'
'     F --> 非心Ｆ分布のＦ値　累積を求める上限F値
'
'     DFN --> DEGREES OF FREEDOM OF NUMERATOR　分子自由度
'
'     DFD -->  DEGREES OF FREEDOM OF DENOMINATOR　分母自由度
'
'     PNONC --> NONCENTRALITY PARAMETER.非心度パラメータ
'
'     CUM <-- CUMULATIVE NONCENTRAL F DISTRIBUTION　非心F分布累積確率(これを返す)
'
'
'
'                              Method(もとプログラムに書いてあったもの)
'
'
'     USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
'     SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
'     (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
'     THE CONVERGENCE CRITERION IS MET.
'
'     FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
'     BY FORMULA 26.5.16.
'
'
'               REFERENCE
'
'
'     HANDBOOD OF MATHEMATICAL FUNCTIONS
'     EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
'     NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
'     MARCH 1965
'     P 947, EQUATIONS 26.6.17, 26.6.18
'     http://odin.mdacc.tmc.edu/anonftp/page_2.html#DCDFLIB 2000/5/20
'
'                              Note
'
'
'     THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
'     TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20).  EPS IS
'     SET TO 1.0E-6 IN A DATA STATEMENT WHICH CAN BE CHANGED.
'
'**********************************************************************

      Dim dsum As Double, dummy As Double, prod As Double, xx As Double, yy As Double
      Dim adn As Double, aup As Double, B As Double, betdn As Double, betup As Double, centwt As Double, dnterm As Double, eps As Double, sum As Double
      Dim upterm As Double, xmult As Double, xnonc As Double, x As Double
      Dim i As Long, icent As long, ierr As Integer
      Dim cum As Double

      Const half As Double = 0.5
      Const done As Double = 1#

      If ((f <= 0#)) Then
        cumfnc = 0#
        Exit Function
      End If
    
      If (pnonc < 0.0000000001) Then
'
'     Handle case in which the non-centrality parameter is
'     (essentially) zero.

        cumfnc = 1 - Application.FDist(f, dfn, dfd)
        Exit Function
      End If
    
      xnonc = pnonc / 2#

'     Calculate the central term of the poisson weighting factor.

      icent = xnonc
      If (icent = 0) Then icent = 1

'     Compute central weight term

      centwt = Exp(-xnonc + icent * Log(xnonc) - Application.GammaLn(CDbl(icent + 1)))

'     Compute central incomplete beta term
'     Assure that minimum of arg to beta and 1 - arg is computed
'          accurately.

      prod = dfn * f
      dsum = dfd + prod
      yy = dfd / dsum
      If (yy > half) Then
          xx = prod / dsum
          yy = done - xx

      Else
          xx = done - yy
      End If

'      CALL bratio(dfn*half+cdbl(icent),dfd*half,xx,yy,betdn,dummy,ierr)
      betdn = Application.BetaDist(xx, dfn * half + CDbl(icent), dfd * half)

      adn = dfn / 2# + CDbl(icent)
      aup = adn
      B = dfd / 2#
      betup = betdn
      sum = centwt * betdn

'     Now sum terms backward from icent until convergence or all done

      xmult = centwt
      i = icent
      dnterm = Exp(Application.GammaLn(adn + B) - Application.GammaLn(adn + 1#) - Application.GammaLn(B) + adn * Log(xx) + B * Log(yy))

      Do Until (qsmall((xmult * betdn), sum) Or i <= 0)
        xmult = xmult * (i / xnonc)
        i = i - 1
        adn = adn - 1
        dnterm = (adn + 1) / ((adn + B) * xx) * dnterm
        betdn = betdn + dnterm
        sum = sum + xmult * betdn
      Loop
      
      i = icent + 1

'     Now sum forwards until convergence

      xmult = centwt
      If ((aup - 1 + B) = 0) Then
    upterm = Exp(-Application.GammaLn(aup) - Application.GammaLn(B) + (aup - 1) * Log(xx) + B * Log(yy))

      Else
          upterm = Exp(Application.GammaLn(aup - 1 + B) - Application.GammaLn(aup) - Application.GammaLn(B) + (aup - 1) * Log(xx) + B * Log(yy))
      End If


      Do
        xmult = xmult * (xnonc / i)
        i = i + 1
        aup = aup + 1
        upterm = (aup + B - 2#) * xx / (aup - 1) * upterm
        betup = betup - upterm
        sum = sum + xmult * betup
      Loop Until (qsmall(xmult * betup, sum))

      cumfnc = sum

End Function

Private Function qsmall(x As Double, sum As Double) As Boolean

     Dim eps As Double
       eps = 0.000001 '1e-6 
'     ..
'     .. Statement Function definitions ..
      qsmall = (sum < 1E-20) Or (x < eps * sum)
End Function