/***********************************************************
p_chi2.c -- カイ2乗分布
***********************************************************/
Const PI = 3.14159265358979323846264
Function p_nor(z As Double) As Double /* 正規分布の下側累積確率 */
Dim i As Long
Dim z2, prev, p, t
z2 = z * z
p = z * Exp(-0.5 * z2) / Sqr(2 * PI):t = p
i = 3
While i < 200
prev = p: t = t * (z2 / i): p = p + t
If (p = prev) Then
p_nor = 0.5 + p
Exit Function
End If
i = i + 2
Wend
If (z > 0) Then p_nor=1 Else p_nor=0
End Function
Function q_nor(z As Double) As Double /* 正規分布の上側累積確率 */
q_nor = 1 - p_nor(z)
End Function
Function q_chi2(df As Long, chi2 As Double) As Double /* 上側累積確率 */
Dim k As Long
Dim s, t, chi
if (df And 1) Then /* 自由度が奇数 */
chi = Sqr(chi2)
if (df = 1) Then
q_chi2 = 2 * q_nor(chi)
Exit Function
End If
t = chi * Exp(-0.5 * chi2) / Sqr(2 * PI): s = t
k = 3
While k < df
t = t * (chi2 / k): s = s + t
k = k + 2
Wend
q_chi2 = 2 * (q_nor(chi) + s)
Else /* 自由度が偶数 */
t = Exp(-0.5 * chi2): s = t
k = 2
While k < df
t= t * (chi2 / k): s = s + t
k=k+2
Wend
q_chi2= s
End If
End Function
Function p_chi2(df As Long, chi2 As Double) As Double /* 下側累積確率 */
p_chi2 = 1 - q_chi2(df, chi2)
End Function
#N88BASIC
Sub main()
Dim i As Long
Dim chi2
Print "***** p_chi2(df, chi2) *****"
Print "chi2 df=1 df=2 df=5 df=20"
For i=0 To 19
chi2 = 0.5 * i
Print chi2, p_chi2(1, chi2), p_chi2(2, chi2), _
p_chi2(5, chi2), p_chi2(20, chi2)
Next
End Sub
main()
