Option Explicit

'***************************************************************************
'*                              ver_2.0.1                      2022.10.10  *
'*                                              presented by sainokawa.com *
'***************************************************************************

'===== 非心t分布関数を計算するユーザ定義関数 =====

Function t_Dist_Nc(x As Double, df As Long, nc As Double, dist As Boolean) As Variant

    Dim CDF1 As Double, CDF2 As Double, dx As Double

    Select Case dist
        Case True
            Call t_Dist_Nc_CDF(x, df, nc, CDF1)
            t_Dist_Nc = CDF1
        Case False
            dx = 0.001
            Call t_Dist_Nc_CDF(x - dx, df, nc, CDF1)
            Call t_Dist_Nc_CDF(x + dx, df, nc, CDF2)
            t_Dist_Nc = (CDF2 - CDF1) / 2 / dx
    End Select
    
End Function


'===== 非心t分布の逆関数を計算するユーザ定義関数 =====

Function t_Inv_Nc(prob As Double, df As Long, nc As Double) As Variant

    Dim CDF As Double, dx As Double, i As Integer, j As Long
    
    If prob <= 0 Or prob >= 1 Then t_Inv_Nc = "#NUM!": Exit Function
    Call t_Dist_Nc_CDF(0, df, nc, CDF)
    If prob < CDF Then
        i = -1
    Else: i = 1
    End If
    For j = i To 9999 * i Step i
        t_Inv_Nc = j
        Call t_Dist_Nc_CDF(t_Inv_Nc, df, nc, CDF)
        If i * prob < i * CDF Then
            t_Inv_Nc = t_Inv_Nc - (i + 1) / 2
            Exit For
        End If
    Next
    If j = 10000 Or j = -10000 Then t_Inv_Nc = "#NUM!": Exit Function
    
    For i = 1 To 8
        dx = 10 ^ -i
        For j = 0 To 9
            t_Inv_Nc = t_Inv_Nc + dx
            Call t_Dist_Nc_CDF(t_Inv_Nc, df, nc, CDF)
            If prob < CDF Then
                t_Inv_Nc = t_Inv_Nc - dx
                Exit For
            End If
        Next
    Next

End Function


'----- 非心t分布の累積分布関数を計算するサブルーチン -----
'               (統計数値表 JSA-1972 p.167)

Sub t_Dist_Nc_CDF(x, df, nc, CDF)

    Dim myA As Double, myB As Double, myM() As Double, myMeven As Double, myModd As Double
    Dim ak As Double, k As Long
    ReDim myM(df)
    
    myA = x / Sqr(df): myB = df / (df + x * x)
    myM(0) = myA * Sqr(myB) * Application.WorksheetFunction.Norm_S_Dist(nc * Sqr(myB), False) * _
        Application.WorksheetFunction.Norm_S_Dist(nc * myA * Sqr(myB), True)
    myM(1) = myB * (nc * myA * myM(0) + myA / Sqr(2 * 3.14159265358979) * Application.WorksheetFunction.Norm_S_Dist(nc, False))
    
    myMeven = myM(0): myModd = myM(1)
    For k = 2 To df - 2
        If k = 2 Then
            ak = 1
            Else: ak = 1 / (k - 2) / ak
        End If
        myM(k) = (k - 1) / k * myB * (ak * nc * myA * myM(k - 1) + myM(k - 2))
        If k Mod 2 = 0 Then
            myMeven = myMeven + myM(k)
        Else: myModd = myModd + myM(k)
        End If
    Next
    
    Select Case df Mod 2
        Case 0
            CDF = Application.WorksheetFunction.Norm_S_Dist(-nc, True) + Sqr(2 * 3.14159265358979) * myMeven
        Case 1
            If df = 1 Then myModd = 0
            CDF = Application.WorksheetFunction.Norm_S_Dist(-nc * Sqr(myB), True) + 2 * myFuncT(nc * Sqr(myB), myA) + 2 * myModd
    End Select
    
End Sub


'----- T関数 -----

Function myFuncT(h, a)

    Dim dx As Double, i As Integer, div As Integer
    div = 10000
    
    dx = a / div
    For i = 1 To div
        myFuncT = myFuncT + 1 / (1 + ((i - 0.5) * dx) ^ 2) * Exp(-h * h / 2 * (1 + ((i - 0.5) * dx) ^ 2))
    Next
    myFuncT = myFuncT / 2 / 3.14159265358979 * dx

End Function