计算非中心 t 分布的累积分布函数之 WPS 下的实现

之前的实现在 WPS 表格下得到的结果和 Excel 中的不一致，经过测试，是 WPS 下 Possion_Dist 和 Beta_Dist 与 Excel 中的实现不一致的原因。
修正后，检测是否运行在 WPS 环境中（通过 Application.Caption），并根据环境使用不同的函数名（如 Poisson vs Poisson_Dist）。

Function NT_DIST(t As Double, df As Double, d As Double, cum As Boolean, _
        Optional iter As Long = 1000, Optional prec As Double = 0.000000000001) As Double
    
    ' return non-central t pdf if cum = FALSE or distribution if cum = TRUE at t with df degrees of freedom and
    ' non-centrality parameter d; iter is the number of iterations and prec is the precision of the answer
    
    If cum Then   ' distribution function
        If t >= 0 Then
            NT_DIST = NT_DIST_CDF(t, df, d, iter, prec)
        Else
            NT_DIST = 1 - NT_DIST_CDF(-t, df, -d, iter, prec)
        End If
        
    Else        ' pdf
        If Abs(t) < 0.00000001 Or (Abs(t) < 0.0000001 And d <= 0.35) Then ' t is close enough to zero
            NT_DIST = Exp(Application.WorksheetFunction.GammaLn((df + 1) / 2) - Application.WorksheetFunction.GammaLn(df / 2) - d ^ 2 / 2) / _
                Application.WorksheetFunction.SqrtPi(df)
        Else
            NT_DIST = df * (NT_DIST(t * Sqr(1 + 2 / df), df + 2, d, True, iter, prec) - NT_DIST(t, df, d, True, iter, prec)) / t
        End If
    End If
    
End Function


Function NT_DIST_CDF(t As Double, df As Double, d As Double, iter As Long, prec As Double) As Double
    ' return non-central t cdf at t >= 0 with df degrees of freedom and non-centrality parameter d
    ' iter is the number of iterations and prec is the precision of the answer
    Dim lambda As Double    ' mean of Poisson distribution
    lambda = d ^ 2 / 2
    Dim i As Long
    Dim df2 As Double       ' df/2
    df2 = df / 2
    Dim dsq As Double       ' d/sqrt(2)
    dsq = d / Sqr(2)
    Dim r As Double         ' r
    r = t ^ 2 / (t ^ 2 + df)
    Dim p As Double         ' p(i)
    Dim ptot As Double      ' total of the p(i)
    ptot = 1#
    Dim ctot As Double      ' total of the cdf values
    ctot = 0#
    Dim hi As Long          ' upper index
    Dim low As Long         ' lower index
    low = Int(lambda)
    hi = low + 1
    On Error GoTo ErrorHandler
    Dim startTime As Double
    startTime = Timer
    '-- 输入参数记录 --
    Debug.Print "===== 函数调用开始 ====="
    Debug.Print "t: " & Format(t, "0.0000"), "df: " & df, _
        "d: " & Format(d, "0.0000"), "iter: " & iter
    Dim isWPS As Boolean
    isWPS = InStr(Replace(UCase(Application.Caption), UCase(Application.ActiveWorkbook.Name), ""), "WPS") > 0
    If isWPS Then
        With Application.WorksheetFunction
            For i = 1 To iter
                If low >= 0 Then
                    p = .Poisson(low, lambda, False)
                    ptot = ptot - p
                    ctot = ctot + p * .BetaDist(r, low + 0.5, df2) + _
                        dsq * p * Exp(.GammaLn(low + 1) - .GammaLn(low + 1.5)) * .BetaDist(r, low + 1#, df2)
                    low = low - 1
                End If
               
                p = .Poisson(hi, lambda, False)
                
                 Debug.Print "hi", r, hi, lambda
                ptot = ptot - p
                 Debug.Print "ptot", r, hi, df2
                ctot = ctot + p * .BetaDist(r, hi + 0.5, df2) + _
                    dsq * p * Exp(.GammaLn(hi + 1) - .GammaLn(hi + 1.5)) * .BetaDist(r, hi + 1#, df2)
               
                If ptot < prec Then Exit For
                hi = hi + 1
                Debug.Print "迭代 " & i & ": low=" & low & " hi=" & hi & _
                    " ptot=" & Format(ptot, "0.000E+00") & _
                    " ctot=" & Format(ctot, "0.0000")
            Next
            NT_DIST_CDF = ctot / 2 + .Norm_S_Dist(-d, True)
        End With
    Else
        With Application.WorksheetFunction
       
            For i = 1 To iter
                If low >= 0 Then
                    p = .Poisson_Dist(low, lambda, False)
                    ptot = ptot - p
                    ctot = ctot + p * .Beta_Dist(r, low + 0.5, df2, True) + _
                        dsq * p * Exp(.GammaLn(low + 1) - .GammaLn(low + 1.5)) * .Beta_Dist(r, low + 1#, df2, True)
                    low = low - 1
                End If
                p = .Poisson_Dist(hi, lambda, False)
                ptot = ptot - p
                ctot = ctot + p * .Beta_Dist(r, hi + 0.5, df2, True) + _
                    dsq * p * Exp(.GammaLn(hi + 1) - .GammaLn(hi + 1.5)) * .Beta_Dist(r, hi + 1#, df2, True)
                If ptot < prec Then Exit For
                hi = hi + 1
                Debug.Print "迭代 " & i & ": low=" & low & " hi=" & hi & _
                    " ptot=" & Format(ptot, "0.000E+00") & _
                    " ctot=" & Format(ctot, "0.0000")
            Next
            NT_DIST_CDF = ctot / 2 + .Norm_S_Dist(-d, True)
        End With
    End If
    Debug.Print "计算耗时: " & Format(Timer - startTime, "0.000") & "秒"
    Debug.Print "最终结果: " & Format(NT_DIST_CDF, "0.0000E+00")
ErrorHandler:
    Debug.Print "错误 " & Err.Number & ": " & Err.Description
    Debug.Print "出错时变量状态: low=" & low & " hi=" & hi
End Function