#/***********************************************************
#    p_chi2.rb -- カイ2乗分布
#***********************************************************/

PI = 3.14159265358979323846264

def p_nor( z)  # 正規分布の下側累積確率 
    z2 = z * z
    t = p = z * Math::exp(-0.5 * z2) / Math::sqrt(2 * PI)
    3.step(200, 2) do |i| 
        prev = p;  t *= z2 / i;  p += t
        if (p == prev); return 0.5 + p; end
    end
    return (z > 0) ? 1 : 0
end

def q_nor( z)  # 正規分布の上側累積確率 
    return 1 - p_nor(z)
end

def q_chi2( df,  chi2)  # 上側累積確率 
    if (df & 1) == 1  # 自由度が奇数 
        chi = Math::sqrt(chi2)
        if (df == 1); return 2 * q_nor(chi); end
        s = t = chi * Math::exp(-0.5 * chi2) / Math::sqrt(2 * PI)
        3.step(df - 1, 2) do |k|
            t *= chi2 / k;  s += t
        end
        return 2 * (q_nor(chi) + s)
    elsif       # 自由度が偶数 
        s = t = Math::exp(-0.5 * chi2)
        2.step(df - 1, 2) do |k|
            t *= chi2 / k;  s += t
        end
        return s
    end
end

def p_chi2( df,  chi2)  # 下側累積確率 
    return 1 - q_chi2(df, chi2)
end

printf("***** p_chi2(df, chi2) *****\n")
printf("chi2   %-16s %-16s %-16s %-16s\n",\
    "df=1", "df=2", "df=5", "df=20")
for i in 0...20
    chi2 = 0.5 * i
    printf("%4.1f %16.14f %16.14f %16.14f %16.14f\n",\
        chi2, p_chi2(1, chi2), p_chi2(2, chi2),\
              p_chi2(5, chi2), p_chi2(20, chi2))
end

exit 0